text:
Merge multiple channels each emits multiple files in nextflow
I need to merge multiple channels each emits multiple file. So far I did something as below:
process prephasing {
publishDir params.out, mode:'copy'
input:
tuple val(prefix), path(bfiles) from plink_data
each chrom from chroms
path hrc_sites_file
output:
tuple val("${chrom}"),path("chr${chrom}*.{bed,bim,fam}") into qc1_chr
"""
PRE_PASHING_QC2.sh "${hrc_sites_file}" "${prefix}" "${chrom}"
"""
}
channel.fromPath('/some/path/genetic_map_chr*_combined_b37.txt').into {shape_map;i
mp_map}
channel.fromPath('/some/path/*hrc.r1-1.ega.grch*.hap.gz').into {shape_ref;
imp_ref}
channel.fromPath('/some/path/*hrc.r1-1.ega.grch37.chr*.legend.gz').into {s
hape_legend;imp_legend}
channel.fromPath('/some/path/*hrc.r1-1.ega.grch37.chr*.samples').into {sha
pe_sample;imp_sample}
qc1_chr.filter{name, bfiles -> bfiles.collect {it.getBaseName()}.every {it.endsWith("step10")}}.into {shape_in;imp_in}

shape_para=shape_in.merge(shape_map).merge(shape_ref).merge(shape_legend).merge(shape_sample)

process phasing {
publishDir params.out, mode:'copy'
input:
each file(penta) from shape_para
output:
path("*step10*") into phased
script:
def (in, map, ref, legend, sample) = penta
"""
shapeit2 -T 32 --input-bed "${in}" --input-map "${map}" --input-ref "${ref}" "${legend}" "${sample}" --duohmm -O "${in}.phased" --seed 54321
"""
}

The "prephasing" process works fine and produce expected output. But the "phasing" process does not produce any output nor any error. Any help?
I have tried according to @Steve as below:
channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/1000Genomes/genetic_map_chr*_combined_b37.txt',size:1).map {group_key, file_list -> tuple(group_key.replaceFirst(/^genetic_map_/,""),file_list.first())}.set {genetic_map}

channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/haplegendsample/*hrc.r1-1.ega.grch37.chr*.haplotypes.{hap.gz,legend.gz,samples}',size:3).map{group_key,file_list -> tuple(group_key.replaceFirst(/^*HRC\.r1-1\.ega\.grch37\./,""),file_list)}.set {ref_panel}

qc1_chr.filter{chrom, bfiles -> bfiles.collect {it.getBaseName()}.every {it.endsWith("step10")}}.join(ref_panel).join(genetic_map).into {shape_in;imp_in}

process phasing {
publishDir params.out, mode:'copy'
input:
tuple val(chrom), path(input_bed), path(intput_ref), path(input_map) from shape_in
output:
tuple val("${chrom}"),path("${chrom}.{phased.haps,phased.sample,log}") into phased
script:
def (bed, bim, fam) = input_bed
def (haplotype, legend, sample) = input_ref
"""
shapeit2 --thread "${task.cpus}" --input-bed "${bed}" "${bim}" "${fam}" --input-map "${input_map}" --input-ref "${haplotype}" "${legend}" "${sample}" --duohmm --output-max "${chrom}.phased" --seed 54321
"""
}

But getting no output so far, (attached image)
I have tried to check individual channel separately using channel.view(), identified the problems, modified the code and checked again. Using the code below I can see individual channel emits 22 lines [chr1-22,[file_list]]. After merging I can see [chr1-22,[joined_file_list].
channel
    .fromFilePairs( "output3/*.{bed,bim,fam}", size:3 )
    .set { plink_data }
plink_data.filter{chrom, bfiles -> bfiles.collect {it.getBaseName()}.every {it.endsWith("step10")}}.map{group_key,file_list -> tuple (group_key.
replaceFirst(/^*\.step10/,""), file_list)}.set {in_pl}

channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/1000Genomes/genetic_map_chr*_combined_b37.txt',size:1).map {
group_key, file_list -> tuple(group_key.replaceFirst(/^*genetic_map_/,""),file_list.first())}.set {genetic_map}
//genetic_map.view()

channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/haplegendsample/*hrc.r1-1.ega.grch37.chr*.haplotypes.{ha
p.gz,legend.gz,samples}',size:3).map {group_key,file_list -> tuple(group_key.replaceFirst(/^.*.hrc\.r1-1\.ega\.grch37\./,""),file_list)}.set {re
f_panel}
//ref_panel.view()

in_pl.join(ref_panel).join(genetic_map).into {shape_in;imp_in}
process phasing {
//publishDir params.out, mode:'copy'
input:
tuple val(chrom), path(input_bed), path(intput_ref), path(input_map) from shape_in
output:
tuple val("${chrom}"),path("${chrom}.{phased.haps,phased.sample,log}") in
to phased
script:
def (bed, bim, fam) = input_bed
def (haplotype, legend, sample) = input_ref
"""
//shapeit="/cluster/projects/p697/projects/moba_qc_imputation/software/shap
eit"
shapeit --thread "${task.cpus}" --input-bed "${bed}" "${bim}" "${fam}" --input-map "${input_map}" --input-ref "${haplotype}" "${legend}" "${samp
le}" --duohmm --output-max "${chrom}.phased" --seed 54321
"""
}

But now I get error:
No such variable: chrom

But I decalred "chrom" in the input
The input bed|bim|fam files from channel from FilePairs works but from a output channel of previous process does not. If I run phasing solo and all the inputs from channel from FilePairs it works. But if I run prephasing and phasing together, the prephasing process works but for phasing process not output and no error; Not sure what I am doing wrong.
hrc_sites_file = file(params.hrc_sites_file)

Channel.fromFilePairs( "${params.input_dir}/*.{bed,bim,fam}", size:3 ).set { plink_data }

Channel.of(1..22).set { chroms }
params.out='output/'

process prephasing {
publishDir params.out, mode:'copy'
input:
tuple val(prefix), path(bfiles) from plink_data
each chrom from chroms
path hrc_sites_file
output:
tuple val("${chrom}"),path("chr${chrom}*.{bed,bim,fam}") into qc1_chr
tuple val("${chrom}"),path("chr${chrom}*log") into qc1_log
tuple val("${chrom}"),path("*-chr${chrom}.step4-HRC.txt") into qc1_hrc
path("chr${chrom}.step1b.snplist.txt") into qc1_snps
path("chr${chrom}.step1a.dupvar") into qc1_dups
"""
PRE_PASHING_QC2.sh "${hrc_sites_file}" "${prefix}" "${chrom}"
"""
}
channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/1000Genomes/genetic_map_chr*_combined_b37.txt',size:1).map {group_key, file_list -> tuple(group_key.replaceFirst(/^*gen
etic_map_/,""),file_list.first())}.set {genetic_map}

channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/haplegendsample/*hrc.r1-1.ega.grch37.chr*.haplotypes.{hap.gz,legend.gz,samples}',size:3).map{group_key,file_list ->
 tuple(group_key.replaceFirst(/^.*.hrc\.r1-1\.ega\.grch37\./,""),file_list)}.set {ref_panel}

qc1_chr.filter{chrom, bfiles -> bfiles.collect {it.getBaseName()}.every {it.endsWith("step10")}}.map{group_key,file_list -> tuple (group_key.replaceFirst(/^*\.step10/,""), file_list)}.join(ref_panel).joi
n(genetic_map).into {shape_in;imp_in}
imp_in.view()
process phasing {
publishDir params.out, mode:'copy'
input:
tuple val(chrom), path(input_bed), path(input_ref), path(input_map) from shape_in
output:
tuple val("${chrom}"),path("${chrom}.*") into phased
script:
def (bed, bim, fam) = input_bed
def (haplotype, legend, sample) = input_ref
"""
shapeit --thread "${task.cpus}" --input-bed "${bed}" "${bim}" "${fam}" --input-map "${input_map}" --input-ref "${haplotype}" "${legend}" "${sample}" --duohmm --output-max "${chrom}.step10" --seed 54321
"""
}

The merge operator is deprecated and will be removed from future Nextflow releases. You'll need instead the join operator to join by chromosome. You can use the fromFilePairs operator to quickly get the group key, formatting it as required:
channel
    .fromFilePairs("/some/path/HRC.r1-1.EGA.GRCh37.chr*.{hap.gz,legend.gz,samples}", size: 3)
    .map { group_key, file_list ->
        tuple( group_key.replaceFirst(/^HRC\.r1-1\.EGA\.GRCh37\./, ""), file_list )
    }
    .set { ref_panel_files }

channel
    .fromFilePairs("/some/path/genetic_map_chr*_combined_b37.txt", size: 1)
    .map { group_key, file_list ->
        tuple( group_key.replaceFirst(/^genetic_map_/, ""), file_list.first() )
    }
    .set { genetic_map_files }

qc1_chr
    .map { chrom, bfiles ->
        tuple( chrom, bfiles.findAll { it.baseName.endsWith("step10") } )
    }
    .join( ref_panel_files )
    .join( genetic_map_files )
    .set { phasing_inputs }

process phasing {

    input:
    tuple val(chrom), path(input_bed), path(input_ref), path(input_map) from phasing_inputs

    script:
    def (bed, bim, fam) = input_bed
    def (haplotypes, legend, sample) = input_ref

    """
    shapeit2 \\
        --thread "${task.cpus}" \\
        --input-bed "${bed}" "${bim}" "${fam}" \\
        --input-ref "${haplotypes}" "${legend}" "${sample}" \\
        --input-map "${input_map}" \\
        --duohmm \\
        --output-max "${chrom}.phased" \\
        --seed 54321
    """
}

Thanks @Steve. Only one thing, the prefixes of hap.gz, legend.gz and samples are not same for each chromosome. That's why in my code there is a * at the beginning and another * at the chromosome label.

Hi @Steve I edited my question adding picture of running script and code snippet.

@zillurrahman: I think, though, you understand the general gist of the answer - that is, to create a tuple with a key as the first element which can then be used to join (multiple) channels. The above is just an example that uses the chromosome name as the join key. You should not expect a drop-in replacement for your code. I don't know what files your glob patterns match, and still don't. Your patterns could match multiple files, but I assumed you had one (hap, legend, sample) per chromosome. If the reference panel files have different prefixes, you may prefer separate channels for these.

Sorry for the inconveniences. I rechecked the code, there were some type. Modified them and can see joined items. But getting a new error. No such variable: chrom But I declared the chrom variable in the input.

The new code snippet was added to my above question.

@zillurrahman I can't reproduce the error you get using your new code snippet. It should have failed earlier because of this typo: path(intput_ref) should be path(input_ref). Do you know what line the error is being produced on?

Thanks @Steve. The error is occurring at line 20 where the output is declared.

In the .nextflow.log file, it is showing that the error is occurring at line 20. And line 20 is tuple val("${chrom}"),path("${chrom}.{phased.haps,phased.sample,log}") into phased where the output was declared.

@Steve you are right. There was a type in my code. The error is Unknown variable 'input_ref. The typo is fixed and the code is running now! Thanks a lot!

If I run phasing separately, it works. But in combination with prephasing  it does not work. I mean bed|bim|fam files from channel from FilePairs works but from the output channel of previous process does not. Even I can not see the joined channel.view() in .command.out in the work directory. No error also. not sure what goes wrong. I added the updated code snippet in the above question.

@zillurrahman: I think your ouput declaration for qc1_chr needs to be: tuple val("chr${chrom}"),path("chr${chrom}*.{bed,bim,fam}") into qc1_chr. HTH.

@Steve Thanks. I updated the output declaration of process prephasing. Still the process phasing does not work nor any error message. In the last screenshot of the above question, you can see that there is a work directory (at the left) assigned for prephasing process but nothing for phasing process. Not sure what I am doing wrong. The same phasing process works if I run it in a separate script where input bed|bim|fam  are assigned from FilePairs.

@zillurrahman: I suspect one of your input channels to join on isn't working as expected. Could you please post the output of qc1_chr.view(), ref_panel.view() and  genetic_map.view()? Thanks!

Thanks @Steve. ref_panel.view() gives [chr,[respective hap.gz,legend.gz,sample.gz]], genetic_map.view() gives [chr,[respective maps]]. The problem with filtering qc1_chr. qc1_chr.filter{chrom, bfiles -> bfiles.collect {it.getBaseName()}.every {it.endsWith("step10")}}.view()  shows nothing. Interestingly, same filtering method works when the bed|bim|fam files were from FilePairs. I also tried qc1_chr.filter { it[1].name.contains('step10') }.view() but nothing. qc1_chr.view() gives [chr1-22,[respective all plink files]]

My mistake, I think you'll need to replace the filter statement with .map { chrom, bfiles -> tuple( chrom, bfiles.findAll { it.baseName.endsWith("step10") } ) }. No filtering is required, we just need to select out the required step10 files.

Thank you very much @Steve. It works now. as there are multiple patterns, I used channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/haplegendsample/*hrc.r1-1.ega.grch37.chr*.{hap.gz,legend.gz,samples}',size:3).map {group_key,file_list -> tuple(gro up_key.replaceAll(/^.*.hrc\.r1-1\.ega\.grch37\.|\.haplotypes\.noibd|\.haplotypes/,""),file_list)}.set {ref_panel}. Do you think replaceAll is ok instead of replaceFirst . I got expected output using using replaceAll. Could I use multiple pattern using replaceFirst?

@zillurrahman: I don't think you can since the first looks like a prefix and the last two look like suffixes. I think you'd need two replaceFirst calls. I'd be tempted to just use the first string that looks like a chromosome name: group_key.find(~/chr\d{1,2}/). HTH.

Thanks @Steve I just tested the find option and it works. Thanks a lot! It seems fewer code and easier to write!

Hi @Steve, Wondering if you get a chance to look at https://bioinformatics.stackexchange.com/questions/15698/parse-csv-text-file-to-create-arrays-in-nextflow . I just need a little help to create an array from a csv file in nextflow.


text:
How to convert txt.gz into fastq.gz?
I'm a newbie to bioinfo. I'm trying to use galaxy.networkanalyst.ca for the analysis of RNA-Seq data and for that I need to convert txt.gz files into fastq.gz files.

Link: GSM4273445 is one sample and the whole series is from here.
I've Googled to find out if there are online tools to convert it, but there aren't any. Besides, there is no discussion addressing this format conversion on researchgate.
I intend to use this RNA-Seq data and investigate the relative expression of some specific transcripts of interest. I have been introduced to 'galaxy.networkanalyst.ca' where the very first step is to upload fastq files. In addition to this I was aware of ncbi GEOdatasets as source of open RnA-Seq data.
Can anyone suggest a method or tool for this?

I've had a quick look at the file you linked to and this is a featureCounts file - this list the number of reads for each gene/transcript (COUTNS.txt in the file name pretty much gave it away that it was going to be a file of, well counts).
I can't help you with Galaxy, but if you want to do an RNA-seq analysis, the count data + metadata is likely all you'll need.
If you really want to do the preprocessing and mapping yourself, the fastq files can be downloaded from the SRA archive (SRP242087) - click on the SRA link in your whole series link.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Why is there no GT-entry in my .vcf file from bacteria
I mapped raw reads obtained from an E. coli genome (.fastq file) to my reference genome. Next, I exported the .vcf (variant calling file), which looks something like the header example shown below (total >50 lines):
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Variants:_LGC19-XL01_S63_L001_R_001_(trimmed)
CFC381_K12_Bw25113  360287  .   GGGAAT  AATTCC  1952.76 .   NS=1;VF=0.380;SB=0.533;SB50=0.60;SB65=1.0;TYPE=Substitution;AVQUAL=30   DP:AO   231:89

How can I back-fill or quickly assign the missing GTs to the .vcf ?
For my downstream process (plink --pca, via bfctools merge), I need that there are GT entries in every .vcf, so I can merge them and the respective GT is shown in the merged .vcf. I only found solutions that need the .bam file, but that is laborious and I don't think I need the .bam in the case of haploid data.

This method needs .bam: https://www.biostars.org/p/319304/
This method needs .fastq: https://training.galaxyproject.org/training-material/topics/variant-analysis/tutorials/non-dip/tutorial.html#examining-the-results

However, the.fastq is actually from a mixed culture, so I wonder do I really need to follow the .bam-strategy? Suggestion for a tool? But in that case, why did it not export the GTs in the fist place (I exported from "geneious" sequence analysis program)

The representation for a variant present in a haploid genome should be 1. A heterozygote in a diploid genome would be 0/1 (or 0|1 or 1|0  if phased). So, to add a 1 as the genotype, you would need to first insert a GT to the beginning of each INFO line and then the 1 at the beginning of the last field. But, this should happen only on lines that are not headers (do not start with #). This should do it:
$ awk -F'\t' -v OFS="\t" '{ if(/^#/){print}else{$9="GT:"$9;$10="1:"$10;print}}' file.vcf
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Variants:_LGC19-XL01_S63_L001_R_001_(trimmed)
CFC381_K12_Bw25113  360287  .   GGGAAT  AATTCC  1952.76 .   NS=1;VF=0.380;SB=0.533;SB50=0.60;SB65=1.0;TYPE=Substitution;AVQUAL=30   GT:DP:AO    1:231:89

So just redirect the output to a fixed file:
$ awk -F'\t' -v OFS="\t" '{ if(/^#/){print}else{$9="GT:"$9;$10="1:"$10;print}}' file.vcf > fixed.vcf

Now, in this specific case, you could also do something like this:
$ sed -i.bak '/^[^#]/{s/\<DP:AO\t/GT:DP:AO\t1:/}' file.vcf

But that assumes that all lines will have the exact same INFO field of DP:AO which is not a safe assumption since each variant can have different INFO fields.

Perfect! It worked super good, I am so grateful to both of you!
So one thing has gotten forgotten, one also needs to adjust the header (it's easy).
Check my post here: https://www.biostars.org/p/9463791/#9467359

You're very welcome, @DavidStreuli. If one of the answers here solved your issue, please take a moment and accept it by clicking on the checkmark on the left. That is the best way to express your thanks on the Stack Exchange sites. Feel very fre to wait for a better answer, no need to accept mine now, I just wanted to let you know how the site works.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Error downloading GSE dataset
I've installed GEOquery with
install.packages("BiocManager")
BiocManager::install("GEOquery")
library(GEOquery)

then I ran the command
GSE68808 <- getGEO("GSE68808")

but it resulted in..
Found 1 file(s)
GSE68808_series_matrix.txt.gz
trying URL 'https://ftp.ncbi.nlm.nih.gov/geo/series/GSE68nnn/GSE68808/matrix/GSE68808_series_matrix.txt.gz'
Content type 'application/x-gzip' length 4138819 bytes (3.9 MB)
downloaded 3.9 MB

error: 'C:\Users\(myname)\AppData\Local\Temp\RtmpUp5Sv7\file8e844a7a59bc' does not exist.

ㅜㅜ
How can I solve this? thanks.

Please consider accepting the answer if it has answered your question

There seems to be some problem with your "platform specific tempdir" that is used by getGEO's destdir parameter. See ?getGEO for more info on the destdir parameter.
Try
GSE68808 <- getGEO("GSE68808", destdir = getwd())


text:
Does SBOL support timing and threshold value parameters?
Like the propagation delay and threshold value parameters in electronic circuits, does SBOL format support any such parameters/glyph to represent the propagation delay (time to trigger the output protein once the input is provided) and threshold value (minimum input concentration required to trigger the model's output) of genetic models?
If currently not supported, can these parameters be incorporated in the SBOL visual representation which would give reader an idea (just by looking at the schematic diagram) about the timing behavior of a model?

Note: I will use SBOL3 in my descriptions; the answer can be equivalently translated to SBOL2
At present, SBOL doesn't have any explicit recommendations for how to represent these types of timing parameters. This is because there is not yet a clear scientific consensus on whether or not such parameters are a good way to represent genetic circuits.
That said, one of the advantages of using SBOL is that its extensible RDF-based model allows additional information like this to be added into a model. Two good ways of doing so are:

If the parameters are simple and can stand alone as measurements, then they can be added to designs as a Measure object with appropriate units and type. For example, if you have biological devices with a 90 minute propagation delay, then on the Component representing the device you can add a Measure with numerical value 90, units of om:minute.  You'd likely want two types, one of SBO:0000225 (Delay), and a custom type for your own notion of propagation delay, since I don't think that's included in any ontology yet.
If the parameters are more complex and/or inter-related, e.g., a full data-sheet, then you'd likely want to instead create one or more extension classes that are used to represent the information. A library that makes this easy is sbol_factory, though one can do it with any of the SBOL libraries.

This becomes more complicated if you want to talk about behavior in different environmental contexts, but the basic principles are still the same.
Finally, as for SBOL visual: again, there's no explicit representation, but the model is extensible, and this information can be added as any sort of text or custom glyph as long as it doesn't conflict with ones already in the standard.


text:
How to trigger next process on completion of previous process?
Simple question but I am stumped. (I am still learning Nextflow)
I have a process which takes as input a channel of files and processes each file separately.
There is a second process which needs to wait until all files are complete, so it can merge the outputs into a single collated file.
How can I do this? I have tried subscribing to "onComplete" event to no avail as I cannot specify a process as the function to trigger.

files = Channel.from params.the_files_list

process process_files {

    input:
    val file_path from files

    output:
    stdout result

    """
    python ./process_one_file.py $file_path
    """
}

process merge_results {

    input:
    ????

    output:
    ????

    script:
    """
    echo "How to trigger this process once, after all individual files are processed?"
    """

}

??? Want to do something like this
result.subscribe onComplete: merge_results

```

I've run into the same problem but I'm interested in trying this method. Did it work for you?

So from this page...
https://groups.google.com/g/nextflow/c/jt77_-uApMs
It is the "scatter-gather" pattern. I think the answer is to collect() the channel as follows - interested to hear if there is a better way
files = Channel.from params.the_files_list

process process_files {

    input:
    val file_path from files

    output:
    stdout result

    """
    python ./process_one_file.py $file_path
    """
}

process merge_results {

    input:
    val go from result.collect()

    script:
    """
    echo "Triggered once after all files complete!"
    """

}

.collect() is the correct way to do this


text:
parse csv/text file to create a channel with array in nextflow
[![example csv file][1]][1]
First I want to start by parsing and doing some operation on a text/csv file. How could I read a csv/text file in nextflow?
Suppose, I have above csv file with 3 columns. Lets say I want to create a channel which emits an array from minimum of column 2 (7.9e5) to maximum of column 2 (8.8e5) with an increment of 1e4. How could I create the following channel and later access it in the script? Any help? The
array=([1]="7.9e5 8.0e5" [2]="8.0e5 8.1e5" [3]="8.1e5 8.2e5" [4]="8.2e5 8.3e5" [5]="8.3e5 8.4e5" [6]="8.4e5 8.5e5" [7]="8.5e5 8.6e5" [8]="8.6e5 8.7e5" [9]="8.7e5 8.8e5" [10]="8.8e5 8.9e5")

I need to do some operations with this array in the script which will generate multiple files for each element in the array. These multiple files need to be merged in a single file.
Sorry, I could not provide any nextflow script that I tried, because I not sure how create an array with splitCsv
I have tried in python:
df1=pd.read_csv('data1/omni-ec-eur-fin-pass-qc.bim',sep='\t',header=None)
arr1=np.arange(df1.iloc[:,3].min()-1,df1.iloc[:,3].max()+1,3000000)
diffs = []
for x in zip(arr1[0::],arr1[1::]):
    diffs.append(x)
arr2=np.arange(1,len(diffs)+1,1)
dict(zip(arr2,diffs))

Which gave me:
{1: (12343, 3012343),
 2: (3012343, 6012343),
 3: (6012343, 9012343),
 4: (9012343, 12012343),...

But I need something like this:
([1]="12343 3012343" [2]="3012343 6012343" [3]="6012343   9012343" [4]="9012343 12012343"....)

It could be arrays (22 arrays for 22 chromosome) or a csv file containing arrays for each chromosome which could be used as an input channel. There are 22 bim files (one for each chromosome) stored in the output channel which would be used to create the above mentioned array.
The original data inside the bim looks like this:
    0   1   2   3   4   5
0   1   rs11240777  0   798959  A   G
1   1   rs4475691   0   846808  A   G
2   1   rs7537756   0   854250  G   A
3   1   rs13302982  0   861808  A   G
4   1   rs1110052   0   873558  C   A

The 4th column would be used to create the arrays.
[1]: https://i.stack.imgur.com/cxqyl.png

PLease [edit] your question and give us an example file and minimal nextflow script to work with.

Is there any book to learn nextflow operations? I was trying to find any book to learn nextflow but failed. At least if there is any example of creating array in nextflow would be great.

Have you read the nextflow documentation? I don't know much about it, but the whole idea is to use it for control flow; it runs other programs and they can have arrays. But next flow isn't a programing language per se. You write your script in Python or peel or whatever and then run the pipeline if scripts with nextflow.

@zillurrahman Do you really need an array? Or do you mean a channel? Your last paragraph suggests the latter would be a better fit. Please edit your question to clarify exactly what is needed.

@zillurrahman Generally, it's better to off-load all work to one or more processes. Writing a process to create a range of values, like in your example, is trivial in Perl/Python. You can use your favorite language in a process' script block. Please see: https://www.nextflow.io/docs/latest/process.html#scripts-a-la-carte

@terdon It's an extension of the Groovy programming language :-)

@steve sure, but it isn't designed to be used as a regular programing language but as a process control system to build and manage pipelines. It seems far more reasonable to manage the data in another script and then have nextflow manage that script instead of trying to use nextflow as your scripting language directly. Which is exacy what you suggested too.

Thanks @Steve and @terdon, I edited the above question. Hope, it would clarify my necessity.  However, I am also trying to create such array in bash . But, I think probably a nextflow channel would be more appropriate in this case.

@terdon wasn't disagreeing with you on that...

@zillurrahman Is there some reason for choosing Bash over Python for example? Bash isn't a good choice given your requirements.

@Steve No reason to choose bash. I tried in python and got an array. np.arange(df1.iloc[:,3].min()-1,df1.iloc[:,3].max()+1,3000000) which gives me: array([    12343,   3012343,   6012343,   9012343,  12012343, .... . How could I get something like: ([1]="12343 3012343" [2]="3012343 6012343" [3]="6012343   9012343" [4]="9012343 12012343"....)

So far, I have done: diffs = [] for x in zip(arr1[0::],arr1[1::]):     diffs.append(x) diffs=np.array(diffs) arr2=np.arange(1,len(diffs)+1,1) arr3=np.column_stack((arr2,diffs)) arr4=arr3.tolist() which gave me [[1, 12343, 3012343],  [2, 3012343, 6012343],  [3, 6012343, 9012343], . Should I convert it to string and use regex to get ([1]="12343 3012343" [2]="3012343 6012343" [3]="6012343   9012343" [4]="9012343 12012343"....) or there are better alternatives?

@zillurrahman Do you really need to create the array string? If all you eventually want to do is create a channel of tuples, your first process should write a simple CSV which can then be split using the splitCsv operator.

@zillurrahman Also, your comments above should be edits to your question. Small snippets of unformatted code isn't especially helpful. Please include a small minimal example and the code that reproduces the issue.

Thanks @Steve. I have edited my above question. Hope it makes sense now.

@zillurrahman Like terdon mentioned earlier, please include some example input (i.e. cat data1/omni-ec-eur-fin-pass-qc.bim). Please format it in a code block.

@Steve . The first few rows of the input data is added to the above question.

Hi, any help/suggestion for the above problem? Did I clarify enough the issue?


text:
Manta SV caller: I can't find a file in the path specified by the user guide
I trying to locally reproduce a Manta run as seen in the log of an illumina basespace RNA-Seq Alignment tool v1.1.1 run I've done. The log reads as follows:
[LOG][2021-03-12 20:58:56.377]03/12/2021,20:58:56.193,Launch process:  \
/path_to_RNA_Tools/bin/python  \
'/path_to_RNA_Tools/Manta/bin/configManta.py'  \
--rna --normalBam='path_to_alignments/alignments.bam' \
--runDir='path_to_samples/CTR/replicates/CTR/MantaFusions/manta.run'  \
--referenceFasta='/path_to_genomes/genome.fa' \
--config '/opt/illumina/Isis/2.6.25.18/Workflows/WholeGenomeRnaSeqWorker/scripts/configManta.py.ini'

The file shown in the line above configManta.py.ini is however no were to be found.
The Manta user guide states in the linked section

Options listed in the file: ${MANTA_INSTALL_PATH}/bin/configManta.py.ini

The configManta.py is in the expected directory, only the configManta.py.ini is absent.
I am using Manta version  1.6.0 from bioconda.
Your insights will be much appreciated!

The --config file option for Manta can be used to override the defaults in the global config file. If you type configManta.py -h, you'll see where the global config file has been installed. For me, it's here:
/home/steve/miniconda3/envs/manta/share/manta-1.6.0-1/bin/configManta.py.ini

So the command in your log file takes a specific config file installed into:
/opt/illumina/Isis/2.6.25.18/...

I suspect 'Isis' is a reference to Illumina's analysis tool 'IsisRNA'. Papers that mention Illumina's 'RnaReadCounter' tool, say that it is part of the internal analysis tool 'IsisRNA'. Because it's one of Illumina's in-house tools, I think you'd need to contact Illumina support to get a copy. I would only do this though if you need to reproduce exactly as per the command in your log file. If you only need to roughly reproduce the output, the defaults in the global config file are probably fine.


text:
SnpEff not correctly annotating multiple-nucleotide polymorphisms
I have some SARS-CoV-2 sequencing data that I'm trying to annotate with SnpEff, however SnpEff doesn't appear to be recognizing multiple adjacent SNPs within the same codon and correctly calling them as a single MNP (multiple-nucleotide polymorphism).
I'm using the GenBank reference MN908947.3. Here are the relevant lines of the annotated VCF file
MN908947.3  28280   .   G   C   2911.06 PASS    AC=2;AF=1.00;AN=2;DP=65;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=31.56;SOR=8.380;ANN=C|missense_variant|MODERATE|N|Gene_28273_29532|transcript|QHD43423.2|protein_coding|1/1|c.7G>C|p.Asp3His|7/1260|7/1260|3/419||,C|upstream_gene_variant|MODIFIER|ORF10|Gene_29557_29673|transcript|QHI42199.1|protein_coding||c.-1278G>C|||||1278|,C|downstream_gene_variant|MODIFIER|S|Gene_21562_25383|transcript|QHD43416.1|protein_coding||c.*2896G>C|||||2896|,C|downstream_gene_variant|MODIFIER|ORF3a|Gene_25392_26219|transcript|QHD43417.1|protein_coding||c.*2060G>C|||||2060|,C|downstream_gene_variant|MODIFIER|E|Gene_26244_26471|transcript|QHD43418.1|protein_coding||c.*1808G>C|||||1808|,C|downstream_gene_variant|MODIFIER|M|Gene_26522_27190|transcript|QHD43419.1|protein_coding||c.*1089G>C|||||1089|,C|downstream_gene_variant|MODIFIER|ORF6|Gene_27201_27386|transcript|QHD43420.1|protein_coding||c.*893G>C|||||893|,C|downstream_gene_variant|MODIFIER|ORF7a|Gene_27393_27758|transcript|QHD43421.1|protein_coding||c.*521G>C|||||521|,C|downstream_gene_variant|MODIFIER|ORF8|Gene_27893_28258|transcript|QHD43422.1|protein_coding||c.*21G>C|||||21| GT:AD:DP:GQ:PL  1/1:0,65:65:99:2925,196,0
MN908947.3  28281   .   A   T   2912.06 PASS    AC=2;AF=1.00;AN=2;DP=78;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=26.91;SOR=8.380;ANN=T|missense_variant|MODERATE|N|Gene_28273_29532|transcript|QHD43423.2|protein_coding|1/1|c.8A>T|p.Asp3Val|8/1260|8/1260|3/419||,T|upstream_gene_variant|MODIFIER|ORF10|Gene_29557_29673|transcript|QHI42199.1|protein_coding||c.-1277A>T|||||1277|,T|downstream_gene_variant|MODIFIER|S|Gene_21562_25383|transcript|QHD43416.1|protein_coding||c.*2897A>T|||||2897|,T|downstream_gene_variant|MODIFIER|ORF3a|Gene_25392_26219|transcript|QHD43417.1|protein_coding||c.*2061A>T|||||2061|,T|downstream_gene_variant|MODIFIER|E|Gene_26244_26471|transcript|QHD43418.1|protein_coding||c.*1809A>T|||||1809|,T|downstream_gene_variant|MODIFIER|M|Gene_26522_27190|transcript|QHD43419.1|protein_coding||c.*1090A>T|||||1090|,T|downstream_gene_variant|MODIFIER|ORF6|Gene_27201_27386|transcript|QHD43420.1|protein_coding||c.*894A>T|||||894|,T|downstream_gene_variant|MODIFIER|ORF7a|Gene_27393_27758|transcript|QHD43421.1|protein_coding||c.*522A>T|||||522|,T|downstream_gene_variant|MODIFIER|ORF8|Gene_27893_28258|transcript|QHD43422.1|protein_coding||c.*22A>T|||||22| GT:AD:DP:GQ:PL  1/1:0,65:65:99:2926,196,0
MN908947.3  28282   .   T   A   2912.06 PASS    AC=2;AF=1.00;AN=2;DP=78;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=28.71;SOR=8.380;ANN=A|missense_variant|MODERATE|N|Gene_28273_29532|transcript|QHD43423.2|protein_coding|1/1|c.9T>A|p.Asp3Glu|9/1260|9/1260|3/419||,A|upstream_gene_variant|MODIFIER|ORF10|Gene_29557_29673|transcript|QHI42199.1|protein_coding||c.-1276T>A|||||1276|,A|downstream_gene_variant|MODIFIER|S|Gene_21562_25383|transcript|QHD43416.1|protein_coding||c.*2898T>A|||||2898|,A|downstream_gene_variant|MODIFIER|ORF3a|Gene_25392_26219|transcript|QHD43417.1|protein_coding||c.*2062T>A|||||2062|,A|downstream_gene_variant|MODIFIER|E|Gene_26244_26471|transcript|QHD43418.1|protein_coding||c.*1810T>A|||||1810|,A|downstream_gene_variant|MODIFIER|M|Gene_26522_27190|transcript|QHD43419.1|protein_coding||c.*1091T>A|||||1091|,A|downstream_gene_variant|MODIFIER|ORF6|Gene_27201_27386|transcript|QHD43420.1|protein_coding||c.*895T>A|||||895|,A|downstream_gene_variant|MODIFIER|ORF7a|Gene_27393_27758|transcript|QHD43421.1|protein_coding||c.*523T>A|||||523|,A|downstream_gene_variant|MODIFIER|ORF8|Gene_27893_28258|transcript|QHD43422.1|protein_coding||c.*23T>A|||||23| GT:AD:DP:GQ:PL  1/1:0,65:65:99:2926,196,0

These are the positions corresponding to codon 3 on the N gene, they should be annotated as a single MNP, p.Asp3Leu, but istead they've been annotated seperately, as three individual SNPs.
I'm using the latest version of SnpEff (v5), and according to the manual SnpEff should have this functionality, what am I doing wrong?
You can see an example of how it's supposed to work here.

Are they expressed as a single variant in the input file? Does the manual say snpEFF can detect MNPs or only that it can handle them?

@terdon no, in the input file they're expressed as seperate SNPs. Based on the example in the link I provided, my take away from the SnpEff manual was that adjacent SNPs on the same codon would be grouped as an MNP, however there is little elaboration so I can't say for sure what is required for SnpEff to call adjacent SNPs as an MNP. I should note that since the variant caller does not have access to the ORFs, I don't think it would be able to call MNPs, so it has to be snpEff that does it.

The ORFs are not relevant here. In order to call a MNP, the variant caller needs to find all alt residues supported by the same read. Given that in your example, the reference allele has 0 support in all three, it seems very likely that the caller would find all three in the same read, but I don't see how snpEFF can be expected to change the list of variants you give it. If you join these yourself, snpEFF will likely deal with it, but I don't see anything in the manual that suggests it can merge them for you.

@terdon The ORFs are absolutely relevant here. The purpose of calling adjacent SNPs as a single MNP is so that the amino acid mutation is correctly annotated, thus you need to know where the breaks are between each codon, and cannot be deduced by the variant caller just from the ref. genome. Again, if you click on the link that I provided, and take a look at the example in the SnpEff manual, you will see how SnpEff annotates separate adjacent SNPs as a single MNP.

There is of course the problem of needing to determine phase (in the case diploid organisms), or even in the case of haploid organisms, that, as you say, SnpEff cannot in general determine from the VCF file whether adjacent SNPs always occur together on the same read. However I can find no info in the manual about how SnpEff handles these types of ambiguities, and my attempts at modifying the INFO field in the VCF file to indicate ploidy has not fixed the issue.

No, the ORFs are not relevant because the variant calling has no concept of genes, let alone translated products or codons. The calling is done using the reference genome and the reads. Of course the ORfs are relevant when you try to figure out the effect of your variants, but they do not affect and cannot affect whether something should be considered one variant or many. That is a function of the reads alone. All I see in the manual is an example of MNP, not a case where the input VCF had them separate and snpEFF joined them. I may be missing it, of course, but that's all I can find.

Let us continue this discussion in chat.


text:
How to plot a graph of Linkage disequilibrium with limited information
I have a snp information in an excel file and I used SNPstat to do some analysis. But I also need a plot of LD for these snps. But my file is not formatted as the most programs required (.ped, hapmap, Plink...). Here is a screen shot of file:
I have also chromosome numbers and rs information for this snps. Due to this data format problem I couldn't use Haploview or other programs. Can you suggest me a solution to create a LD graph?

https://academic.oup.com/bioinformatics/article/28/19/2537/288671

you need an LD matrix, which your table is not.  The LD matrix is a pairwise matrix.  See this page for an example: https://www.rdocumentation.org/packages/gaston/versions/1.4.9/topics/LD.plot


text:
Inflate operation failed: progress temporarily not possible, or in() / out() returned an error
I download a BAM file from ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam (171G), and I am trying to view the comparison result of the BAM file using the following command:
[yyn@localhost All]$ samtools flagstat NA12878_GRCh38.bam

The result shows that the file is truncated, the result as follows:
[E::bgzf_uncompress] Inflate operation failed: progress temporarily not possible, or in() / out() returned an error
[E::bgzf_read] Read block operation failed with error 1 after 0 of 4 bytes
[bam_flagstat_core] Truncated file? Continue anyway.
72639 + 0 in total (QC-passed reads + QC-failed reads)
896 + 0 secondary
0 + 0 supplementary
2529 + 0 duplicates
72639 + 0 mapped (100.00% : N/A)
71743 + 0 paired in sequencing
36865 + 0 read1
34878 + 0 read2
48600 + 0 properly paired (67.74% : N/A)
69949 + 0 with itself and mate mapped
1794 + 0 singletons (2.50% : N/A)
10233 + 0 with mate mapped to a different chr
2869 + 0 with mate mapped to a different chr (mapQ>=5)

The total number of reads is only 72639, Obviously, this is the wrong result. And I checked the file again with the following command
[yyn@localhost All]$ samtools quickcheck -qv NA12878_GRCh38.bam && echo 'all ok' || echo 'failed check'

The result is:
all ok

I checked the contents of the file again, and the result showed that there was no problem with the first 72620 lines in the file, but there was a problem with the first 72630 lines, so the file should be truncated in the 72620 to 72630 lines.
[yyn@localhost All]$ samtools view NA12878_GRCh38.bam | head -n 72610 >> head72610.sam

That's ok.
[yyn@localhost All]$ samtools view NA12878_GRCh38.bam | head -n 72620 >> head72620.sam
[E::bgzf_uncompress] Inflate operation failed: progress temporarily not possible, or in() / out() returned an error
[E::bgzf_read] Read block operation failed with error 1 after 0 of 4 bytes
[main_samview] truncated file.

samtools sort and samtools index are the same error results as above. Does anyone know what's going on? Thanks!

As you say, the file is truncated. Why are you investigating this instead of just downloading the entire file?

I downloaded it three times using lftp, all of which are like this, and I can't restore it with bamUtil. But downloading with Aspera is slower than lftp.

My best guess would be a corruption during download as:
samtools quickcheck -qvvvv ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam
verbosity set to 4
checking ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam
opened ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam
ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam is sequence data
ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam has 196 targets in header.
ftp://ftp-trace.ncbi.nlm.nih.gov:/giab/ftp/data/NA12878/10Xgenomics_ChromiumGenome_LongRanger2.0_06202016/NA12878_GRCh38.bam has good EOF block.

Plus I can successfully run idxstats from remote and fetch individual chromosomes without error via samtools view. Try redownloading before digging deeper.

Thanks for your advice, I will redownload it by using Aspera instead of lftp.


text:
Holistic enzyme activity determination with computation
I'm currently working on a program which will determine the activity difference between an original enzyme and a variant with just 1 or 2 variants. Of course, I'm not talking about the "kinetic activity", but some indications of a probable activity difference should be given. The system has to be sensitive enough to determine the probable activity difference between unthought-of variants.
From the information that I've got from my research, there is two option. I can use one of 3 quantitive enzyme reaction analysis tools (such as EVB, FEB and LRA) or their hybrids.
EVB => Empirical Valence Bond
FEB => Free Energy Perturbation
LRA => Linear Response Approximation
Source:
https://www.researchgate.net/post/Can-the-activity-of-proteins-be-studied-using-molecular-dynamic-simulations
Now, my questions are: - Does these application be used for such purpose? - For an accurate prediction of  enzyme activity difference, are these programs good enough or require developement? - Which one of them can give the best result? - And lastly, are QM/MM methods a must for an accurate answer?
I'm a high school student so sorry if my questions sound dumb.
If you can answer my inquiry's I would appreciate it:)
Thanks

Your title says holistic. This is a tad problematic as there's layers upon layers. Say, post-translation regulation, inhibiting metabolites, interacting protein etc. That is why when talking of an enzyme inhibitor, at the biochemical level one speaks in terms of k_i (inhibition constant), while at the cellular level ("holistic") one talks of IC50.
This could be a biocatalysis/synthetic biology problem (make an enzyme better) or translational medicine problem (why is this mutation pathogenic).
In terms of in silico enzyme design, one can have all sorts of issues that need to be addressed:

catalysis
side reactions
stability
aggregation
solubility

Also there's the expression problems, which may be due to subpar promoter, insufficient cofactors and poor codon usage.
For clinical variant assessment, a variant may:

destabilise leaving to degradation
remove a post-translational modification
disrupt the active site, affecting catalysis
lose a binding partner
result in incorrect localisation
result in solubility issues

So force-field scores are not the sole solution. Furthermore, a forcefield method gives one a ∆∆G (difference in relative Gibbs free energy), not a K_M (Michaelis constant) and a k_cat (catalytic rate), although it is somewhat possible to calculate these, or a relative variant of these.

a variant with just 1 or 2 variants

A variant protein can have one or two mutations (e.g. G10R) in its sequence. A mutant and a variant are near synonyms, but variant is generally intended to be a mutant in a population.

Does these application be used for such purpose?

For enzyme design MD simulations are not generally the first step. Rosetta pmut_scan for example is a nice starting step as the mutational landscape of a protein is N^19, while a MD simulation takes weeks and 6+ core per single variant.
This calculates only an increase in stability (relative difference in Gibbs folding energy), not reactivity. For that you can use a transition state, whose results are inaccurate, but majorly affected by the set-up, i.e. junk-in junk-out plays a big part. There's a saying 'a week in the lab saves you an hour at the library': this is backwards for sarcasm as it happens too often, namely that is background research is always a lifesaver. An QM-MD simulation on the wild type reaction definitely can help in elucidating things if it is not clear.
Online tools
Generally, however, there are resources out there to do many stability calculations already.
SynBio If you are interested in picking a residue or two to make an enzyme more stable, without resorting to a thermophilic enzyme, PROSS is a very nice tool.
Clinical (Caveat: If your variant is a variant that affects you, a friend or a relative, you/they should speak to a genetic councillor and not do research on it yourself) There are several tools that can describe a variant at the protein level. Do remember that ∆∆G is a potential so more negative is more stable —except in SDM where the score is inverted.

Missense3D
Venus
Miscast
SDM


text:
Error in FindIntegrationAnchors- Seurat package
I am working on integrating a labelled single cell RNA seq cell atlas with an unlabelled one. I am wondering how do I determine the value of max.features when integrating two single cell RNA seq datasets?
Below is the code I am using:
epithelial = Read10X(data.dir="/projects/b1025/sdi0596/Covid19_Single_Cell_RNASeq_project/filtered_matrices_count/epithelial/filtered_feature_bc_matrix/",gene.column=1,cell.column=1)
habermann = Read10X(data.dir="/projects/b1025/sdi0596/Covid19_Single_Cell_RNASeq_project/filtered_matrices_count/Habermann_Single_Cell_Data/",gene.column=1,cell.column=1)

data.list=c(epithelial,habermann)

normalize and identify variable features for each dataset independently
ifnb.list <- lapply(X = data.list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

Perform integration
immune.achors= FindIntegrationAnchors(object.list = ifnb.list,max.features = 200,verbose = TRUE)

Error in FindIntegrationAnchors(object.list = ifnb.list, max.features = 200, :
Max dimension too large: objects 1, 2 contain fewer than 30 cells.
Please specify a maximum dimensions that is less than the number of cells in any object (4).

Regarding the dimension of these objects, the epithelial one is 33538 x 49109 and the habermann is 33694 x 220213.
Any insights would be highly helpful.
Thank you.

There's a few problems with your code, first, when you do Read10X() it returns you a sparse matrix, and you need to put this into a Seurat object with meta data, before doing the integration.
So for example, i create some example data that is similar to your output from Read10X() :
library(Matrix)
epithelial = Matrix(matrix(rnbinom(100000,mu=200,size=1),ncol=50),sparse=TRUE)
rownames(epithelial) = paste0("gene",1:2000)
colnames(epithelial) = paste0("e",1:50)

habermann = epithelial + rnbinom(100000,mu=20,size=1)
colnames(habermann) = paste0("h",1:50)

Then convert to a seurat object:
data.list=c(epithelial = CreateSeuratObject(epithelial),
            habermann = CreateSeuratObject(habermann))

Then find variable features:
data.list <- lapply(data.list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

Should look like this:
$epithelial
An object of class Seurat
2000 features across 50 samples within 1 assay
Active assay: RNA (2000 features, 2000 variable features)

$habermann
An object of class Seurat
2000 features across 50 samples within 1 assay
Active assay: RNA (2000 features, 2000 variable features)

Now perform integration, below I have to reduce k.filter because I have very little cells in this example. In your case, you can simply use the default settings.
features <- SelectIntegrationFeatures(object.list = data.list)
anchors= FindIntegrationAnchors( data.list,max.features = 200,
k.filter=50,k.anchor = 3,verbose = TRUE)

It is possible that some of your samples have too few cells, so it is recommended to check the data and remove samples with low cell counts.

As it’s currently written, your answer is unclear. Please [edit] to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center.


text:
Nextflow: Can I parameterize similar processes inside a for loop?
I have set of ~5 processes which perform very similar tasks such that their differences can be parameterized. Rather than putting 5 items through some channel (which seems complicated), or defining 5 separate processes, it appears to be possible in Nextflow to declare the processes in a for loop, and output into an array of channels.
Below is a prototypical example and it seems to work as intended. (I understand this could be done differently with some channel cross product operator but I am just giving this as an example of a process in a for loop.)
Is there anything wrong with doing this? Is there a better way to do it?
Bonus question - the name of the process below is "foo". I would like to have the index reflected in the name i.e. "foo_$i" but this doesn't seem possible by default - any way to do this??
some_inputs = ['prot','dna', 'rna']

some_parameters = [5,3,2,7,1]

// Set up an array of channels to output to
receivers = new ArrayList<Channel>()
for (i = 0; i < 5; i++) {
    receivers.add(Channel.create())
}

// Iterate 5 times to create 5 processes each with different parameters
for (i = 0; i < 5; i++) {

    // WOW - process can be defined inside a for loop!!
    process "foo" {

      // Extract parameters here so correct $i value is used
      p = some_parameters[i]

      input:
      val x from Channel.from(some_inputs)

      output:
      val "$x $p" into receivers[i] // Note that correct index $i seems to be used here

      """
      """
    }
}

// Confirm that each of the processes outputted independently with the correct parameter
receivers[0].view { "Received[0]: $it" }
receivers[1].view { "Received[1]: $it" }
receivers[2].view { "Received[2]: $it" }
receivers[3].view { "Received[3]: $it" }
receivers[4].view { "Received[4]: $it" }

Output - note that the parameter is consistent within each input channel i.e. working correctly:
N E X T F L O W  ~  version 20.10.0
Launching `play.nf` [fabulous_keller] - revision: c89e2b2281
executor >  local (15)
[10/b477ec] process > foo (3) [100%] 3 of 3 ✔
[16/8282ba] process > foo (2) [100%] 3 of 3 ✔
[db/f7327e] process > foo (1) [100%] 3 of 3 ✔
[fb/886f6f] process > foo (1) [100%] 3 of 3 ✔
[46/c34ed3] process > foo (1) [100%] 3 of 3 ✔
Received[1]: rna 3
Received[1]: prot 3
Received[4]: rna 1
Received[2]: rna 2
Received[0]: dna 5
Received[2]: dna 2
Received[2]: prot 2
Received[0]: prot 5
Received[3]: rna 7
Received[0]: rna 5
Received[1]: dna 3
Received[3]: dna 7
Received[4]: dna 1
Received[3]: prot 7
Received[4]: prot 1

Use channels, this is what nextflow's entire design is based on. I'm honestly surprised that spawning 5 processes with the same name in a loop even works.
Turn your parameters into regular channels, use the combine operator (https://www.nextflow.io/docs/latest/operator.html#combine) to transform the 2 parameter channels into a channel that has each combination of params/inputs. Feed this channel into the  foo  process. Afterward use branch (https://www.nextflow.io/docs/latest/operator.html#branch) to split your foo output channel into 5 different channels.
typechan = Channel.from(5,3,2,7,1)
numchan = Channel.from('prot','dna', 'rna')

cartesian = numchan.combine(typechan)

process foo {
  input:
  val numtype from cartesian
  output:
  val numtype into bar
  script:
  """
  echo $numtype
  """
}

bar.branch {
  one: it[1] == 1
  two: it[1] == 2
  three: it[1] == 3
  five: it[1] == 5
  seven : it[1] == 7
}.set {result}
result.one.view()
```

You don't need to create processes in a loop when you have channels. The each qualifier allows you to repeat the execution of a process for each item in a collection. For example, the contents of 'test.nf':
some_inputs = Channel.of( 'prot', 'dna', 'rna' )
some_parameters = Channel.of( 5, 3, 2, 7, 1 )

process foo {

    tag { "task: ${task.index} - ${x},${y}" }

    input:
    val x from some_inputs
    each y from some_parameters

    output:
    tuple val(x), val(y) into receivers

    """
    echo "x: ${x}; y: ${y}"
    """
}

results = receivers.branch { x, y ->

    five: y == 5
        return tuple( x, y )

    three: y == 3
        return tuple( x, y )

    two: y == 2
        return tuple( x, y )

    seven: y == 7
        return tuple( x, y )

    one: y == 1
        return tuple( x, y )
}

results.five.view { "My result for param 5 is: $it" }
results.three.view { "My result for param 3 is: $it" }
results.two.view { "My result for param 2 is: $it" }
results.seven.view { "My result for param 7 is: $it" }
results.one.view { "My result for param 1 is: $it" }

Then running like:
nextflow run -ansi-log false test.nf

Results:
N E X T F L O W  ~  version 20.10.0
Launching `test.nf` [fabulous_franklin] - revision: 14c3290d2b
[6e/c12ce7] Submitted process > foo (task: 7 - dna,3)
[93/616f93] Submitted process > foo (task: 4 - prot,7)
[27/0bc3c4] Submitted process > foo (task: 3 - prot,2)
[0d/0d0b51] Submitted process > foo (task: 9 - dna,7)
[d4/f56d18] Submitted process > foo (task: 6 - dna,5)
[b0/d1535f] Submitted process > foo (task: 8 - dna,2)
[13/0a660a] Submitted process > foo (task: 5 - prot,1)
[f0/1c2f57] Submitted process > foo (task: 2 - prot,3)
[3f/3e893d] Submitted process > foo (task: 1 - prot,5)
[19/2921de] Submitted process > foo (task: 10 - dna,1)
My result for param 3 is: [dna, 3]
My result for param 7 is: [prot, 7]
My result for param 7 is: [dna, 7]
[a7/953a2b] Submitted process > foo (task: 11 - rna,3)
[1c/0fd4ef] Submitted process > foo (task: 12 - rna,5)
My result for param 1 is: [prot, 1]
[39/d21adf] Submitted process > foo (task: 14 - rna,7)
My result for param 5 is: [dna, 5]
My result for param 2 is: [prot, 2]
My result for param 2 is: [dna, 2]
[f6/8509a7] Submitted process > foo (task: 15 - rna,1)
My result for param 5 is: [prot, 5]
[6f/cee0ae] Submitted process > foo (task: 13 - rna,2)
My result for param 3 is: [rna, 3]
My result for param 1 is: [dna, 1]
My result for param 2 is: [rna, 2]
My result for param 3 is: [prot, 3]
My result for param 1 is: [rna, 1]
My result for param 5 is: [rna, 5]
My result for param 7 is: [rna, 7]

The answer to your bonus question is to use the task.index variable in your script. I used this to tag your job, but it can be used anywhere in the process declaration. You likely won't need this variable though.
Usually, using the branch operator is enough for most use cases. However,it's possible to implement your own branching function. The following uses the new DSL, DSL2, to automatically fork the 'receivers' channel:
nextflow.enable.dsl=2

process foo {

    tag { "task: ${task.index} - ${x},${y}" }

    input:
    val x
    each y

    output:
    tuple val(x), val(y)

    """
    echo "x: ${x}; y: ${y}"
    """
}

def branchByElement( input_ch, elements, by=0 ) {

    def targets = [:]
    for( element in elements ) {
        def key = element
        targets[ key ] = { input_ch.filter { it[by] == key } }
    }
    return targets
}

workflow {

    some_inputs = [ 'prot', 'dna', 'rna' ]
    some_parameters = [ 5, 3, 2, 7, 1 ]

    inputs_ch = Channel.fromList( some_inputs )
    parameters_ch = Channel.fromList( some_parameters )

    foo( inputs_ch, parameters_ch )

    parameter_results = branchByElement( foo.out, some_parameters, by=1 )
    input_results = branchByElement( foo.out, some_inputs, by=0 )

    parameter_results[7].call().view { "Param 7: $it" }
    input_results.dna().view { "Input dna: $it" }
}

Results:
N E X T F L O W  ~  version 20.10.0
Launching `test.nf` [distraught_meninsky] - revision: 04df251752
[bb/32325f] Submitted process > foo (task: 4 - prot,7)
[a0/19fdf5] Submitted process > foo (task: 8 - dna,2)
[e9/4c1d09] Submitted process > foo (task: 6 - dna,5)
[a3/dcb002] Submitted process > foo (task: 2 - prot,3)
[db/7e612e] Submitted process > foo (task: 5 - prot,1)
[cd/dfa84e] Submitted process > foo (task: 9 - dna,7)
[d4/6a85f2] Submitted process > foo (task: 7 - dna,3)
[f0/8e8bf4] Submitted process > foo (task: 1 - prot,5)
[89/f20a28] Submitted process > foo (task: 3 - prot,2)
Input dna: [dna, 2]
Param 7: [prot, 7]
Input dna: [dna, 5]
[53/056fe4] Submitted process > foo (task: 10 - dna,1)
[f3/a61abf] Submitted process > foo (task: 12 - rna,3)
Input dna: [dna, 7]
[f9/a568a8] Submitted process > foo (task: 13 - rna,2)
Param 7: [dna, 7]
[06/72bb6c] Submitted process > foo (task: 11 - rna,5)
Input dna: [dna, 3]
[a3/6d80d8] Submitted process > foo (task: 14 - rna,7)
[ac/acdb5a] Submitted process > foo (task: 15 - rna,1)
Param 7: [rna, 7]
Input dna: [dna, 1]

I think each is even more elegant than combine, do you also have a more elegant way to split the resulting channel than branch ?

@Pallie I think there is using a function and the new DSL2 - please see my update above.


text:
correlation coefficient versus DEGs analysis: what's the best approach for low expressed genes?
I have 6 folders. Each one contains 7 datasets of a specific type of cancer (RNA-Seq) and 7 datasets of normal tissue (healthy control for that type of cancer). A total of 84 datasets.
I want to investigate my gene of interest in each type of cancer.
I believe the best strategy is to merge the FPKM columns and perform correlation analyses of my gene of interest (GOI) versus all other genes in both cancer and normal tissues. As a result, I will have 14 correlation tables, one for each type of cancer and one for each control tissue. Then, I can investigate the pathways related to the genes most correlated with my GOI in each type of cancer and write my biological interpretation.
Am I in the right direction? Is this a good strategy to answer to my question?
I could also compare the control and cancer samples to investigate the significant DEGs (fold change > 1) and the pathways activated by these DEGs, but my gene of interest doesn't have a strong pattern of gene expression, so I'm afraid it won't be in the list of statistically significant DEGs, that's why I'm preferring the correlation analysis approach.
Any other idea?
The controls were collected from the same patient and processed at the same time. I have the gene counts and can use EdgeR, but what should I do if my GOI doesn't appear among the statistically significant DEGs (fold change > 1)? How can I investigate the role of one single gene in a RNA-Seq dataset if it is not differentially expressed? It can have a role in cancer, even if it is not an important DEG. That's why I thought correlation analyses could give me interesting information.

First of all, RNASeq is extremely sensitive to batch effects.  Are these matched controls processed at the same time as the tumors?  IF they are not, lots of your differences will be caused by batch effect, and not cancer.

I believe the best strategy is to merge the FPKM columns and perform
correlation analyses of my gene of interest (GOI) versus all other
genes in both cancer and normal tissues.

What you propose isn't publishable.  Use the software that people have been using and validating for years.  Use either DESeq2 or EdgeR, or limma if you don't have access to raw gene counts.
If the underlying counts are low, the correlation won't mean much, because your real values could be quite off. If your gene has 2 counts in one sample, and 3 counts in another, that does not mean it really correlates with a gene that has 500 reads in the one sample and 750 in the other. I don't know that you can trick your way around it; you need more reads, or more samples for more power.

thank you very much for you answer, @swbarnes2


text:
Why am I getting empty expression data from GEO?
I am trying to analyze the scRNAseq data from this study. In their Method section they write:

The accession number for the RNA and DNaseq data reported in this
paper is GEO:
GSE116237.

When I go ahead and try to pull this data in R using GEOquery, I am faced with a data frame with no features (i.e. no genes). The following code
library("GEOquery")
gse <- getGEO("GSE116237")
eset1 <- gse[[1]]
eset2 <- gse[[2]]

print(dim(eset1))
print(dim(eset2))

yields
Features  Samples
       0      184
Features  Samples
       0      681

I also didn't manage to download the data manually; I get an error message when trying to unzip the .zip files or get redirected to SRA Run Selector from which I also didn't manage to figure out how to download stuff.
Am I missing something or are these files simply corrupted? The paper is published in Cell, so I would assume this data is, to some extend, legit...

I am sure the data is legit, you are just approaching it incorrectly. getGEO is an application for microarray data, not for digital count data such as (sc)RNA-seq, therefore what you aim to do is simply not possible by design. Unless you want to start from the raw reads why not taking the at the file named GSE116237_scRNAseq_expressionMatrix.txt.gz provided at the bottom of https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE116237?


text:
Timeout when downloading the ncbi nr blast database
I am experiencing timeout problems when downloading the NCBI nr preformatted blast database using the update_blastdb script (version 504861).
I run the script with the following paramters
update_blastdb --decompress --passive --verbose nr

and I get the following error message (in verbose mode)
Downloading nr (45 volumes) ...
Downloading nr.00.tar.gz...Net::FTP=GLOB(0x5610fb59b8f8)>>> PASV
Net::FTP=GLOB(0x5610fb59b8f8)<<< 227 Entering Passive Mode (165,112,9,229,195,144).
Net::FTP=GLOB(0x5610fb59b8f8)>>> RETR nr.00.tar.gz
Net::FTP=GLOB(0x5610fb59b8f8)<<< 150 Opening BINARY mode data connection for nr.00.tar.gz (18745730730 bytes)
Net::FTP: Net::Cmd::getline(): timeout at /usr/share/perl/5.26/Net/FTP/dataconn.pm line 82.
Unable to close datastream at /usr/bin/update_blastdb line 202.
Net::FTP=GLOB(0x5610fb59b8f8)>>> PASV
Net::FTP: Net::Cmd::getline(): unexpected EOF on command channel: Connection reset by peer at /usr/bin/update_blastdb line 203.
Failed to download nr.00.tar.gz.md5!
Net::FTP: Net::Cmd::_is_closed(): unexpected EOF on command channel: Connection reset by peer at /usr/bin/update_blastdb line 101.
Net::FTP: Net::Cmd::_is_closed(): unexpected EOF on command channel: Connection reset by peer at /usr/bin/update_blastdb line 101.

The timeout happens after ~35 minutes and a file that is approximately 18GB big is being downloaded, which matches the expected filesize. The checksum file (nr.00.tar.gz.md5) is not downloaded. So I'm not sure which of the two files is actually the problem.
I tested downloading the nt database and everything seems to work fine, so I don't think the script is the problem. For comparison, this is the nt download output for the first file
Downloading nt.00.tar.gz...Net::FTP=GLOB(0x562575a73168)>>> PASV
Net::FTP=GLOB(0x562575a73168)<<< 227 Entering Passive Mode (165,112,9,229,196,51).
Net::FTP=GLOB(0x562575a73168)>>> RETR nt.00.tar.gz
Net::FTP=GLOB(0x562575a73168)<<< 150 Opening BINARY mode data connection for nt.00.tar.gz (4065912989 bytes)
Net::FTP=GLOB(0x562575a73168)<<< 226 Transfer complete
Net::FTP=GLOB(0x562575a73168)>>> PASV
Net::FTP=GLOB(0x562575a73168)<<< 227 Entering Passive Mode (165,112,9,229,195,107).
Net::FTP=GLOB(0x562575a73168)>>> RETR nt.00.tar.gz.md5
Net::FTP=GLOB(0x562575a73168)<<< 150 Opening BINARY mode data connection for nt.00.tar.gz.md5 (47 bytes)
Net::FTP=GLOB(0x562575a73168)<<< 226 Transfer complete

Any help would be appreciated.

The answer of the support stated that it is likely that the size of the first download is the problem and that I should try using another tool like rsync. Because I'm unfamiliar with rsync I decided to write a python script that did the job.
import urllib, tarfile, json

base_url = "https://ftp.ncbi.nlm.nih.gov/blast/db"

# download manifest file to get the filenames
manifest_file = 'blastdb-manifest.json'
urllib.request.urlretrieve(f"{base_url}/{manifest_file}", manifest_file)

with open(manifest_file) as f:
    manifest_data = json.load(f)

# download everything
for file in manifest_data['nr']['files']:
    # download checksum file
    checksum_file = f"{file}.md5"
    urllib.request.urlretrieve(f"{base_url}/{checksum_file}", checksum_file)

    # download the archive file
    urllib.request.urlretrieve(f"{base_url}/{file}", file)

    # check that the checksums match
    calculated_checksum = get_md5_for_file(file)

    with open(checksum_file) as f:
        for line in f:
            line = line.strip()
            if line.endswith(file):
                checksum = re.split("\s+", line)[0]
                if checksum == calculated_checksum:
                    checksum_matches = True
                else:
                    raise Exception(f"Checksum doesn't match expected for file {file}")

    # unpack the archive
    tar = tarfile.open(archive_file)
    tar.extractall(file)
    tar.close()

The function to calculate the md5 checksum is the following (borrowed from this question)
def get_md5_for_file(file):
    hash_md5 = hashlib.md5()
    with open(file, 'rb') as f:
        for chunk in iter(lambda: f.read(4096), b""):
            hash_md5.update(chunk)

    return hash_md5.hexdigest()

Note that this is not the exact script I use, so I can't guarantee it works exactly like this, but I tried to include all the important stuff.
In my version I include checks to whether my existing blast database needs to be updated (I keep the manifest files and compare the manifest['nr']['last_updated'] values) and I also try to skip already unpacked archives while iterating (check if the nr.xx.phd file exists) to save time in case the script fails due to problems in the download.

Same here. I used axel to download it a week ago, but fail to get finished. Also, the downloaded taxa zipped file (60Mb) doesn't match md5sum. Apparently, they update the FTP files every 30 minutes or 1 hour.

Ok, good to know. I wrote them an email and will update when I get an answer

They're usually pretty responsive. Have you heard back yet?

Sorry for the late reply, I posted an answer.


text:
how to sort and join two files based on First Column Id
I have got two files large files like this in Tab-delimited format, trying to merge in R, do we need to sort the file before merging?
 otu <- read.table("asv.txt", sep="\t", header=T, blank.lines.skip=F, check.names=F)
> head(otu)
                ID                   Soil.9.S42 Soil-27.S33.L001 Soil-45.S54.L001
1 d2ec9f3b77975c0f457e4b7413b217ff         34               51               84
2 3147790f0d5a78316fb9dd64f53b9473         28               32               52
3 97aecc1f35cc1f50db507ad71dd22367         33                0                7
4 bfad6370d28182cc6304844e9bec7fb6        545              237              127
5 5fa2a987221a1d9ca416148570c18086          0                0                0
6 8e60d301122d7aa359eb6b0b00f37f62        124              143              128

> taxa <- read.table("taxonomy.txt", sep="\t", header=F ,stringsAsFactors=F,quote="")
> head(taxa)
                  ID                     V2             V3                  V4
d2ec9f3b77975c0f457e4b7413b217ff Bacteria  Cyanobacteria    Oxyphotobacteria
3147790f0d5a78316fb9dd64f53b9473 Bacteria Proteobacteria Alphaproteobacteria
97aecc1f35cc1f50db507ad71dd22367 Bacteria Proteobacteria Gammaproteobacteria
bfad6370d28182cc6304844e9bec7fb6 Bacteria Proteobacteria Alphaproteobacteria
5fa2a987221a1d9ca416148570c18086 Bacteria Proteobacteria Alphaproteobacteria
8e60d301122d7aa359eb6b0b00f37f62 Bacteria Proteobacteria Deltaproteobacteria

Now having issue with sorting file but merge command work perfectly;
newotu <- otu[order(ID),]
Error in order(ID) : object 'ID' not found
> total <- merge(otu,taxa,by="ID")

Thanks for your help.

What do you mean 'it didn't work'?

Also I'm confused - you seem to be reading the files into R and then using unix commands to merge? I prefer to do this sort of things in R - this should help https://www.statmethods.net/management/merging.html

I just use R showing the file content but wanted to do using unix command or R command.

OK - it's better to try in R IMO. Did you try follow the link I posted?

Error while doing in R, updated the command and error in the post.

Do we need to sort file in R before merging?

There was some issue with the command like header True or False, now fix it and trying again, but do we need to sort the file in R before merging?

no issue now while merging but when trying to sort file it shows error;

it doesn't need to be sorted

Thank you so much for your help. this is a very easy way to merge files in R. I also checked the file and it merges perfectly.

@user438383 or bioinfonext could one of you please post an answer with the solution you found so this question can be marked as answered?

You should check left_join() method in dplyr package. In this way you do not need to sort the data frames, as the function automatically take care of this.
For example: merged_df <- left_join(otu, taxa, by = "ID").


text:
How to access XL-MS raw data and databases?
I'm looking for databases and/or raw data of protein mass spectrometry experiments.
After reviewing xlinkdb and Proxl, I was curious to learn about other alternatives and options that allow users to run the database locally and browse raw experimental data. NCBI and proteomexchange.org don't seem to have enough data either.


text:
Use conditional in workflow in Nextflow dsl2
How should I do:
...
if ($(params.aligner) == "bowtie") {
    align_bowtie2(get_fq_info.out.fq_info)
}
else {
    align_bwa(get_fq_info.out.fq_info)
}
...

You can move that into your workflow definition. There's an example in the docs of how best to use workflow parameters:

A workflow component can access any variable and parameter defined in
the outer scope:
params.data = '/some/data/file'

workflow my_pipeline {
    if( params.data )
        bar(params.data)
    else
        bar(foo()) }

Your example workflow would look like:
workflow {

    fq_info = get_fq_info.out.fq_info

    if( params.aligner == "bowtie" ) {
        align_bowtie2( fq_info )
    }
    else {
        align_bwa( fq_info )
    }

    ...
}
```

Is this inside the main.nf file or inside the process.nf file? Anyway, the param.anything variable should already be available in the nextflow context when executed, thus, you could simply do:
if (params.aligner == "bowtie") {
    align_bowtie2(get_fq_info.out.fq_info)
} else {
    align_bwa(get_fq_info.out.fq_info)
}

Check out this main.nf of a nextflow pipeline of mine, it contains some working if/else snippets.

Here it is another StackOverflow question on conditional nextflow pipelines that might help.

Hope it helps.


text:
splitCsv map with text, join with another channel, nextflow
I have a csv file which describes number of chunks per chromosome, like;
chrom,chunks
chr1,"CHUNKS=([1]=""0e6 3e6"" [2]=....)"
chr2,"CHUNKS=([1]=""0e6 3e6"" [2]=....)"
chr3,"CHUNKS=([1]=""0e6 3e6"" [2]=....)"
.....

I could crate a channel as follows,
channel
    .fromPath('chunks3.csv')
    .splitCsv(header:true)
    .map{row -> tuple(row.chrom,row.chunks)}
    .set{chunk_ch}

It gives me chromosome with respective chunk definition.
But when I try to join this channel with another channel, it does not work.
For example I have a channel "phased". It is a output channel created from a process.
channel
    .fromFilePairs('/some/path/genetic_map_chr*_combined_b37.txt',size:1)
    .map {group_key, file_list -> tuple(group_key.replaceFirst(/^*genetic_map_/,""),file_list.first())}
    .set {genetic_map}

channel
    .fromFilePairs('/some/path/*hrc.r1-1.ega.grch37.chr*.{hap.gz,legend.gz,samples}',size:3)
    .map {group_key,file_list -> tuple(group_key.replaceAll(/^.*.hrc\.r1-1\.ega\.grch37\.|\.haplotypes\.noibd|\.haplotypes/,""),file_list)}
    .set {ref_panel}

qc1_chr
    .map{chrom,bfiles -> tuple(chrom,bfiles.findAll{it.baseName.endsWith('step10')})}
    .join(ref_panel).join(genetic_map)
    .into {shape_in;imp_in}

process phasing {
    publishDir params.out, mode:'copy'

    input:
    tuple val(chrom), path(input_bed), path(input_ref), path(input_map) from shape_in

    output:
    tuple val("${chrom}"),path("${chrom}.*") into phased

    script:
    def (bed, bim, fam) = input_bed
    def (haplotype, legend, sample) = input_ref
    """
    shapeit --thread "${task.cpus}" --input-bed "${bed}" "${bim}" "${fam}" --input-map "${input_map}" --input-ref "${haplotype}" "${legend}" "${sample}" --duohmm --
    output-max "${chrom}.step10" --seed 54321
    """
}

Which emits something like:
phased.view()
[chr1,[file1]
chr2,[file2]
chr3,[file3]....]

If I try chunk_ch.join(phased) it gives nothing.
Any help?

Below is a reproducible example. I tried to simulate the same channels using chr1 as an example.
#!/usr/bin/env nextflow

// Simulate the chunk channel
channel
    .fromPath('chunk.csv')
    .splitCsv(header:true)
    .map{row -> tuple(row.chrom,row.chunks)}
    .set{chunk_ch}

// Simulate chr and chr_path as a single entry list
channel.from(["chr1", "chr1.txt"]).toList().map{tuple(it[0], it[1])}.set{shape_in}

// uncomment to view the content
// shape_in.view()

// process of phased_ch. I modified it to output the chr and chr_path as given in your example
process phased_ch {
  echo true
  input:
  tuple val(chr), val(chr_path) from shape_in
  output:
  tuple val("${chr}"), path("${chr_path}") into phased
  script:
  """
  echo ${chr}
  touch ${chr_path}
  """
}

// joining the two channels to produce the correct join
chunk_ch.join(phased).view()

// will produce
// [chr1, "CHUNKS=([1]=""0e6 3e6"" [2]=....)", /media/Documents/Nextflow/work/8e/2fcd9a00b04f0f5d44711068d68a0f/chr1.txt]

chunk.csv contains
chrom,chunks
chr1,"CHUNKS=([1]=""0e6 3e6"" [2]=....)"
chr2,"CHUNKS=([1]=""0e6 3e6"" [2]=....)"
chr3,"CHUNKS=([1]=""0e6 3e6"" [2]=....)"

In your case, you should verify what are the contents of chunk_ch and phased on output. The first key of join() should be similar otherwise it won't give anything. From the documentation of join:

The join operator creates a channel that joins together the items
emitted by two channels for which exists a matching key. The key is
defined, by default, as the first element in each item emitted.


text:
Unexpected error using "stdin" in bedtools intersect as a piped command
From the bedtools intersect man (for version 2.30.0):
-a BAM/BED/GFF/VCF file “A”. Each feature in A is compared to B in search of overlaps.
   Use “stdin” if passing A with a UNIX pipe.

The stdin part, however, does not apply to version 2.27.1. The help on this version gives no advice for piping.
bam mem ref fq1.fastq fq2.fastq | samtools view -b | bedtools intersect -a stdin -b blacklist.bed -v > align_black.sam
Error: unable to open file or unable to determine types for file stdin

- Please ensure that your file is TAB delimited (e.g., cat -t FILE).
- Also ensure that your file has integer chromosome coordinates in the
  expected columns (e.g., cols 2 and 3 for BED).

I also tried:
(map to bam) | bedtools intersect -a "stdin" -b blacklist.bed -v > align_black.sam
(map to bam) | bedtools intersect -a -b blacklist.bed -v > align_black.sam
(map to bam) | bedtools intersect -b blacklist.bed -v > align_black.sam

But these don't work either, either for the same or a different reason. Here is a site that tells me to use the above syntax: https://bedtools.readthedocs.io/en/latest/content/example-usage.html
What is the solution to this if I want to pipe everything together for my users?

With v2.27.1 (or any version I ever used) instead of stdin use -
(...) | bedtools intersect -a - -b blacklist.bed -v > out.bam

Be aware that the output of this is a BAM, not SAM file. It would be easier though to simply make a list of regions you want to keep, e.g. the complement of the backlist file with the entire genome (bedtools complement), and then use samtools view -L option to only keep overlapping reads. That saves you the run through bedtools.

Thanks, I wonder why this isn't in the docs. I had my own workaround before you commented, but I'll accept this as the answer. Is there a reason you suggest not using bedtools?

@Jeff As samtools can do the same and you anyway run it you save one additional tool/pipe/core. That is the only reason. Not sure how the performance compares.

Here's a way to skip conversion to BAM and pipe between tools, which performs the same set operations:
$ bam mem ref fq1.fastq fq2.fastq | sam2bed - | bedops -e 1 - blacklist.bed > blacklisted_reads.bed

With this toolkit, the hyphen (-) follows a Unix convention of serving as a placeholder for stdin.

While I have chosen ATpoint as the official answer, another possible syntax to use is
bedtools intersect -a stdin -b blacklist.bed -v < `(map to bam)` > align_black.bam

The main issue with this syntax is that it obfuscates the direction of flow for the pipeline, piping step 2 after step 1 is better in the official answer.


text:
Find common SNPs in multiple VCF files
I have 3 VCF files.
A1.vcf
A2.vcf
A3.vcf

I want to get the common SNPs that are present in all these three files. And output must be in vcf format.
output:
common_A.vcf

Do you mean SNPs in terms of shared variants or shared positions?

google intersect or overlap between VCF files, -> bcftools isec

I have tried isec, but the produced files do not contains the VCF file contains common file

A potential solution is to compress and index the files using bgzip and tabix, then use bcftools isec to find SNPs common to all three files, e.g.
for files in A1.vcf A2.vcf A3.vcf
  do
  bgzip "$files"
  tabix "$files"
done

bcftools isec -n=3 -c all -o all_common_SNPs.vcf A1.vcf.gz A2.vcf.gz A3.vcf.gz
```

Hi @Jared_mamrot - can you please explain -c all I've seen -c none -c all/both I'm unable to understand/follow this flag.

As I understand it @Death-Metal, if you have two vcfs, one with the variant T->CCA at chr1 pos10001 and the other with the variant T->A at chr1 pos10001 , the -c all flag will treat these as a 'common to the two files' (i.e. there is a variant at chr1:10001 and it will be 'kept' when you isec). If you specify -c none, the REF/ALTs need to match exactly (i.e. both vcfs need to be e.g. T->A to be kept). With -c both, the variant must be the same type e.g. both SNPs or both INDELs, regardless of the ALTs. More details here: https://samtools.github.io/bcftools/bcftools.html#common_options


text:
Sequence allignment with suffix array?
I'm writing the code from scratch and not using libraries for the actual indexing/search since it's a project for school. Any blackbox advice is fairly useless.
I'm planning on aligning some short reads with the D. melanogaster genome but, since I was thrown into this project with very little info/ experience I'm stuck on a few things.
Firstly, for preprocessing the fasta file:

What should I do with Ns (ambiguous bases), I've seen some people removing them or replacing them by gaps but this does not seem like the same information? Should they be left in?

As far as I understand, the lowercase letters stand for potential repeats. However, If I just ignore this, it is possible my sequences might not align to the repeat because it is not actually stated as a repeat. E.g: "ATATATAT" will not align to 'at'.

As far as I understand the difference between using suffix array and suffix tree: suffix tree will search in linear time (but is memory intensive), while suffix array will be space efficient (but searches through ALL sequences each time). I'm leaning towards suffix array because I'm working on my laptop and memory is kind of limited. I'm also thinking that, since you sort the array lexicographically, you really only have to check the sequences if they start with the same first base as the read you're searching for?
Another question on the arrays: if my genome has 100M base pairs, the array will have 100M entries, some of them being close to 100M bases long? Should I be using the seed and extend method?
On the other hand, since the fruit fly genome is one of the smaller genomes, maybe it would be feasible to use a tree?
Is there any way to validate the alignment quickly at the end? How much error should I be allowing?
This were a lot of question, I've been trying to read up on things but I can mostly only find applications using biopython etc.
Any help would be hugely appreciated!

STAR uses suffix arrays to map RNA-seq reads. It is both very fast and quite accurate.  Perhaps their manuscript may help with your questions.

AH, sorry I did not notice the "it's a homework" part. However, note that practically the mapping tools out there are open source, so you can actually browse the code and check what exactly they did.

What should I do with Ns (ambiguous bases), I've seen some people removing them or replacing them by gaps but this does not seem like the same information? Should they be left in?
As far as I understand, the lowercase letters stand for potential repeats. However, If I just ignore this, it is possible my sequences might not align to the repeat because it is not actually stated as a repeat. E.g: "ATATATAT" will not align to 'at'.

Good questions. I think this depends on the scope of your class project. Frankly, I think it would be best to just start with unambiguous sequences (only ATGC) and not worry about potential repeats. It's also a good software engineering exercise to solve the easy case first in a way that's flexible and allows you to expand to the additional complexities above. For example, the original Bowtie software considered references with ambiguous characters invalid. You can see for yourself here how Bowtie2 handles it.

As far as I understand the difference between using suffix array and suffix tree: suffix tree will search in linear time (but is memory intensive), while suffix array will be space efficient (but searches through ALL sequences each time). I'm leaning towards suffix array because I'm working on my laptop and memory is kind of limited. I'm also thinking that, since you sort the array lexicographically, you really only have to check the sequences if they start with the same first base as the read you're searching for?

Yes, the most naive implementation of a search with a suffix array will be $O(m\log n)$ time, where $m$ is the read length and $n$ is the genome length. However, with enhanced suffix arrays you can also obtain the same $O(n)$ search complexity of suffix trees. There is also the use of suffix arrays to construct Burrows-Wheeler transforms and even further compression with the FM index. There are open source implementations of all of these data structures that may be interesting to you.

Another question on the arrays: if my genome has 100M base pairs, the array will have 100M entries, some of them being close to 100M bases long? Should I be using the seed and extend method?

Suffix arrays store the indices of the start of each suffix in order to save space, so each entry is a single integer.

On the other hand, since the fruit fly genome is one of the smaller genomes, maybe it would be feasible to use a tree?

Its worth trying to compute this yourself (: You can determine the number of bytes required per base pair in your input and see how that scales with your available memory. For example, if your suffix tree  implementation required $20n$ bytes, you'd end up w/ roughly 2Gb needed for a fruit fly genome.

Is there any way to validate the alignment quickly at the end? How much error should I be allowing?

Again, this depends on the scope of the project. Allowing for errors in mapping reads is an active area of research. See the supplemental of the Bowtie2 paper for a dynamic programming algorithm that handles a certain number of mismatches.

This were a lot of question, I've been trying to read up on things but I can mostly only find applications using biopython etc.

All good questions, and I think many answers depend on the scope of the project. Reimplementing a read mapper is no small feat, especially when you start considering all of the different complexities such as sequencing error and ambiguous bases. I would strongly consider starting with looking at just the pure string problem of read mapping: Only consider ATCG in both reference and reads, and don't even consider reads who do not have exact matches to the reference. From there, time permitting, you can start adding on to your solution. As far as reading goes, there's a deep rabbit whole of these topics, from tools like bwa-mem, bowtie2, mini-map, to the more theory-oriented papers discussing optimized versions of the data structures. I can attach some here if you would like. Hope these answers helped!

Thank you! This was really helpful

@nibs feel free to accept the answer if it sufficiently answers your questions (:

Suffix arrays vs Suffix tree is a problem quite often discussed in bioinformatics classes, but it's not that much used in practice. Most of people use already optimized mappers bowtie2 or bwa-mem.
I would advise you against your intention of mapping the reads using your own implementation of anything. You don't need to reinvent the wheel, just use an existing solution.
By the way, the Ns in the reference are either hard-masked nucleotides (that would be repeats of some sorts) or (more often) simply unknown nucleotides. So there is no way to map your reads there anyway (and the mentioned mappers will handle the Ns for you).

Good advice, however op stated that this is homework, so the wheel must be reinvented (:

Yeah, i did not notice that tag when I wrote the answer. Well, I guess it does not hurt it's here :-).


text:
Normalize RNA seq data from multiple runs for expression analysis
I have RNA samples sequenced with TruSeq Stranded Total RNA kit protocol in Illumina HiSeq (2x125bp) and NovaSeq platforms (2x150bp) - almost 100 samples altogether. I have to use the samples data for expression analysis. Question is how do we normalize samples from different runs and of read length? What is best way (with available packages) we can normalize the count (to overcome the bias) and do a DE analysis? ANy suggestion would be helpful.

There are a number of methods. If you're doing DE, you have ComBatSeq, SVAseq, RUVseq, BUSseq. You could also try Z-score normalization ($(x - \overline{x}) / \sigma$) or even quantile transformation. For the latter two, make sure you work on each batch individually, not the whole dataset at once. See more on that here.
To visually compare if transformation works, plotting PCA/UMAP/t-SNE on raw and transformed data can perhaps be of some insight.

Thanks Roman for the response. Very useful. Do you recommend to trim the reads to 125bp for the NovaSeq to compare with the other batches from HiSeq? Or do you suggest to include read length as a factor during DE analysis? What do you suggest.

This sounds like another question. :) Feel free to post it as such.

Okay, let me do that. Don't know if this kind of experiment has been done previously (happens data is used from multiple sources for integrated analysis)

I agree with swbarnes that you don't need to explicitly model for different read length, but you need to normalize for different production of sequencing runs.
Among the several available packages, DEseq2 is one of the most widely used.
The normalization method computes the geometric mean of expression of each gene across samples; then, expression of the gene in each sample is divided by the geometric mean. Once this has been don for all genes, in each sample it computes the median of values, and this is the "size factor". Ok, it looks convoluted, but it is indeed easy, and it is implemented in DESeq2 and in other software packages.
You can find methodological details here.
Finally, DESeq2 also offer the possibility to account for additional experimental details (e.g. read length) and this should reduce the risk of spurious results.

Thanks for your response, Fabio. If I use ether median-to-ratios (geometric mean in DEseq2) or upper-quartile normalization, I think it should work good despite different sequencing machines? Also do you have any articles where the read length is used as a factor in the design for DE analysis? How robust is this method compared to trimming of reads as trimming on large samples is a huge additional task.

If you have different read length then trim all samples to the shortest length across batches, that will avoid mappability bias. Including it into the DE design is possible but as you can completely eliminate this bias in silico by trimming this is what I'd do.

Sequencing on different machines is usually not much of a confounding factor, but the fact that two machines were involved suggest different batches of experiments. I hope that batch is not confounded by group so each batch should contain samples of all groups to avoid a general batch effect that cannot be corrected for. Library prep should (in a well-designed experiment) be fully identical.

Use any of the established methods, be it RLE or TMM (they usually perform ery similar), and then perform PCA based on the top variable log2-transformed counts to explore whether there are batch effects that are worth correcting, so whether you see clustering or differences based on batch rather than experimental design in the early PCs. If so (given it is not confounded by experimental groups) you can include it into the DE design as a covariate.

As said by others, normalization and batch correction are independent processes. Also as I said in on top, you can eliminate read length bias (even though I doubt it is a major factor here) in silico which is what I'd recommend so you do not have to bother with it downstream.

I don't think RNASEq is that sensitive to running different things on different instruments.  Differing read lengths will change things.  You could always trim the reads down so they all match.  Or include read length as a factor in the design when looking for DE genes.

Normalization and batch-effect correction are two different things.
You always have to do normalization of RNAseq data in order to be able to compare the runs and do DE analysis, for instance. Plenty of methods and you have to search for them as normalization methods indeed.
What you are asking and concerned about is "batch-effects" in your data as a function of having used two different Illumina machines. This is a valid concern, but before you try to correct for it or mitigate it, you first have to assess whether it is actually there. Maybe is not at all. To figure out if you have batch-effects you have to first do data exploration, this is why the use of PCA/UMAP/t-SNE has been already suggested. In this plots you will be trying to see whether your data cluster together by machine type instead of by the biology it represents. Read about how to assess the presence of batch effects.
If you do observe it or anyway want to correct for its presence, then indeed you can try one of the methods also already mentioned: ComBatSeq, SVAseq, RUVseq, BUSseq.
Rather than the two different machines, the batch-effect might actually be due to different people doing the lib preparation, experimental settings if experiments were done at different moments, etc.


text:
Clustering with Seurat
My data consists of two groups control and stim. I'm concerned that my clustering in the stim group is dominated by the unregulated genes. I can remove these genes and recluster but is there a way to then project the removed genes onto the new umap space so see if they are indeed all clustered together. They are no longer in the seurat object so I can't think of how to do this. I hope this makes sense.

Yes, this can be done, particularly if your concern is with the clustering / UMAP rather than the transformation process; I've done this recently with some Rhapsody single cell data that we generated. The main trick is adding 'return.only.var.genes = FALSE', as well as calculating a variable set that is a bit larger than the desired set. Depending on the number of cells, this will balloon the size of the resulting Seurat object pretty significantly.
Note that UMAP is a technique to place cells, rather than genes. If you have the information on where the cells should go, then plotting gene expression on top of that via FeaturePlot should work without any reprojection needed.
Here's what I did (with experiment-specific variables excluded):
sct.data <- SCTransform(sct.data, vars.to.regress="percent.mt",
                        variable.features.n = 3000 + length(genesToExclude),
                        return.only.var.genes = FALSE);

sct.data <- RunPCA(sct.data, verbose = TRUE,
                   features=setdiff(VariableFeatures(sct.data), genesToExclude));
sct.data <- RunUMAP(sct.data);
sct.data <- FindNeighbors(sct.data);
sct.data <- FindClusters(sct.data);

All transformed genes should be available for querying in the results.


text:
Can an alternate allele be more common than a reference allele?
I'm trying to analyze data from a genotype-phenotype association study. The genotype data is stored in csv files that record for each variant site and each subject, the alleles that the subject has at that site. A quick online search reveals that 0 means the reference allele, 1 means the most common alternate allele, 2 means the 2nd most common alternate allele, etc.
But I noticed that at some variant sites, most (>90%) subjects have 1/1, i.e. two copies of the alternate allele. I asked around and some suggested I just swap the roles of reference and alternate alleles whenever I find that the alternate allele frequency > 0.5. My questions are:

Can there be a good reason to classify a less common allele as the reference allele? I thought that the reference allele was defined as the most common one in the population.
Is swapping the roles of reference and alternate alleles whenever I find that the alternate allele frequency > 0.5 a standard practice?

Thanks

You are mixing up a bit the terminology here. One thing is a reference allele and another is a major allele. The first refers to the genotype present at a given position in the "reference" genome (e.g. if you work with human data this can be the reference version hg19 or hg38, etc). The second refers to the most common allele in a population. This definition can lead to problems since what is the most common allele can differ among different populations. In this respect, it is perfectly fine to have a reference allele with a lower frequency than another alternate allele.

No for the reasons above. If you mean whether swapping major to minor is standard practice whenever you find out the major allele is actually minor as per the definition (less frequent), then I guess you can but this will lead to inconsistencies whenever you try to merge your data with other data sets, or other people merge their data with your data. Since this for sure happens and issues about strandness are well known, one has to make sure that you make explicit how you treat your data once you release it, and also one has to make sure that major/minor are defined the same way when merging data, and do the fix whenever needed.

Yes, depending on how the reference allele was determined. If, for example, you want to refer to hg38 as the reference genome, then it's important to use the hg38 allele as reference.

Without stating the reference allele convention? Not that I'm aware of. However many people refer to minor/major alleles rather than non-reference alleles, in which case it is assumed that the major allele is the most common allele.

As long as you state what your reference genome and/or population is, then it's fine to choose your own convention about reference alleles. But please be consistent and explicit about what you've done. One of my biggest struggles in GWAS has been working out the correct strand orientation for stated variants (important when merging datasets), and it doesn't help for complementary variants (i.e. A/T, C/G) when the ref/alt alleles are arbitrarily chosen.


text:
Number of sing-copy genes identified by orthofinder is far less than BUSCO assessment?
Dear Orthofinder users and BUSCO users,
Recently, I have assessed 30 genome annotataions using BUSCO with lineage dataset ascomycota_odb10, I got 1100+ complete and single-copy genes. But the orthofinder  only identified 300+ single-copy genes. Aneone else have encountered this similar problem?
Thanks,
Yanpeng

I know what caused the difference. BUSCO uses more strict criterion to determine single-copy genes but manly of the result genes are identified as multicopy genes by OrthoFinder using less strict parameters.

It's hard to say for sure without seeing the exact commands you used. However, Orthofinder identifies/defines single-copy orthologues as those genes present in a single copy in all species being compared in the analysis. If a BUSCO orthologue is present in more than one copy in one of the species being compared (could be multi-copy, could be fragmented etc.), then it won't be recovered as single copy by Orthofinder.

Thanks for your reply. I run OrthoFinder with default parameters. The BUSCO identified 1000+ complete single-copy genes but only 387 identified by OrthoFinder. I am trying to confirm whether BUSCO single-copy genes are true single-copy.