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.
- 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.
- 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.
- 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.
