2
$\begingroup$

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.

$\endgroup$
5
  • 1
    $\begingroup$ Oh yes, I do upvote every response. and after implementing I accept one of the solutions. $\endgroup$
    – pubsurfted
    Commented Dec 23, 2022 at 8:27
  • 1
    $\begingroup$ Thanks @pubserfted and Merry Christmas! $\endgroup$
    – M__
    Commented Dec 23, 2022 at 13:38
  • 1
    $\begingroup$ @M__ Merry Christmas! Hope you have a good one. ^-^ $\endgroup$
    – pubsurfted
    Commented Dec 23, 2022 at 14:02
  • 1
    $\begingroup$ Thank you @pubsurfted! I'm worked to the bone and will definitely enjoy the vacation. $\endgroup$
    – M__
    Commented Dec 23, 2022 at 14:44
  • $\begingroup$ Having to generate a sample sheet is just wasting time when you can detect it from the directory contents. $\endgroup$
    – raygozag
    Commented Apr 2, 2023 at 1:00

2 Answers 2

2
$\begingroup$

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]
$\endgroup$
2
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Welcome to the site @ATpoint $\endgroup$
    – M__
    Commented Dec 22, 2022 at 22:07

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.