First of all, thanks so much for making this pipeline. It's incredibly helpful!
Location
subworkflows/differential_expression.nf — map_transcriptome process
Current code:
"""
minimap2 -t ${task.cpus} -ax splice -uf -p 1.0 "${index}" "${fastq_reads}" \
| samtools view -Sb > "output.bam"
samtools sort -@ ${task.cpus} "output.bam" -o "${meta.alias}_reads_aln_sorted.bam"
"""Problem
samtools sort (without -n) produces coordinate-sorted BAM.
Salmon's documentation explicitly states:
Salmon uses a streaming inference method ... One of the fundamental assumptions of such inference methods is that observations (i.e., reads / alignments) are randomly ordered. Alignments should not be sorted by target or position.
When alignments are sorted by target/position, reads from the same transcript arrive in bursts. This causes path-dependent bias in Salmon's online EM — the first transcript to receive a burst of reads can over-claim multi-mapped fragments, producing artifactually integer-like counts instead of the correct fractional (probabilistic) assignments.
We confirmed this experimentally:
| Sort order |
Read count behavior |
Cause |
Name-sorted (samtools sort -n) |
Fractional counts |
Correct probabilistic EM assignment |
| Coordinate-sorted (current default in wf-transcriptomes) |
Artifactual integer counts |
Streaming EM path-dependent bias |
Proposed fix
Either approach works:
Option A — change to name-sorted:
samtools sort -n -@ ${task.cpus} "output.bam" -o "${meta.alias}_reads_aln_sorted.bam"Option B — remove the samtools sort entirely. Minimap2's native output order follows FASTQ read order, which is naturally random with respect to transcript position and already satisfies Salmon's requirement. This also avoids the extra compute and I/O:
minimap2 -t ${task.cpus} -ax splice -uf -p 1.0 "${index}" "${fastq_reads}" \
| samtools view -Sb -o "${meta.alias}_reads_aln_unsorted.bam"(Note: samtools index on a name-sorted or unsorted BAM will fail with NO_COOR reads not in a single block at the end, so the indexing step — if present downstream — should also be removed for Salmon's input.)
Additional context
The reference_assembly.nf subworkflow has the same pattern (samtools sort without -n), though its downstream tools (StringTie, etc.) may legitimately require coordinate-sorted BAM — so the fix only applies to the Salmon input path in differential_expression.nf.
First of all, thanks so much for making this pipeline. It's incredibly helpful!
Location
subworkflows/differential_expression.nf—map_transcriptomeprocessCurrent code:
Problem
samtools sort(without-n) produces coordinate-sorted BAM.Salmon's documentation explicitly states:
When alignments are sorted by target/position, reads from the same transcript arrive in bursts. This causes path-dependent bias in Salmon's online EM — the first transcript to receive a burst of reads can over-claim multi-mapped fragments, producing artifactually integer-like counts instead of the correct fractional (probabilistic) assignments.
We confirmed this experimentally:
samtools sort -n)Proposed fix
Either approach works:
Option A — change to name-sorted:
Option B — remove the
samtools sortentirely. Minimap2's native output order follows FASTQ read order, which is naturally random with respect to transcript position and already satisfies Salmon's requirement. This also avoids the extra compute and I/O:(Note:
samtools indexon a name-sorted or unsorted BAM will fail withNO_COOR reads not in a single block at the end, so the indexing step — if present downstream — should also be removed for Salmon's input.)Additional context
The
reference_assembly.nfsubworkflow has the same pattern (samtools sortwithout-n), though its downstream tools (StringTie, etc.) may legitimately require coordinate-sorted BAM — so the fix only applies to the Salmon input path indifferential_expression.nf.