Nextflow: Not all items in channel used by process

Viewed 282

I've been struggling to identify why a nextflow (v20.10.00) process is not using all the items in a channel. I want the process to run for each sample bam file (10 in total) and for each chromosome (3 in total).

Here is the creation of the channels and the process:

ref_genome = file( params.RefGen, checkIfExists: true )
ref_dir    = ref_genome.getParent()
ref_name   = ref_genome.getBaseName()
ref_dict   = file( "${ref_dir}/${ref_name}.dict", checkIfExists: true )
ref_index  = file( "${ref_dir}/${ref_name}.*.fai", checkIfExists: true )

// Handles reading in data if the previous step is skipped
if( params.Skip_BP ){
  Channel
    .fromFilePairs("${params.ProcBamDir}/*{bam,bai}") { file -> file.name.replaceAll(/.bam|.bai$/,'') }
    .ifEmpty { error "No bams found in ${params.ProcBamDir}" }
    .map { ID, files -> tuple(ID, files[0], files[1]) }
    .set { processed_bams }
}
// Setting up the chromosome channel
if( params.Chroms == "" ){
  // Defaulting to using all chromosomes
  chromosomes_ch = Channel
                      .from("AgamP4_2L", "AgamP4_2R", "AgamP4_3L", "AgamP4_3R", "AgamP4_X", "AgamP4_Y_unplaced", "AgamP4_UNKN")
  println "No chromosomes specified, using all major chromosomes: AgamP4_2L, AgamP4_2R, AgamP4_3L, AgamP4_3R, AgamP4_X, AgamP4_Y_unplaced, AgamP4_UNKN"
} else {
  // User option to choose which chromosome will be used
  // This worked with the following syntax nextflow run testing.nf --profile imperial --Chroms "AgamP4_3R,AgamP4_2L"
  chrs = params.Chroms.split(",")
  chromosomes_ch = Channel
                    .from( chrs )
  println "User defined chromosomes set: ${params.Chroms}"
}


process DNA_HCG {
  errorStrategy { sleep(Math.pow(2, task.attempt) * 600 as long); return 'retry' }
  maxRetries 3
  maxForks params.HCG_Forks

  tag { SampleID+"-"+chrom }

  executor = 'pbspro'
  clusterOptions = "-lselect=1:ncpus=${params.HCG_threads}:mem=${params.HCG_memory}gb:mpiprocs=1:ompthreads=${params.HCG_threads} -lwalltime=${params.HCG_walltime}:00:00"

  publishDir(
    path: "${params.HCDir}",
    mode: 'copy',
  )

  input:
  each chrom from chromosomes_ch
  set SampleID, path(bam), path(bai) from processed_bams
  path ref_genome
  path ref_dict
  path ref_index

  output:
  tuple chrom, path("${SampleID}-${chrom}.vcf") into HCG_ch
  path("${SampleID}-${chrom}.vcf.idx") into idx_ch
  
  beforeScript 'module load anaconda3/personal; source activate NF_GATK'

  script:
  """
  if [ ! -d tmp ]; then mkdir tmp; fi
  taskset -c 0-${params.HCG_threads} gatk --java-options \"-Xmx${params.HCG_memory}G -XX:+UseParallelGC -XX:ParallelGCThreads=${params.HCG_threads}\" HaplotypeCaller \\
    --tmp-dir tmp/ \\
    --pair-hmm-implementation AVX_LOGLESS_CACHING_OMP \\
    --native-pair-hmm-threads ${params.HCG_threads} \\
    -ERC GVCF \\
    -L ${chrom} \\
    -R ${ref_genome} \\
    -I ${bam} \\
    -O ${SampleID}-${chrom}.vcf ${params.GVCF_args}
  """
}

But for reasons I cannot figure out, nextflow only creates 3 jobs: [d8/45499b] process > DNA_HCG (0_wt5_BP-CM029350.1) [ 0%] 0 of 3

I thought maybe it was because it only took the first sample and then one process for each chromosome. Though I doubted this since the code works for a different reference genome correctly. Regardless, I adjusted the input channels:

processed_bams
  .combine(chromosomes_ch)
  .set { HCG_in }

and

input:
set SampleID, path(bam), path(bai), chrom from HCG_in

But this resulted in only a single job being created: [6e/78b070] process > DNA_HCG (0_wt10_BP-CM029350.1) [ 0%] 0 of 1

Confusingly, when i use HCG_in.view() there are 30 items. And to further confuse me the correct number of jobs comes from the following code:

chrs = params.Chroms.split(",")
      chromosomes_ch = Channel
                         .from(chrs)

Channel
  .fromFilePairs("${params.ProcBamDir}/*{bam,bai}") { file -> file.name.replaceAll(/.bam|.bai$/,'') }
  .ifEmpty { error "No bams found in ${params.ProcBamDir}" }
  .map { ID, files -> tuple(ID, files[0], files[1]) }
  .set { processed_bams }

process HCG {
  executor 'local'

  input:
  each chrom from chromosomes_ch
  set SampleID, path(bam), path(bai) from processed_bams
  //set SampleID, path(bam), path(bai), chrom from HCG_in

  script:
  """
  echo "${SampleID} - ${chrom}"
  """
}

Output: [75/c1c25a] process > HCG (27) [100%] 30 of 30 ✔

I'm hoping I've just missed something obvious, but I cannot see it at the moment. Thanks in advance for the help.

1 Answers

Issues like this almost always involve the use of multiple input channels:

When two or more channels are declared as process inputs, the process stops until there’s a complete input configuration ie. it receives an input value from all the channels declared as input.

Your initial assessment was correct. However, the reason only three processes were run (i.e. one sample for each of the three chromosomes), is because this line (probably) returned a list (i.e. a java LinkedList) containing a single element, and lists behave like queue channels:

ref_index  = file( "${ref_dir}/${ref_name}.*.fai", checkIfExists: true )

You might have expected this to return a UnixPath. Ultimately, the solution is to ensure ref_index is value channel.

Related