2
$\begingroup$

I'd like to merge files from different lanes for forward and reverse reads. For e.g., if a sample A has three lanes in forward, and three reverse, ultimately one reverse and one forward is to be generated. I have normal and tumor for each sample. N as normal and T as tumor.
The lanes files are for e.g. as:

MM-265-DNA-N-01-01_L003_R1_001.fastq.gz
MM-265-DNA-N-01-01_L003_R2_001.fastq.gz
MM-265-DNA-N-01-01_L004_R1_001.fastq.gz
MM-265-DNA-N-01-01_L004_R2_001.fastq.gz

MM-265-DNA-T-01-01_L003_R1_001.fastq.gz
MM-265-DNA-T-01-01_L003_R2_001.fastq.gz
MM-265-DNA-T-01-01_L004_R1_001.fastq.gz
MM-265-DNA-T-01-01_L004_R2_001.fastq.gz

MM-465-DNA-N-01-01_L003_R1_001.fastq.gz
MM-465-DNA-N-01-01_L003_R2_001.fastq.gz
MM-465-DNA-N-01-01_L004_R1_001.fastq.gz
MM-465-DNA-N-01-01_L004_R2_001.fastq.gz

MM-465-DNA-T-01-01_L003_R1_001.fastq.gz
MM-465-DNA-T-01-01_L003_R2_001.fastq.gz
MM-465-DNA-T-01-01_L004_R1_001.fastq.gz
MM-465-DNA-T-01-01_L004_R2_001.fastq.gz  

The final output needs to be as:

MM-465_normal_R1.fastq MM-465_normal_R2.fastq
MM-465_tumor_R1.fastq MM-465_tumor_R2.fastq
MM-265_normal_R1.fastq MM-265_normal_R2.fastq
MM-265_tumor_R1.fastq MM-265_tumor_R2.fastq

I've code below broken in three parts:

  1. regular expression
  2. process that iterates over files
  3. workflow that calls the above two code blocks

regex:

def get_sample_id( file ) {
def matcher = (file =~ /(MM-\d+)-DNA-N-.+/)
    if (matcher) {
        return matcher[0][1]
    } else {
        return ''
    }
}

process:

process file_print {

debug true

input:

tuple val(sample_id), path(datasetFile)
val(sample_status) 
val(read_value)

output:
publishDir '/arion/learn_new_concepts/', mode:'copy'

tuple val (sample_id), path("${output_file_name}")

script:

output_file_name = "${sample_id}_${sample_status}_${read_value}.fastq"

for (g in datasetFile){
println(g)

"zcat ${g} >> ${output_file_name} "

}

}

Workflow:

workflow {

normal_reads=Channel.fromFilePairs("/sc/arion/MM-*-N-*_L00*_R{1,2}_*.fastq.gz")
//MM-2658-DNA-N-01-01_L004_R2_001.fastq.gz  //

normal_forward_reads=normal_reads
        .map{ file -> return tuple( get_sample_id(file), file[1][0] )}.groupTuple(sort:true ) 

file_print(normal_forward_reads,"normal","R1").view() /// similar for reverse [1][1], R2 .. for tumor pass reads, "tumor", "R1"/R2"

}

I have not put code for reverse, but I get error as:

Caused by:   Process `file_print (1)` terminated with an error exit status (127) Command executed:   null

Command exit status:   127

I tried many things such as ${ output_file_name }, { output_file_name }, $output_file_name but in vain.

I'd like to improve few things in code:

a) The iteration/loop over file in datasetFile I think something like each could help, but I need order to be maintained

b) Order in which loop of datasetFile is created. For instance, I'd like to keep lanes L003, L004, L005 in this order. Or, L001, L002, L003 order, I hope this is clear.

I cannot move ahead or able to sort this by myself.

I do not quite understand how and when to use .map{}

Edited: after output added tuple for val and path. I still get same error.

$\endgroup$

1 Answer 1

2
$\begingroup$

I think you are on the right track using groupTuple(sort: true). Another way to get what you want would be to instead parse the fromFilePairs grouping key. Unless you're looking to do more than just concatenate (and decompress) FASTQ files, I think all the process really requires is an output prefix (in addition to the FASTQ files of course), for example:

params.reads = "/sc/arion/MM-*_R{1,2}_001.fastq.gz"


