How does one account for both single end and paired end reads as input in a Nextflow pipeline?

Question

Suppose you have a pipeline that takes reads as input from the user and trims them.

// Declare syntax version
nextflow.enable.dsl=2

/* 
 * pipeline input parameters 
 */


params.reads = "$projectDir/data/*_{1,2}.fastq"


/*
 * define the `TRIM` process that creates trimmed.reads.fastq files
 * given the reads.fastq file. 
 */
 
process TRIM {

    tag "$pair_id"

    input:
    tuple val(pair_id), path(reads)

    output:
    path(pair_id)

    script:
    """
    atria -r ${reads[0]} -R ${reads[1]} -t 1 -l 12 -q 30 -o $pair_id
    """
}

workflow {

    read_pairs_ch = Channel.fromFilePairs( params.reads, checkIfExists:true )

    trim_ch=TRIM(read_pairs_ch)

}

This particular one takes paired end reads as input. What if the user wanted to trim both single end and paired end reads, how to account for single end reads as the input?

Best wishes.

Answer 1

There are many ways to accomplish this. Assuming you don't need to process single-end and paired-end reads simultaneously (i.e. within the same run), one way using the fromFilePairs factory method is to just use an additional flag to define the number of files each emitted item is expected to hold. And then adjust your process definition accordingly:

params.reads = "${projectDir}/data/*_{1,2}.fastq.gz"
params.singleEnd = false


process TRIM {

    tag "${pair_id}"

    cpus 1

    input:
    tuple val(pair_id), path(reads)

    output:
    tuple val(pair_id), path("${pair_id}/*.atria.fastq{,.gz}")

    script:
    def single = reads instanceof Path

    def read1 = !single ? /-r "${reads[0]}"/ : /-r "${reads}"/
    def read2 = !single ? /-R "${reads[1]}"/ : ''

    """
    atria \\
        -t ${task.cpus} \\
        ${read1} \\
        ${read2} \\
        -l 12 \\
        -q 30 \\
        -o "${pair_id}"
    """
}

workflow {

    Channel
        .fromFilePairs( params.reads, size: params.singleEnd ? 1 : 2 )
        .ifEmpty { error "Could not find any reads matching pattern: ${params.reads}" }
        .set { reads }

    TRIM( reads ).view()
}

Although we've used an additional flag, it does let the pipeline fail immediately if it can't find the files it expects. I think this approach is more explicit and provides a better user experience:

$ nextflow run main.nf --reads './data/*.fastq'
N E X T F L O W  ~  version 22.10.0
Launching `main.nf` [nasty_allen] DSL2 - revision: 2cbc311974
Could not find any reads matching pattern: ./data/*.fastq

 -- Check script 'main.nf' at line: 38 or see '.nextflow.log' file for more details

Doh! I forgot that these are single ended reads:

$ nextflow run main.nf --reads './data/*.fastq' --singleEnd
N E X T F L O W  ~  version 22.10.0
Launching `main.nf` [clever_gutenberg] DSL2 - revision: 2cbc311974
executor >  local (3)
[8e/f90c90] process > TRIM (X) [100%] 3 of 3 ✔
[Y, /path/to/work/e6/48df1485bb7de8355ed8ce0f1b53da/Y/Y.atria.fastq]
[Z, /path/to/work/eb/aca3389344ac03ff283ded0c46a7a6/Z/Z.atria.fastq]
[X, /path/to/work/8e/f90c90210103c03a3fbde697e01227/X/X.atria.fastq]

Answer 2

My personal strategy is to always read data from a samplesheet which then auto-detects whether a sample is paired-end or single-end and stores this information in a so-called meta map (inspired by nf-core) which then passes this information to the modules.

Example:

Say you have a samplesheet that is a three-column csv file with a header:

$ cat samplesheet.csv 
sample,r1,r2
sample1,/Users/atpoint/sample1_R1.fq.gz,/Users/atpoint/sample1_R2.fq.gz
sample2,/Users/atpoint/sample2_R1.fq.gz,

Read this into Nextflow and detect whether r2 has an entry. If so it is paired-end, else it is single-end:

# main.nf

nextflow.enable.dsl = 2

// Channel for the samplesheet
ch_samplesheet = Channel.fromPath("/Users/atpoint/samplesheet.csv")

// Parse it line by line
ch_reads = ch_samplesheet.splitCsv(header:true).map {

    // This is the read1 and read2 entry
    r1 = it['r1']
    r2 = it['r2']

    // Detect wiether single-end or paired-end
    is_singleEnd = r2.toString()=='' ? true : false
    
    // The "meta" map, which is a Nextflow/Groovy map with id (the sample name) and a single_end logical entry
    meta = [id: it['sample'], single_end: is_singleEnd]
    
    // We return a nested map, the first entry is the meta map, the second one is the read(s)
    r2.toString()=='' ? [meta, [r1]] : [meta, [r1, r2]]

}

If you .view() the ch_reads you get:

[[id:sample1, single_end:false], [/Users/atpoint/sample1_R1.fq.gz, /Users/atpoint/sample1_R2.fq.gz]]
[[id:sample2, single_end:true], [/Users/atpoint/sample2_R1.fq.gz]]

A process that can use that could be:

process A {

    tag "$meta.id"

    input:
    tuple val(meta), path(reads)

    script:
    
    // The meta map holds the single/paired information that can now easily be parsed in the module in a standardized way
    if(meta.single_end){

        """
        do_single_end_things with ${reads}
        """

    } else {

        """
        do_paired_end_things ${reads}[0] ${reads}[1] 
        """
        
    }

}

A(ch_reads)

The meta map idea (cudos nf-core) is simple, yet versatile and powerful as it allows standardized parsing of such metadata-like information (single_end, libtype, id/sample_name, any other value you can think of) and most importantly, it allows to process single-end paired end data at the same time. All the decision making is based on the meta map and therefore automated if the modules are coded appropriately as above.

Hope that helps.