4
$\begingroup$

I have a workflow that needs to:

  1. Generate .fastq files from .bams (while preserving the @RG group from the original bam)
  2. Split the .fastqs
  3. Align
  4. etc.

The only way I could think to preserve the @RG group was to print it to a file during step 1 and emit the file. I can't figure out how to synchronize the emitted @RG file with the .fastqs, however. I've tried several variations without luck.

Questions:

  1. How do I synchronize the emitted @RG file (sample_RG_file) with the split .fastqs?
  2. Is there a better way?
  3. Bonus question (may be worth a separate post): does Nextflow not stage files until the actual .bash script starts? I tried to operate on the sample_RG_file within the script using Groovy, but got an error that the file didn't exist. (The code was working before I found I needed to synchronize sample_RG_file with the

Other variations:

  1. The code was 'working' before I found I needed to synchronize sample_RG_file with the split .fastqs. I was originally passing GENERATE_FASTQS_PROC.out.sample_RG as a separate parameter to ALIGN_FASTQ_BWA_PROC, but it was always passing the same (first) file.
  2. I also tried GENERATE_FASTQS_PROC.out.sample_RG.flatten() in the .map operator, but that had the same error.

Current error:

The code below fails with the following error:

Error executing process > 'REALIGN_SAMPLES_WF:ALIGN_FASTQ_BWA_PROC (1)'

Caused by: Not a valid path value type: groovyx.gpars.dataflow.DataflowBroadcast (DataflowBroadcast around DataflowStream[?])

Semi-pseudo code

Workflow

workflow my_WF {
    take:
        sample_input_file_ch
        original_ref

    main:
        GENERATE_FASTQS_PROC(sample_input_file_ch, original_ref)
        
        reads_per_job = 1_000_000
        fq_split_reads_and_RG_ch = GENERATE_FASTQS_PROC.out.fastqs
                        .splitFastq(by: reads_per_job, pe: true, file: true)
                        .map {
                            fq1, fq2 ->
                                tuple( fq1, fq2, GENERATE_FASTQS_PROC.out.sample_RG )
                        }
                        .view()

        ALIGN_FASTQ_BWA_PROC(fq_split_reads_and_RG_ch, ...)
}

GENERATE_FASTQS_PROC

process GENERATE_FASTQS_PROC {

    input:
        path(sample_input_file)
        val(original_ref)

    output:
        tuple path("**/*_R1.fastq"), path("**/*_R2.fastq"), emit: fastqs
        path("**/*_RG.txt"), emit: sample_RG

    script:
    """
    generate_fastqs.sh $sample_input_file $original_ref $task.cpus $mem_per_thread
    """
}

ALIGN_FASTQ_BWA_PROC

process ALIGN_FASTQ_BWA_PROC {

    /*
     * Publish results. 'mode: copy' will copy the files into the publishDir
     * rather than only making links.
     */
    // publishDir("${params.results_dir}/01-REALIGN_SAMPLES", mode: 'copy')

    label 'REALIGN_SAMPLES'

    input:
        tuple path(fq1), path(fq2), path(sample_RG_file)
        val(align_to_ref)
        val(align_to_ref_tag)
        val(original_ref)
        val(output_format)


    output:
        path '*.{bam,cram}', emit: final_alignment

    script:

    /*
     * Tried to read the RG file using Groovy, but
     * got error that the file didn't exist.
     *
     * Does Nextflow not stage files until the actual .bash script starts?
     */
//    println "sample RG file: $sample_RG_file"
//    file_lines = sample_RG_file.readLines()
//    sample_name = file_lines[0]
//    sample_RG = file_lines[1]

    """
    lines=()
    while read line
    do
        lines+=("\\\$line")
    done < $sample_RG_file

    sample_name="\\\${lines[0]}"
    sample_RG="\\\${lines[1]}"

    bash realign_bwa-split.sh $fq1 $fq2 \\\$sample_name \"\\\$sample_RG\" ${align_to_ref} ${output_format} $task.cpus
    """
}