process concat_reads {

    tag { prefix }

    input:
    tuple val(prefix), path(reads, stageAs: 'reads/*')

    output:
    tuple val(prefix), path("${prefix}.fastq")

    """
    zcat -f ${reads} > "${prefix}.fastq"
    """
}
workflow {

    Channel.fromFilePairs( params.reads ) \
        | branch { key, fastqs ->

            def matcher = key =~ /^(MM-\w+)-DNA-(\w+)-.*_L(\d+).*$/

            matches: matcher
                def (_, sample_id, sample_type, lane) = matcher[0].findAll()

                def new_sample_id = "${sample_id}_${sample_type}"

                return tuple( new_sample_id, lane as int, *fastqs )

            other: true
                return tuple( key, *fastqs )
        } \
        | set { samples }

    samples.matches \
        | groupTuple( sort: true ) \
        | flatMap { sample_id, lanes, fq1s, fq2s ->
            [
                tuple( "${sample_id}_R1", fq1s ),
                tuple( "${sample_id}_R2", fq2s ),
            ]
        } \
        | concat_reads \
        | branch { prefix, fastq ->

            def idx = prefix.lastIndexOf('_')

            def sample_id = prefix.substring(0, idx)
            def suffix = prefix.substring(idx + 1)

            read1: suffix == 'R1'
                return tuple( sample_id, fastq )
            read2: suffix == 'R2'
                return tuple( sample_id, fastq )
        } \
        | set { result }

    result.read1 \
        | join( result.read2 ) \
        | map { sample_id, fastq1, fastq2 ->

            tuple( sample_id, [ fastq1, fastq2 ] )
        } \
        | view()
}

Results:

$ nextflow run main.nf
N E X T F L O W  ~  version 23.04.1
Launching `main.nf` [jolly_noether] DSL2 - revision: 4a81f7a875
executor >  local (8)
[b8/a0b534] process > concat_reads (MM-265_T_R2) [100%] 8 of 8 ✔
[MM-465_T, [/path/to/work/9c/8b6fd328c848379423a0aba098cb74/MM-465_T_R1.fastq, /path/to/work/34/913780b1da902b5035f3bf56aac038/MM-465_T_R2.fastq]]
[MM-465_N, [/path/to/work/69/595120a1ee0a42217d06855ca9d1bc/MM-465_N_R1.fastq, /path/to/work/8f/31c406c5180073de577502fdb88fbf/MM-465_N_R2.fastq]]
[MM-265_N, [/path/to/work/1b/9c17cefb5aa51b3ee69832dc65d915/MM-265_N_R1.fastq, /path/to/work/d0/aac9f46e4a7f7165693b2537ad5630/MM-265_N_R2.fastq]]
[MM-265_T, [/path/to/work/d5/6295ff109df3ccb829a6f0a97b6e45/MM-265_T_R1.fastq, /path/to/work/b8/a0b53499ca13a523d811ecf772449b/MM-265_T_R2.fastq]]

Note that unless the FASTQ files are small, I would try to avoid decompressing the FASTQ files if possible. Most bioinformatics tools that deal with large FASTQ files will likely have gzip support, meaning you can replace the zcat -f with a cat in the example above (and remember to add the .gz extension to the output filename).

$\endgroup$
10
  • $\begingroup$ Thank you for your reply and code. I'm sorry, I cannot understand what you did in the code as it is scary and complicated for me. May I know what is wrong in my code? I've .fastq extension as a suffix. May be a little explanation of the steps you did would be highly helpful to me. For example, when to branch branch. Or, how do I put normal in the concatenated files? $\endgroup$ Jul 20, 2023 at 11:42
  • 1
    $\begingroup$ @DeathMetal I suspect your error message is being truncated because it tells us exactly why it failed: .command.sh: line 2: null: command not found. Try running your workflow instead using nextflow run main.nf -ansi-log false. Since your script block returns null, this becomes the command Nextflow tries to run for better or worse. I think you could replace this with a script block containing: zcat -f ${datasetFile} > "${output_file_name}" to get it running. Also, use normal_reads.map{ key, files -> to unpack the fromFilePairs output. $\endgroup$
    – Steve
    Jul 20, 2023 at 12:55
  • 1
    $\begingroup$ @DeathMetal Note that processes do not necessarily process items in the order that they are received. For example, concat_reads can return results in any order. So, if we want to get back the reads where the R1 and R2 files are the first and second in the pair, respectively, then we need some way to separate the outputs so they can be joined correctly. In the above, I just used branch to create two channels: read1 and read2 and then a join to put them back together. $\endgroup$
    – Steve
    Jul 20, 2023 at 13:06
  • 1
    $\begingroup$ @DeathMetal Also, if you wanted to replace N with normal and T with tumor, I think the place to do this would when defining the new_sample_id in the above. You could use a Map for example. $\endgroup$
    – Steve
    Jul 20, 2023 at 13:15
  • 1
    $\begingroup$ @DeathMetal If you just need the output of concat_reads, the code after and including the last branch statement is not necessary. But usually you would want to do something with these new pairs of files, which is why I included some code to show how this could be done. In reality, there are lots of ways to do this though. The underscore is just a placeholder for the full match. In Python and Groovy multi-assignment, at least, it is convention is to use an underscore for un-used variables. Use samples.matches.view() or insert a println(your_variable) into the branch statement to debug $\endgroup$
    – Steve
    Jul 21, 2023 at 0:16

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.