1
$\begingroup$

I have written a Nextflow script with three process:

  • The first process takes a pair of fastq files and aligns with reference genome. The process writes the resulting SAM file into sam channel.
  • Second process takes input from the sam channel and creates a BAM file from it, and writes it into bam channel.
  • Final process reads from bam channel and sorts the BAM files.

Following is the entire code:

#!~/bin nextflow

reads_ch = Channel.fromFilePairs('raw/*_{1,2}.fastq.gz', flat: true)
ref_genome = file('reference/genome/index/ref_genome')

process mapToRef {
    memory '20 GB'
    cpus 16
    input:
    tuple val(sample_id), file(read1), file(read2) from reads_ch
    
    output:
    tuple val(sample_id), file("${sample_id}.sam") into sam_ch

    script:
    """
    ~/biotools/bwa/bwa mem -t 60 ${ref_genome} ${read1} ${read2} > ${sample_id}.sam
    """   
}

process samToBam {
    memory '20 GB'
    cpus 16
    input:
    tuple val(sample_id), file("${sample_id}.sam") from sam_ch

    output:
    tuple val(sample_id), file("${sample_id}.bam") into bam_ch

    script:
    """
    ~/biotools/SAMTOOLS/samtools-1.12/samtools view -S -b ${sample_id}.sam > ${sample_id}.bam 
    """
}

process sortBam {
    memory '20 GB'
    cpus 16
    input:
    tuple val(sample_id), file("${sample_id}.bam") from bam_ch

    output:
    tuple val(sample_id), file("${sample_id}.sorted.bam") into sorted_bam_ch

    script:
    """
    ~/biotools/SAMTOOLS/samtools-1.12/samtools sort ${sample_id}.bam -o ${sample_id}.sorted.bam
    """
}

All three processes successfully execute their commands. However, they are executed sequentially. For example, the first process mapToRef runs sequentially for all read pairs, followed by samToBam, which too, runs sequentially after all mapToRefs have finished execution. I was under impression that parallelization would "work" out of the box, but I must be doing something wrong. I am a complete beginner in Nextflow and would greatly appreciate any help regarding this.

EDIT

Output of reads_ch.view():

[ERR2512393, /home/work/raw/ERR2512393_1.fastq.gz, /home/work/raw/ERR2512393_2.fastq.gz]
[ERR2512394, /home/work/raw/ERR2512394_1.fastq.gz, /home/work/raw/ERR2512394_2.fastq.gz]
[ERR2512391, /home/work/raw/ERR2512391_1.fastq.gz, /home/work/raw/ERR2512391_2.fastq.gz]
[ERR2512392, /home/work/raw/ERR2512392_1.fastq.gz, /home/work/raw/ERR2512392_2.fastq.gz]
[ERR2512390, /home/work/raw/ERR2512390_1.fastq.gz, /home/work/raw/ERR2512390_2.fastq.gz]
$\endgroup$
8
  • $\begingroup$ Before launching any process, it would be wise to check what does your input channel reads_ch contains. On the next row after defining ref_genome variable, write reads_ch.view() and check the output on the console to see what inputs you have. $\endgroup$
    – user324810
    Sep 4 at 19:16
  • $\begingroup$ @Law I have updated my question with the output of reads_ch.view() $\endgroup$
    – PPrasai
    Sep 4 at 19:26
  • $\begingroup$ Thank you @PPrasai, I can see that your input channels contain for each sample, 2 fastqs where as in your input for mapToRef, you are indicating file(sample_reads). In your Script, the command line should be ~/biotools/bwa/bwa mem -t 60 ${ref_genome} ${sample_reads[0]} ${sample_reads[1]} > ${sample_id}.sam $\endgroup$
    – user324810
    Sep 4 at 19:32
  • 1
    $\begingroup$ As a side note, why don't you make just one process with all three steps in the strips. Pipe bwa into samtools sort, and there is even now an option to write an index on the fly. Three processes are quite an overkill for such a simple task. Basically bwa mem (...) | samtools sort --write-index -o out.bam in a single process. No need for all these intermediate files and tons of log and cache files that are produced by three processes. $\endgroup$
    – ATpoint
    Sep 5 at 16:43
  • 1
    $\begingroup$ @PPrasai You have to declare files with the .bai suffix in output: because that is an additional file being generated so nextflow must be told about it first, like any other output that is produced. In DSL2 I always output this as a tuple so index and bam are kept together: github.com/ATpoint/atac_chip_preprocess/blob/main/modules/… Don't take my workflows as inspiration, I am still rather noobish with nf. $\endgroup$
    – ATpoint
    Sep 6 at 14:17
