2
$\begingroup$

I have a pipeline that has an input list of samples. I would like to split those samples by chromosomes, process each chromosome apart, and then merge them back by sample name.

Here's an example:

samplesList = Channel.from(["sample1", "sample2"])

// A process that splits the samples by chromosome as an example
// In my actual case, I have a VCF that is split by `snpSift Split` command
// and generates files like [sample]_[chr].vcf
// In the example below, I only gave sample1
// But take into consideration that I have multiple samples

process bar {
   echo true

  output:
    file("*.tsv") into bar_output

  script:
  """
  echo 1 > sample1_chr1.tsv
  echo 2 > sample1_chr2.tsv
  echo 3 > sample1_chr3.tsv
  echo 4 > sample1_chr4.tsv
  echo 5 > sample1_chr5.tsv
  """
}

// Process each chromosome apart for each sample
// This means the output of the process foo
// is only ONE chromosome per sample
process foo {
  input:
    each x from bar_output  // to treat each file separately from the channel
  output:
    file("*.tsv") into xasxa
  script:
    """
    echo \\\$(cat $x)\\\$(cat $x) > sample1_chr\\\$(cat $x)_2.tsv
    """
}

// How should I write the next process to group back them together by sample name?
// I do not wish to do a collect() on channel foo because it will force 
// all the processes of foo to finish first.
// In this example, foo would emit 5 different channels of chromosome for sample1
// This process should merge sample1 only when all chromosomes are processed
// independently of other samples. 
process mergeChromosomesBySamples {
   input:
     ???
}

As the example above, the process mergeChromosomesBySamples should return something similar/close to

["sample1", ["/path/to/workdir/sample1_chr1.tsv",
             "/path/to/workdir/sample1_chr2.tsv",
             "/path/to/workdir/sample1_chr3.tsv",
             "/path/to/workdir/sample1_chr4.tsv",
             "/path/to/workdir/sample1_chr5.tsv"]
]

EDIT1: more explication for the process mergeChromosomesBySamples is given in the comments and as follows ...

The channel xasxa processes each chromosome apart, this means that the output of the channel xasxa contains only one chromosome of the sample like "/path/to/workdir/sample1_chr1.tsv" only. My goal is then to merge the processes out xasxa for each sample together.

I do not wish to use collect() in this case because it would have to force all the processes of foo to finish before. I want them to work independently for each sample.

Where each sample has its own VCF generated and I process them independently.

EDIT2: Here's a graph of what I am trying to achieve:

enter image description here

If I use .collect(), then all the processing on the chromosomes of all the samples has to finish first before being able to merge them, which is something I do not want.

I appreciate your help.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
5
$\begingroup$

Use the groupBy operator. If the content of channel xasxa is:

["/path/to/workdir/sample1_chr1.tsv", "/path/to/workdir/sample1_chr2.tsv", "/path/to/workdir/sample2_chr1.tsv", "/path/to/workdir/sample2_chr2.tsv

You could probably do:

xasxa.groupBy{ it -> it.fileName.toString().split("_")[0] }

The result would be a channel with:

[sample1:[/path/to/workdir/sample1_chr1.tsv, /path/to/workdir/sample1_chr2.tsv, sample2:[/path/to/workdir/sample2_chr1.tsv, /path/to/workdir/sample2_chr1.tsv]]

This breaks if you have an underscore in your samplename before _chr*. More robust would be to initialise all channels with a sample value, so that each channel is in the form of [samplename, chromosome]. Then all you need is a groupTuple on your channel.

Improve this answer
$\endgroup$
5
  • $\begingroup$ The channel xasxa processes each chromosome apart, this means that the output of the channel xasxa contains only one chromosome of the sample like "/path/to/workdir/sample1_chr1.tsv" only. My goal is then to merge the processes out xasxa for each sample together. Sorry if that wasn't clear in my explanation. $\endgroup$
    – user324810
    Dec 16, 2020 at 17:07
  • $\begingroup$ Still not exactly sure what you need, does a .collect() on xasxa do what you wish? $\endgroup$
    – Pallie
    Dec 16, 2020 at 17:11
  • $\begingroup$ I do not wish to use collect() in this case because it would have to force all the processes of foo to finish before. I want them to work independently for each sample. $\endgroup$
    – user324810
    Dec 16, 2020 at 17:13
  • $\begingroup$ Do you start with a single vcf that contains multiple samples? $\endgroup$
    – Pallie
    Dec 17, 2020 at 14:29
  • $\begingroup$ No, each sample has his own VCF generated and I process them independantly. BTW, thanks for the proposed solution. $\endgroup$
    – user324810
    Dec 17, 2020 at 15:11
2
$\begingroup$

I was able to find a solution that answers my problem finally.

So, here are two different methods below, one which does exactly what I want while the other has to wait for all the previous processes to finish:

First method: waiting for previous process to finish using .groupTuple() like what @Pallie mentionned. This method is simpler than the second method. The problem is that it has to wait for the previous processes to finish. I want my samples to be treated in parallel. If I missed something with the logic of the method, please feel free to edit.

#!/usr/bin/env nextflow

// my samples
samplesList = Channel.from(["sample1", "sample2"])

// the process that is supposed to split the sample by chromosomes
process foo {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
         val(x) from samplesList
      output:
         file("*.vcf") into foo_out
         val(x) into output_sample
      script:
      """
         echo 1 > ${x}_chr1.vcf
         echo 2 > ${x}_chr2.vcf
      """
}

//Combine the samples by chromosome and filter out only the samples that match the files output
foo_out.flatten().combine(output_sample).filter{ it[0].toString().contains(it[1].toString()) }.set{maybework}

// the process that is supposed to treat each chromosome apart
process bar {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
         set file(y), val(x) from maybework
      output:
         set val(x), file("*chr*vcf") into bar_out
      script:
      """
         echo ${x}
         filename=\\\$(basename -- ${y})
         chromosome=\\\$(echo "\\\${filename%.*}" | cut -f2 -d'_')
         cat ${y} > ${x}_"\\\$chromosome"_"\\\$chromosome".vcf
         if [ ${x} == "sample2" ]; then sleep 20; fi
      """
}

// the process that merges back the chromosomes by sample
process mergefoobar {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
      set val(sample), file(x) from bar_out.groupTuple()

      script:
      """
         echo $sample
         echo ${x}
      """
}

will output: FirstMethod As you can see above, the process mergefoobar is still waiting for all the processes of bar to finish.

Second method: does not wait for the previous processes to finish, uses .join() for each chromosome... Yes, the downside is that you would have to write the join() for all the chromosomes. Additionally, if you use snpSift Split, one sample may generate 3 files out of 5 (because only 3 chromosomes had variants listed in them), you would need to create the other empty VCF files because the output excepts them.

#!/usr/bin/env nextflow

// my samples
samplesList = Channel.from(["sample1", "sample2"])

// the process that is supposed to split the sample by chromosomes
process foo {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
         val(x) from samplesList
      output:
         file("*.vcf") into foo_out
         val(x) into output_sample
      script:
      """
         echo 1 > ${x}_chr1.vcf
     echo 2 > ${x}_chr2.vcf
         echo X > ${x}_chrX.vcf
     echo Y > ${x}_chrY.vcf
      """
}

//Combine the samples by chromosome and filter out only the samples that match the files output
foo_out.flatten().combine(output_sample).filter{ it[0].toString().contains(it[1].toString()) }.set{maybework}

// the process that is supposed to treat each chromosome apart
process bar {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
         set file(y), val(x) from maybework
      output:
         set val(x), file("*chr1*vcf") optional true into bar_out1
         set val(x), file("*chr2*vcf") optional true into bar_out2
         set val(x), file("*chrX*vcf") optional true into bar_outX
         set val(x), file("*chrY*vcf") optional true into bar_outY
      script:
      """
         echo ${x}
     filename=\\\$(basename -- ${y})
         chromosome=\\\$(echo "\\\${filename%.*}" | cut -f2 -d'_')
         cat ${y} > ${x}_"\\\$chromosome"_"\\\$chromosome".vcf
         if [ ${x} == "sample2" ]; then sleep 20; fi
      """
}

// the process that merges back the chromosomes by sample
process mergefoobar {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
      set val(sample), file(x), file(y), file(i), file(t) from bar_out1.join(bar_out2).join(bar_outX).join(bar_outY)

      script:
      """
         echo $sample
     echo ${x}
         echo ${y}
     echo ${i}
         echo ${t}
      """
}

SecondMethod Above, you see that the process mergefoobar has launched even if the previous bar has not finished yet.

I hope this also helps anyone in the future who might have this issue.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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