This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How should I define output for a Nextflow demultiplexing rule?

I'm trying to write a Nextflow rule which takes in a PacBio HiFi reads BAM and demultiplexes it using the tool lima ( https://lima.how ).

The input will be a single HiFi reads BAM, a FASTA containing barcode sequences, and a sample sheet:

Barcode, Bio Sample ...
bc2001--bc2001,BioSample1
bc2002--bc2002,BioSample2
...
bc200n--bc200n,BioSampleN

I'm familiar with the lima tool and I know how I want the command to look:

lima \
  --output-missing-pairs \
  --hifi-preset SYMMETRIC-ADAPTERS \
  --store-unbarcoded \
  --split-named \
  --biosample-csv ${BIOSAMPLE_CSV} \
  ${INPUT_BAM} \
  ${BARCODE_FASTA} \
  ${OUTPUT_LOCATION}

The output of this will be a collection of demultiplexed HiFi reads BAM files, one for each barcode plus another for any unassigned reads, as well as index files and XML files for each, e.g.:

redo.hifi_reads.bc2001--bc2001.bam
redo.hifi_reads.bc2001--bc2001.bam.pbi
redo.hifi_reads.bc2001--bc2001.consensusreadset.xml
redo.hifi_reads.bc2002--bc2002.bam
redo.hifi_reads.bc2002--bc2002.bam.pbi
redo.hifi_reads.bc2002--bc2002.consensusreadset.xml
redo.hifi_reads.bc200n--bc200n.bam
redo.hifi_reads.bc200n--bc200n.bam.pbi
redo.hifi_reads.bc200n--bc200n.consensusreadset.xml
redo.hifi_reads.unbarcoded.bam
redo.hifi_reads.unbarcoded.bam.pbi
redo.hifi_reads.unbarcoded.consensusreadset.xml

I have written the rule successfully such that it demultiplexes a given input BAM. I can find the results of this is in the relevant work directory.

Rule definition

I am struggling to define what output Nextflow should provide from my lima rule. I have tried to write this rule and I include it here:

process limaNew {
    conda 'bioconda::lima==2.12'

    publishDir 'results', mode: 'symlink'

    input:
        val barcodes
        path undoneBAM
        path sampleSheet
        path barcodeFASTA

    //output:

    script:
    """
    lima -j32 --output-missing-pairs \
        --hifi-preset SYMMETRIC-ADAPTERS \
        --store-unbarcoded \
        --split-named \
        --biosample-csv "$sampleSheet" \
        "$undoneBAM" \
        "$barcodeFASTA" \
        redo.hifi_reads.bam
    """
}

As you can see, the output definition is blank. I tried to do this a certain way which perhaps overcomplicated things. I wrote this channel in the workflow definition which successfully generates a list of barcode names from the input samplesheet:

    barcodeList = channel.fromPath(params.sampleSheet)
        .splitCsv(header:true)
        .map { row-> row.Barcode.replaceFirst(/--.*/, '') }
        .collect()

However, trying to then use this barcodeList in the output definition kept throwing errors saying that the variable does not exist. Apparently it is not allowable to use such a variable in the output definition?

Regardless, I think trying to solve this is not useful, I believe my whole approach is wrong. I've read back through documentation and tutorials but I can't find anything which seems to explain how you handle this more unpredictable form of output. I looked for analogous rules in other people's workflows but haven't really found anything I can generalise to my problem.

The below does work:

    output:
        file('*.bam')
        file('*.pbi')

Is this sufficient/recommended?

I'm currently doing as many tasks as I can using Nextflow to try to get to grips the it, I really want to learn as much as I can here. But try as I might, this particular problem eludes me and I'd appreciate advice on what is philosophically the most Nextlfow way to do this. Other than that, I apologise for the overlong post and I thank you for reading if you got this far!

hifi lima demultiplexing nextflow
• 380 views • link
written by sholt6 • sholt6

1 answer 

output:
        path('*.bam'), emit: bams
        path('*.pbi'), emit: pbis

I think your approach sounds fine if you'd like to do it like this. I'd change the deprecated file to path, and use the emit tag. The emit helps to subdivide unrelated output channels.

Then you can use the files if you need to with a subsequent process workflow description like

process2(process1.out.bams)
• 0 views • link
written by colindaven colindaven

Thanks for your reply, appreciate the improvements. Is it standard practice to generalise outputs in such a way? My original approach was to try to anticipate exactly what output filenames will be produced and I wasn't able to find a way to do that.

• 0 views • link
written by sholt6 • sholt6

I guess that depends on the program. If you are writing the outputs yourself you can decide on the naming. If lima decides dynamically (depending on number of barcodes or whatever), you would seem to be stuck.

But thats ok, these are output channels. So if there is one bam output, your process2 (say calculating bam stats) only runs once. If 10 bams, ten (parallel) runs of process2 would be started.

• 0 views • link
written by colindaven colindaven

This is great, sure enough they do it with a simple wildcard:

    output:
    // Demultiplexed file name changes depending on the arg '--output-types'
    tuple val(meta), path('output/*.fq.gz')                         , emit: sample_fastq

So I can safely say I was overthinking it originally! And I guess I can still generate a list of sample names/barcode IDs from the TSV if I want to be explicit.

• 0 views • link
written by sholt6 • sholt6

Log in to answer this question.

More posts like this