2
$\begingroup$

I have tested it on my end with the following script:

#!/usr/bin/env nextflow

reads_ch = Channel.fromFilePairs('/root/test_nextflow/*_{1,2}.fastq.gz', flat: true)


process mapToRef {
  memory '1 GB'
  echo true
  cpus 1

  input:
  tuple val(sample_id), file(r1), file(r2) from reads_ch

  output:
  tuple val(sample_id), file("${sample_id}.sam") into sam_ch

  script:
  """
  echo bwa mem ${r1} {r2}
  if [ ${sample_id} == "ERR2512393" ]; then sleep 10; fi
  touch ${sample_id}.sam
  """
}

process samToBam {
  memory '1 GB'
  cpus 1
  echo true

  input:
  tuple val(sample_id), file(sam) from sam_ch

  output:
  tuple val(sample_id), file("${sample_id}.bam") into bam_ch

  script:
  """
  echo samtools view -S -b ${sample_id}.sam
  touch ${sample_id}.bam
  """
}

process sortBam {
  memory '1 GB'
  cpus 1
  input:
  tuple val(sample_id), file(bam) from bam_ch

  output:
  tuple val(sample_id), file("${sample_id}.sorted.bam") into sorted_bam_ch

  script:
  """
  echo samtools sort ${sample_id}.bam
  touch ${sample_id}.sorted.bam
  """
}

It's producing the following result:

enter image description here

As you can see, the two processes are running in parallel (2 of 2)

It could be that your are limited by the ressources? If you PC doesn't have 40GB of memory and 32 CPUs, he will run them sequentially.

EDIT: By adding sleep 10 to mapToRef, we can clearly see both of them are running in parallel now:

enter image description here

EDIT 2: Added condition to "sleep" on one of the samples:

enter image description here

Final EDIT: be mindful of the ressources you allocate to your processes. You could use ${task.cpus} and ${task.memory.toGiga()} to facilitate the configuration of your processes so that they are not different from what you have defined and what is actually used. See CPUs

$\endgroup$
10
  • $\begingroup$ Thanks @Law. My execution screenshot also looks almost exactly the same. Can you replace touch ... on your script with sleep 5 to emulate long running tasks and confirm if the processes are running parallelly? I have well over 20GB memory and 16 CPUs on this machine, so limited resource should not be an issue. $\endgroup$
    – PPrasai
    Sep 4 at 19:47
  • $\begingroup$ @PPrasai, check again my edit, i was doing it ;) $\endgroup$
    – user324810
    Sep 4 at 19:49
  • 1
    $\begingroup$ Check my second edit. I have added a condition to mapToRef to "sleep" if one of the sample_id match one of the samples. As you can see, both of them runs in parallel first, since one of the samples doesn't match the condition, the other processes are executed afterward. $\endgroup$
    – user324810
    Sep 4 at 19:58
  • 1
    $\begingroup$ Thanks @Law! Reducing cpus to 4 on process and mapToRef, and memory to 4 GB make it all work parallelly! Thank you so much. Could you update your answer so that I can approve it $\endgroup$
    – PPrasai
    Sep 4 at 20:11
  • 1
    $\begingroup$ Well, I did say toward the end the following It could be that your are limited by the ressources? If you PC doesn't have 40GB of memory and 32 CPUs, he will run them sequentially. Also, you're welcome. $\endgroup$
    – user324810
    Sep 4 at 20:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.