Update

Per @Steve's suggestion, I'm trying to have GENERATE_FASTQS_PROC create the split .fastq files using Unix split, but I'm still struggling with various concepts.

Supposed I have this emitted from my process:

[sample_A_RG.txt, [s1_R1.split_00.fq, s1_R2.split_00.fq, ..., s1_R1.split_XX.fq, s1_R2.split_XX.fq]]`

In my mind, I want to pull out s1_RG.txt file from the Channel and extract the @RG. Then I want to call:

[<all_fastqs_for_sample_A>]
    .map{ 
        // pull out split ID numbers 
    }
    .groupTuple()

I believe this would group _R1 and _R2 for each split set, and then call

[[00.fq1, 00.fq2], ..., [XX.fq1, XX.fq2]]
    .combine(sample_RG_var)

to get [[@RG, 00.fq1, 00.fq2], ..., [@RG, XX.fq1, XX.fq2]].

Seems like it would be so much easier if I could use Channel factories and operators (e.g., .fromFilePairs().combine( rg )) in the output:.

Getting files into a process is relatively straightforward, but getting them out in a reasonable manner is less so.

$\endgroup$

2 Answers 2

3
$\begingroup$

Nextflow channels guarantee that items are delivered in the same order as they are sent. So if a process declares two (or more) output channels, the items that are emitted should already be synchronized. This is true unless of course one of the channels sets the optional output attribute.

  1. How do I synchronize the emitted @RG file (sample_RG_file) with the split .fastqs?

Since the splitFastq operator can create n outputs for each input, you need some way to combine the @RG files using some key. This key, for example, could be the basename of the sample input file. This, I think, would work nicely given the input file should just have a '.bam' extension. Another option, would be to supply some other value in the input declaration and just use that in your output declaration. If we went with the former, your output declaration might look like:

output:
tuple val(sample_input_file.baseName), path("sample_R1.fastq"), path("sample_R2.fastq"), \
    emit: fastqs
tuple val(sample_input_file.baseName), path("sample_RG.txt"), \
    emit: sample_RG

Then, the body of your workflow might look like:

workflow {

    reads = GENERATE_FASTQS_PROC(sample_input_file_ch, original_ref)
        
    reads.fastqs
        .splitFastq( by: params.reads_per_job, pe: true, file: true )
        .combine( reads.sample_RG, by: 0 )
        .map { sample, fq1, fq2, rg_txt -> tuple(fq1, fq2, rg_txt) }
        .set { fq_split_reads_and_RG_ch }

    ...
}

You'll notice that I removed your fuzzy glob patterns from your output declaration above. Your workflow definition suggests that the generate_fastqs.sh script produces a single pair of FASTQ files and a single readgroup file for a given set of inputs. However, as you know, glob patterns can match one or more output files. When declaring output files, it's always best to name these guys explicitly. If you do use a glob pattern and for some reason your process was to produce two files that match a given pattern, you'll instead get back a tuple, which would break the pipeline.

  1. Is there a better way?

Since you need the @RG tags with the split FASTQ files, I don't think there is. However, I generally avoid using the splitting operators on large files. My preference is always to pass "work" off (including IO) to a separate process. A separate process can also provide increased flexibility if/when requirements change.

  1. Does Nextflow not stage files until the actual .bash script starts? I tried to operate on the sample_RG_file within the script using Groovy, but got an error that the file didn't exist.

Nextflow evaluates the code in the script definition before returning the script string to be executed. This is why it's possible to implement a conditional script. Basically, you want to avoid all IO operations here. The script string to be executed has not even been decided at this point.


Update

I'm not sure how you plan to operate per sample. IIUC, [<all_fastqs_for_sample_A>] is basically a Channel of lists, where each list contains all split FASTQ files for a given sample. Using map here to pull out the split ID numbers followed by the groupTuple operator might not give the expected results. You might also run into the problems with calling groupTuple, since, without a size attribute, it will basically wait for all outputs before emitting the grouped result:

You should always specify the number of expected elements in each tuple using the size attribute to allow the groupTuple operator to stream the collected values as soon as possible. However there are use cases in which each tuple has a different size depending on the grouping key. In this case use the built-in function groupKey that allows you to create a special grouping key object to which it’s possible to associate the group size for a given key.

Part of the problem is that we don't know what the generate_fastqs.sh script does and how the other scripts work. If you are aligning using BWA MEM and have samtools available in your $PATH, it might even be easier to produce interleaved output (since BWA MEM can read interleaved FASTQs) using the samtools fastq command. This would save some work having to group each list of FASTQs into pairs, ensuring they are ordered as expected. This would be reasonably involved.

Not sure what BAMs you need to process, or if you're looking for a more generic solution, but I thought I'd throw in an alternative solution you might find works well for you or that you may not have yet considered. Since BAMs often contain multiple readgroups, I find that if I split by readgroup (using samtools split, flatten the output Channel, then samtools collate followed by samtools fastq, I can achieve somewhat uniform mapping walltimes (assuming the input BAMs are similar) without needing to chunk the reads and then merge the alignments back together. If your input BAMs are not similar at all, then you will likely need to chunk to generate uniform alignment jobs.

Below is an example, that shows how to create a groupKey, followed by the transpose, join/combine and groupTuple pattern. There are some caveats and there is some room for optimization:

nextflow.enable.dsl=2

params.original_ref = './ref.fasta'
params.reads_per_job = 10000


process samtools_view_header {

    tag { "${sample_name}:${bam.name}" }

    cpus 1
    memory 1.GB

    input:
    tuple val(sample_name), path(bam)

    output:
    tuple val(sample_name), path("${bam.baseName}.header.txt")

    """
    samtools view \\
        -H \\
        -o "${bam.baseName}.header.txt" \\
        "${bam}"
    """
}

process samtools_name_sort {

    tag { "${sample_name}:${bam.name}" }

    cpus 4
    memory 8.GB

    input:
    tuple val(sample_name), path(bam)

    output:
    tuple val(sample_name), path("${bam.baseName}.nsorted.bam")

    script:
    def avail_mem = task.memory ? task.memory.toGiga().intdiv(task.cpus) : 0
    def mem_per_thread = avail_mem ? "-m ${avail_mem}G" : ''

    """
    samtools sort \\
        -@ "${task.cpus - 1}" \\
        ${mem_per_thread} \\
        -n \\
        -o "${bam.baseName}.nsorted.bam" \\
        -T "${bam.baseName}.nsorted" \\
        "${bam}"
    """
}

process samtools_fastq {

    tag { "${sample_name}:${bam.name}" }

    cpus 1
    memory 2.GB

    input:
    tuple val(sample_name), path(bam)

    output:
    tuple val(sample_name), path("${bam.baseName}.fastq")

    """
    samtools fastq \\
        -O \\
        -T RG,BC \\
        -0 /dev/null \\
        "${bam}" \\
        > \\
        "${bam.baseName}.fastq"
    """
}

process split_fastq {

    tag { "${sample_name}:${fastq.name}" }

    cpus 1
    memory 1.GB

    input:
    tuple val(sample_name), path(fastq)

    output:
    tuple val(sample_name), path("${fastq.baseName}.${/[0-9]/*5}.fastq")

    """
    split \\
        --suffix-length=5 \\
        --additional-suffix=".fastq" \\
        -d \\
        -l "${params.reads_per_job*4}" \\
        "${fastq}" \\
        "${fastq.baseName}."
    """
}

process bwa_index {

    tag { fasta.name }

    cpus 1
    memory 12.GB

    input:
    path fasta

    output:
    tuple val(fasta.name), path("${fasta}.{amb,ann,bwt,pac,sa}")

    """
    bwa index "${fasta}"
    """
}

process bwa_mem {

    tag { "${sample_name}:${fastq.name}" }

    cpus 8
    memory 12.GB

    input:
    tuple val(sample_name), path(fastq), path(header)
    tuple val(idxbase), path(bwa_index)

    output:
    tuple val(sample_name), path("${fastq.baseName}.bam")

    script:
    def task_cpus = task.cpus > 1 ? task.cpus - 1 : task.cpus

    """
    bwa mem \\
        -p \\
        -t ${task_cpus} \\
        -C \\
        -H <(grep "^@RG" "${header}") \\
        "${idxbase}" \\
        "${fastq}" |
    samtools view \\
        -1 \\
        -o "${fastq.baseName}.bam" \\
        -
    """
}

process samtools_coordinate_sort {

    tag { "${sample_name}:${bam.name}" }

    cpus 4
    memory 8.GB

    input:
    tuple val(sample_name), path(bam)

    output:
    tuple val(sample_name), path("${bam.baseName}.csorted.bam")

    script:
    def avail_mem = task.memory ? task.memory.toGiga().intdiv(task.cpus) : 0
    def mem_per_thread = avail_mem ? "-m ${avail_mem}G" : ''

    """
    samtools sort \\
        -@ "${task.cpus - 1}" \\
        ${mem_per_thread} \\
        -o "${bam.baseName}.csorted.bam" \\
        -T "${bam.baseName}.csorted" \\
        "${bam}"
    """
}

process samtools_merge {

    tag { sample_name }

    cpus 1
    memory 4.GB

    input:
    tuple val(sample_name), path(bams)

    output:
    tuple val(sample_name), path("${sample_name}.bam{,.bai}")

    """
    samtools merge \\
        -c \\
        -p \\
        "${sample_name}.bam" \\
        ${bams}
    samtools index \\
        "${sample_name}.bam"
    """
}

workflow {

    original_ref = file( params.original_ref )
    bwa_index( original_ref )

    Channel.fromPath('./data/*.bam')
        | map { tuple( it.baseName, it ) }
        | set { sample_bam_files }

    samtools_view_header( sample_bam_files )

    sample_bam_files
        | samtools_name_sort
        | samtools_fastq
        | split_fastq
        | map { sample_name, fastq_files ->
            tuple( groupKey(sample_name, fastq_files.size()), fastq_files )
        }
        | transpose()
        | combine( samtools_view_header.out, by: 0 )
        | set { realignment_inputs }

    bwa_mem( realignment_inputs, bwa_index.out )
        | samtools_coordinate_sort
        | groupTuple()
        | map { sample_name, bam_files ->
            tuple( sample_name.toString(), bam_files )
        }
        | samtools_merge
        | view()
}

There are some caveats however, in no particular order:

  • It unravels your shell scripts. I think this is a better solution though because it gives you control over each of the components.

  • No CRAM support. But this should be easy to add if needed.

  • A single process to read a BAM header is wasteful, especially if the input BAM is large. More time will be spent localizing the file than actually doing work.

  • Replace samtools sort -n with samtools collate. There's really no need to name sort. You may prefer to pipe the output of samtools collate into samtools fastq for avoid writing an extra file. This might be useful if you need to operate on an extremely large number of large BAMs.

  • The 'samtools_fastq' step should write gzipped output with low compression. This means replacing coreutils split with something else. A solution using Biopython could be useful here to read and write gzipped output in a single pass.

  • The last 'map' command in the workflow is unnecessary. I just thought I'd show that calling .toString() on the GroupKey gives you back an actual String object.

  • Also, it's not necessary to pass around a tuple of 'sample name' and paths as process inputs and outputs. There's no harm in doing this though. When you are more confident with Nextflow Channels, refactoring these processes should be easy.

  • My samtools sort processes ensure that they do not oversubscribe the nodes that they land on. They use approximately 320% - 380% CPU, depending on the size of the BAM. The resource requests for the other processes are not optimized at all.

  • With newer versions of samtools, BAM indexing can be done at the same time as merging, with the --write-index option.

$\endgroup$
8
  • $\begingroup$ Thanks, @Steve. You've made several important points. Question about your suggested workflow code: wouldn't .splitFastq fail since reads.fastqs contains a tuple including the sample_file_input.baseName? I've been battling how to work with the tuples of tuples in various scenarios (e.g., in this situation and after a .groupTuple()). $\endgroup$ Feb 12 at 2:46
  • $\begingroup$ Also, per your suggestion, I've been working to have GENERATE_FASTQS_PROC split the .fastqs itself and emit those. That's causing me a headache. The level of abstraction in Nextflow has really made these details quite elusive for me. $\endgroup$ Feb 12 at 2:47
  • $\begingroup$ I added an update to my question explaining what I feel like I should be able to do and what I ultimately want to get. $\endgroup$ Feb 12 at 3:43
  • $\begingroup$ Thanks, @Steve. I think the interleaved option might be the simplest. Great suggestion. Still, I'm shocked by how difficult it is to get files out of a process without mixing them between samples. What kills me is, I know all files for a single sample are safely contained in a unique directory, but I'm forced to mix them up with files from other samples in the same Channel, so I'm left with convoluted methods to try to distinguish them from each other. If I emit the path to the directories and then use Channel.fromPath(), it would be sooooo easy. But that's not possible, AFAIK. $\endgroup$ Feb 12 at 21:31
  • 1
    $\begingroup$ This was a huge help, @Steve. You've also made several additional helpful suggestions. It had never dawned on me that samtools was oversubscribing because -@ is for additional threads, though that should have been completely obvious. Really appreciate your help. I'm still not a fan of NF output Channels, but maybe I'll get there. :-) $\endgroup$ Feb 16 at 14:57
2
$\begingroup$

Update 3

While the code in Update 2 works, @Steve pointed out that my approach may not be ideal in cloud environments (see comments). He also provided a Nextflow solution that works. I'm still not a fan of how Nextflow output: Channels work, but that may be from my inexperience.

Update 2

Pretty sure I got it this time. The following code works. I'm getting all of the tuples I was expecting.

Here's the updated code for my workflow to create the tuples from each line in the sample_tuple_file I created from a bash for loop. I spent quite a while trying to get a combination of .splitText() and .splitCSV() to work, but for some reason, you cannot chain as follows: it.splitText().splitCSV().

     GENERATE_FASTQS_PROC.out.sample_tuple_file
         .map {
            it.splitText()
         }
         .flatten()
         .map {
             toks = it.trim().split(',')
             sample_name = toks[0]
             fq = toks[1]
             rg = toks[2]
             tuple(sample_name, fq, rg)
             }
         .view()
        .set { sample_fq_tuples_ch }

Update 1

I thought I had it, but I'm actually only getting ONE of the tuples from my GENERATE_FASTQS_PROC.out.sample_tuple_file.map { } statement (the last one). And I have no idea why...


Original posted answer

Despite @Steve's exceptional help and patience with me, I have really battled to figure out how to work with the Nextflow output: Channels, but simply could not figure out any solution using them.

I finally figured out a super simple method that is much easier than dealing with Nextflow Channels, which I share here in case it will help others. I personally will probably use this method from now on, unless someone can convince me that's a bad idea. In my view, using my approach gives me complete control over the way outputs and inputs are handled, and I don't have to worry whether Nextflow is doing what I hope it's doing.

generate_fastqs.sh generate_fastqs.sh takes a single .bam -> sorts by name -> converts to .fastq (now w/ samtools fastq, thanks to @Steve), and splits them with Unix split. Then, at the end of the script, I create a single file with all of the tuples that I want to pass to my next process.

Here are the key lines for these steps (including exit codes to ensure Nextflow receives error notifications). The last step is the key (for me, anyway):

...

#####################
# Sort .bam by name #
#####################
if ! time samtools sort -@ $n_threads -n -m ${mem_per_thread}G $sample_input > $sample_sorted_by_name; then
    echo "ERROR: sort failed"
    exit 1
fi

...

#############################
# Create interleaved .fastq #
#############################
if ! time samtools fastq -o $fq -s /dev/null -0 /dev/null --reference $original_ref $sample_sorted_by_name; then
    echo "ERROR: bamtofastq failed"
    exit 1
fi

...

################
# Split fastqs #
################
if ! time split --additional-suffix=".fastq" -dl $n_lines $fq "${split_dir}/${sample_name}.interleaved_R1_R2.split."; then
    echo "ERROR: split failed"
    exit 1
fi

...

###########################################################
# Store sample name, .fastq, and RG to file for alignment #
###########################################################
RG_file="${sample_name}.tuples.txt"

# Create empty file
> $RG_file
for fq_file in $PWD/${split_dir}/*
do
echo "${sample_name},$fq_file,$RG" >> $RG_file    # <-- Creating future tuples in file, here (comma separated).
done

GENERATE_FASTQS_PROC

In GENERATE_FASTQS_PROC, emit one simple file--the tuple file I created myself with a simple bash loop (so easy and feels so good):

    emit:
        path("*.tuples.txt"), emit: sample_tuple_file

Workflow

In my workflow, I call GENERATE_FASTQS_PROC, access the file from the Channel, read each line, and create a Nextflow tuple for each line (again, the lines I so easily created in a bash for loop).

    /*
     * Generate fastq files from the original .(cr|b)am file and split into
     * mini .fastq files with `n_reads_per_run` in each mini .fastq
     */
    GENERATE_FASTQS_PROC(sample_input_file_ch, n_reads_per_run, original_ref)


    /*
     * Get sample name, fastq, and read group (RG) from the tuple file using Nextflow
     */

     GENERATE_FASTQS_PROC.out.sample_tuple_file
        .map{
            sample_tuple_file ->
                sample_tuple_file.withReader {          // Open file
                    it.eachLine {                       // Loop over each line in file
                        line ->
                            toks = line.split(',')      // Split by comma
                            sample_name = toks[0]       // Get each item
                            fq = toks[1]
                            rg = toks[2]

                            tuple(sample_name, fq, rg)  // Create the simple tuple that has eluded me for two days!
                    }
                }
        }
        .set { sample_fq_tuples_ch }                    // Assign it to a new Channel name and move along!


    /*
     * Align fastq file pairs
     */

    ALIGN_FASTQ_BWA_PROC(sample_fq_tuples_ch, align_to_ref,
                             align_to_ref_tag, original_ref, output_format)
```
$\endgroup$
7
  • $\begingroup$ Ugh. I thought I had it, but I confused myself with my own print statements. I'm actually only getting ONE of the tuples from my GENERATE_FASTQS_PROC.out.sample_tuple_file.map { } statement. And I have no idea why... $\endgroup$ Feb 13 at 1:51
  • $\begingroup$ Pretty sure I got it working with update #2. Hopefully others can learn from my mistakes and final victory. $\endgroup$ Feb 13 at 3:11
  • $\begingroup$ I can see you have put in quite a bit of effort here, but this is not something I would recommend to others. This is a solution that will work on a local filesystem only, where all files exist under a common root. This is what makes it possible to pick them up later on if you have a pointer to them (your tuples.txt contains file pointers). Since the required files are not declared in the output block, Nextflow will not attempt to "unstage" them (upload them to a remote object store, like a bucket) when you deploy your workflow to the cloud. They will just be lost. $\endgroup$
    – Steve
    Feb 13 at 14:52
  • $\begingroup$ Thanks for including your generate_fastqs.sh. I have a better picture now of what you are trying to do. I will try to update my answer with a working solution soon. I just need to grab a couple of test BAMs. I should have a couple somewhere.. $\endgroup$
    – Steve
    Feb 13 at 15:02
  • $\begingroup$ True. I hadn’t considered cloud file systems. I was worried you’d burst my bubble. :-) $\endgroup$ Feb 14 at 2:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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