feat(quant): --spliced, mapping output in genome coordinates (stacked on #1141, #1142) - #1123
Open
BenjaminDEMAILLE wants to merge 12 commits into
Open
BenjaminDEMAILLE wants to merge 12 commits into
BenjaminDEMAILLE wants to merge 12 commits into
Conversation
BenjaminDEMAILLE
marked this pull request as draft
August 12, 2026 13:32
BenjaminDEMAILLE
marked this pull request as ready for review
August 12, 2026 19:10
rob-p
added
the
postponed
Collaborator
|
Marking this postponed rather than closed, and converting to draft. The features are genuinely nice — qualities, real CIGARs, NM/MD, read groups, unmapped records, spliced genome-coordinate output would make What would change the calculus is demonstrated need: users showing up with concrete workflows that require salmon's mapping output to be a full SAM/BAM (rather than using a dedicated aligner for that job, which remains the natural tool). If those requests accumulate, this PR is the right starting point and the work here won't be wasted — that's why it's parked (labeled |
rob-p
marked this pull request as draft
August 21, 2026 03:14
BenjaminDEMAILLE
force-pushed
the
feat/bam-enrichment
branch
from
August 21, 2026 20:23
b2fd2d3 to
1531ebd
Compare
BenjaminDEMAILLE
changed the title
Contributor
Author
|
Thanks, and the scope objection is fair, so I acted on it rather than arguing it: this PR is split. #1141 is the half that is not new territory.
None of that is salmon promising to be an aligner. It is the diagnostic view misdescribing salmon's own mappings, which is worse than a minimal view. #1141 fixes exactly that and adds no new output mode. Its cost is zero when no mapping output is requested (a run that writes no records never enters the realization module), and about 2x on This PR keeps the part that genuinely is new surface: the genome projection, On demonstrated need, I have added the concrete workflow to the description rather than claiming a queue of users: Picard |
BenjaminDEMAILLE
force-pushed
the
feat/bam-enrichment
branch
from
August 21, 2026 20:52
1531ebd to
305f80d
Compare
BenjaminDEMAILLE
changed the title
Contributor
Author
|
Follow-up on the split: it is now three PRs rather than two, because two of the things I had bundled as "correctness" are not the same kind of claim.
The distinction matters for your scope argument: a diagnostic view that omits fields is defensible, and I am not arguing #1142 is obligatory. A diagnostic view that describes an alignment salmon did not make is not, whatever its scope, which is the whole of #1141. |
Closed
5 tasks
Mapping runs ksw2 in KSW_EZ_SCORE_ONLY: it needs the score to rank and filter placements, not the path that produced it, and skipping the traceback matrix is the right trade for a kernel that runs once per candidate across the library. That is also why no CIGAR survives to the record writer, and why the writer has been synthesizing one from the read length. This module re-aligns a placement with traceback, on demand, producing the CIGAR plus the edit distance and MD string that describe it. Nothing calls it yet. Three things keep it affordable, since it will run once per written record: * Provably-ungapped fast path. Any alignment containing a gap pays at least gap_open + gap_extend and can do no better than match everywhere else, so it cannot exceed perfect - (gap_open + gap_extend). The ungapped alignment scores perfect - mismatches * (match + mismatch). When the second is at least the first, no gapped alignment can win and no DP runs at all. At salmon's defaults that admits up to one mismatch, which covers most reads at realistic error rates. * Traceback skipped when the banded optimum equals the ungapped score over the full read: filling and walking the traceback matrix is the expensive half. * Anchor repair by comparison rather than by DP. The mapper's position comes from a seed chain, which a mismatch near the 5' end pushes rightwards; sliding the read and counting mismatches finds the true start directly, unbounded by any band. An exact probe locates candidates first, with the exhaustive sweep as the fallback for the fraction of a percent the probe cannot place. The `ambiguous` slice names reference offsets whose base indexing had to replace, so NM and MD can report the base the user supplied rather than the substitute. Callers pass an empty slice until the index records them. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Indexing replaces non-ACGT bases with ACGT, as salmon has always done, because the k-mer structure cannot hold them. Nothing recorded where, so an alignment over such a position counted a match where the input had an N, and MD named the substitute instead of the base the user supplied. The build now writes ambiguous.bin, the offsets it replaced: 40 bytes for the whole human transcriptome. Alignment still runs against the cleaned sequence, the only thing it can do, but a record can now report what was actually there. An index built before this file existed loads as "none recorded" and behaves exactly as it did, so nothing has to be rebuilt to keep working; rebuilding is what earns the correction. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
--writeMappings has always been a diagnostic view of salmon's mappings rather than a general-purpose aligner output, and that is the right scope. But the view stated things that were not true of the mappings it described. * The CIGAR was synthesized from the read length: <readLen>M plus a clip at a transcript end. A placement with an indel was written as if it had none, so every base past the gap sat at the wrong reference offset, and the record contradicted its own placement. * AS was the mapper's chain score, and nothing else in the record could be checked against it. Auditing AS against a recomputed NM over a whole output found 2.4% of records self-contradicting: a 15M1I1D15M1I1D... sawtooth at NM:i:80 on a 100 base read, beside an AS claiming about three mismatches. * MAPQ was the constant 1, so `samtools view -q 255`, the standard "uniquely placed" filter, discarded the entire file. Records now carry the CIGAR of an alignment actually realized for that placement, with real I/D operations, the NM and MD that describe it, the score of that same alignment, and MAPQ derived from the placement count on STAR's scale (255 unique, 3, 1, 0), which is what the tools downstream already filter on. Three separate causes behind the self-contradicting records, each of which alone produces a false one: * The anchor search reached only indel_margin. That bounds indel size and says nothing about how wrong a seed-derived position can be: a chain on a repeat or a paralogous stretch can imply a start most of a read away. One offending read belonged 69 bases earlier, where it matched with three mismatches. Out of reach, the banded DP still returned an alignment: a sawtooth of one-base indels walking the read diagonally across the reference. The search now covers a read length either side. * The search gave up entirely when the read did not fit inside the reference at its anchor, which is exactly the case most worth looking at: an anchor near a transcript end becomes a soft-clipped record that throws bases away when the read often belongs further in and matches in full. "Does not fit" is now treated as "matches nothing" so the search runs. * ksw2's extension pins the read's first base to the window's first base, so a realized alignment can grow rightwards but never shift left, and a mismatch near the 5' end left the DP no way to say "the read started earlier" except by opening a leading insertion: 26I34M at NM 26 where the answer is 60M at NM 1. A leading insertion of n bases is itself the statement that the read began n bases earlier, so it is also the correction. A repaired anchor is a hypothesis, not a verdict: a read with a genuine deletion also looks badly placed to an ungapped scan, so both anchors are realized and the better score wins. An alignment scoring below the mapper's own min_accepted_score for the bases it claims is refused, and the record falls back to the synthesized CIGAR with no NM or MD, which at least says "no base-level alignment here". Cost, measured on 1M fragments at -p 8, best of 3: a run without mapping output is unchanged (0.89s), because nothing is realized unless records are written; --writeBam goes from 1.04s to 1.94s on a single-mapping fixture and 1.36s to 3.28s at 2.5 placements per fragment, the cost being per record. NM and MD match samtools calmd across 2M records, including for the transcript containing an N. QUAL, M5/UR, MC/MQ, ZW, read groups and unmapped records are all still absent. Those are omissions the format allows, unlike a CIGAR that does not describe the alignment beside it, and they are the subject of a separate change. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Everything here is an omission rather than a falsehood: SAM permits `QUAL` to be `*`, and `M5`, `UR`, `MC`, `MQ` and `ZW` to be missing. The file was still less useful than it had to be, and salmon already had every value in hand. * QUAL, taken from the FASTQ and reversed alongside SEQ on the reverse strand, so `samtools fastq` can round-trip the library back out of the BAM. It is not behind a flag: whenever the input carries qualities the record does. A quality string whose length disagrees with the read is dropped rather than written, since that mismatch is what makes a record unreadable, and it warns once. * `@SQ M5`, the MD5 of each reference, and `UR`, the index it came from. M5 is what ties a BAM to the exact bases it was made against, and what CRAM needs. Where indexing had to replace a base the checksum hashes the `N` that was there, so the header names the reference the user supplied. The whole transcriptome is hashed in parallel at startup: 0.16s to 0.08s on a 44 MB transcriptome, and it matches `samtools dict`. * `MC` and `MQ`, the mate's CIGAR and mapping quality, which only became meaningful once the CIGAR was real and MAPQ was not a constant. A mate-aware tool no longer has to sort the file to find them. * `ZW`, the placement's equivalence-class weight, which salmon already computed and nothing else in the file reports. * `@RG` and `RG:Z` via `--rgLine`, in the spelling bwa and STAR users already type, parsed before any mapping work so a malformed line fails immediately rather than after the run. Read groups are how a merged BAM keeps track of which sample each record came from. * `@HD VN:1.6 SO:unsorted GO:query` and a `@CO` stating the one ordering guarantee. `SO:unknown` was not false, but it hid a property the output does have: records for one fragment are contiguous, so a reader can pair mates in a single streaming pass instead of sorting first. `--writeQualities` becomes accepted-and-default rather than warned-as-missing, since qualities are now written whenever the input has them. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
A BAM containing only the reads that mapped is a record of one filtering decision, not of the experiment, and the reads it drops are exactly the ones anybody debugging a low mapping rate wants. The format has a place for them: FLAG 0x4, no reference, no position, sequence and qualities intact. --sampleUnaligned no longer requires --sampleOut (which remains unimplemented). On its own it adds those records to --writeMappings/--writeBam, and warns when there is no mapping output for them to go into. It also covers an orphan's unmapped mate, which is a read that did not align just as much as one whose whole fragment failed. The mapped mate carries MATE_UNMAPPED, so without the mate's record the file described a read it did not contain: samtools fastq could not round-trip it and anything following the flag chased nothing. The mate is written at the mapped mate's position, so a sort keeps the pair together, and the two records point at each other. Without --sampleUnaligned the mapped record keeps RNEXT '*', since naming a mate that was not written is worse than admitting there is none. With this, samtools fastq round-trips a whole real library out of the BAM: all 457366 pairs of DRR028171 returned, not one sequence or quality string altered. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…rdinates A transcript has no introns, so a read across an exon junction is contiguous in transcript coordinates and its CIGAR never contains an N. That is an accurate description of the alignment salmon made, but it is not one a genome browser or a junction counter can read. Add salmon_quant::splice, which cuts a transcript-coordinate alignment at the exon boundaries an annotation declares and fills the gaps with N operations, producing the spliced genome alignment a genome aligner would have made. Reverse-strand transcripts run against the genomic axis, so their alignments are turned around: the CIGAR is reversed, the reported position becomes the alignment's leftmost genomic base, and the projection reports the flip so the caller can apply it to SEQ, QUAL and the strand bits together. Chromosome lengths for @sq come from the genome FASTA index when one exists, since an annotation never states them and an understated LN makes readers reject the file. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Wire the genome projector into the record path. With --spliced, @sq names chromosomes instead of transcripts, positions are genomic, and CIGARs carry N operations across exon junctions, so the output is the spliced genome BAM the rest of the RNA-seq toolchain expects rather than a transcriptome one. A projected reverse-strand transcript turns the alignment around, so both mates' strand bits swap, their stored bases are reverse-complemented once more, and TLEN is recomputed over the genomic span, which introns make longer than the transcript span. NM and MD survive projection unchanged: they describe the read against the transcript's bases, which are the genome's bases inside an exon. A placement on a reference the annotation does not describe (a decoy, or a transcript absent from the GTF) is dropped rather than half-projected, since the surviving record would claim a mate that is not in the file. --spliced requires --annotation and --genome, checked before mapping starts rather than when the first record is written. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
bramble's GTF reader is commented "1-based inclusive -> 0-based half-open" but keeps the 1-based start and adds one to the end, so exon lengths are right while every coordinate sits one base too high. That shift is invisible in the intron lengths, since both exon ends move together, and only shows up in the reported position: an end-to-end run put a read at chr1:272 that belonged at chr1:271. Subtract the one where the exon blocks are built. Verified against the genome with samtools calmd: the NM and MD it recomputes from the projected positions and CIGARs now match the ones salmon wrote, which they cannot do if the positions are off by a base. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…ment A reverse-strand transcript runs against the genomic axis, so projecting its alignments turns them around: the CIGAR reverses and the stored sequence is reverse-complemented. MD was left in transcript orientation, so it described bases that were no longer the ones under the CIGAR beside it. We wrote MD:Z:71T28 where the genome says 28A71. The existing spliced fixture could not show this. It has a reverse-strand transcript, but its reads carry no mismatches, and an MD of a single run length is unchanged by reversing it. Only a mismatch makes the orientation visible. Building a fixture that could show it, 500 genes on a 5 Mb chromosome, three to eight exons each, both strands, 200k fragments at a 1% error rate, found 105686 of 399524 records disagreeing with samtools calmd recomputed against the genome. Every one is now identical, including the 161836 that span a junction. NM needs no such treatment: it counts edits, and reversing an alignment does not change how many there are. Also bound the probe search. A probe is useful because it is rare, so one that matches dozens of places within a few hundred bases is sitting in a tandem repeat and says nothing about which of them the read came from. Collecting them all turned a linear scan quadratic on exactly the low-complexity sequence that produces the most hits, which real transcripts are full of. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…ations The projector cut alignments at exon boundaries only within an operation. A junction landing exactly between two of them was missed, so a deletion ending on an exon edge left the following match crossing the intron with no N at all: the bases after the boundary silently claimed reference from the wrong side of it, and the record read 98M1D2M where the truth is 98M1D...N...2M. Real data found it. Human chr22 with GENCODE v44, 300k fragments, 1.85M records: 431 of them disagreed with samtools calmd recomputed against the chromosome. No synthetic fixture could show it, because it needs an indel to land exactly on an exon edge, which random exon boundaries almost never do. Six records still differ, in form rather than in bases: a deletion cut in two by an intron is spelled ^GG where calmd writes the canonical ^G0^G. Same NM, same reconstructed reference. Documented in the module rather than papered over. Also confirms the transcript-coordinate path is exact: on the same library without --spliced, every record carrying NM and MD matches calmd, and the only records that differ are the ones salmon deliberately writes without them. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
MD was computed while the alignment was still in transcript coordinates and then carried across projection, adjusted only for strand. That holds for every operation except a deletion the projection has to cut in two: two deletions are spelled with an empty match between them, ^G0^G, not ^GG. The merged form disagrees with the CIGAR beside it, and samtools calmd writes the canonical one. Re-derive the string from the projected operations instead. The reference bases are unchanged by projection, only their grouping is, so expanding the transcript's MD to one entry per reference base and re-emitting it against the projected operations gets the grouping right by construction: the separating run length falls out of the normal rule that one sits between every pair of tokens. The strand flip then applies to the re-derived string. This closes the last discrepancy on real data. Human chr22 with GENCODE v44, 300k fragments, 1.85M records of which 685k span a junction: every record carrying NM and MD now matches samtools calmd recomputed against the chromosome, where six did not before. The 500-gene synthetic fixture is likewise exact. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Two pieces of the genome-coordinate output that the correctness split (COMBINE-lab#1141) deliberately left out, since neither means anything without a projector. XS:A reports the transcript's genomic strand. StringTie and Cufflinks read an intron the wrong way round without it, and silently mis-assemble. A projection can refuse every placement of a fragment when the annotation does not describe its transcript, and the fragment then produced no records at all: on a fixture where half the transcripts were absent from the GTF, half the fragments disappeared with no warning and no accounting. Such a fragment is now written as unaligned when --sampleUnaligned is on, and the projector counts what it refuses so the run ends with a warning naming the number instead of a silently short file. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
BenjaminDEMAILLE
force-pushed
the
feat/bam-enrichment
branch
from
August 30, 2026 11:44
305f80d to
20b0e01
Compare
BenjaminDEMAILLE
marked this pull request as ready for review
September 26, 2026 11:16
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Rebased and narrowed, now into a three-PR stack. Following the scope objection below, this PR was split:
ASdescribing a different alignment from the CIGAR beside it,MAPQpinned to 1. 2.4% of records were self-contradicting.QUAL,@SQ M5/UR,MC/MQ,ZW,@RG/--rgLine, unmapped records. Omissions, not falsehoods, so a separate decision.splice.rs,--spliced,XS, the projection accounting. Seven commits, ~1.4k lines.GitHub cannot point a base branch at a fork, so until #1141 and #1142 merge the diff shown here still includes them.
What
--spliceddoesA transcript has no introns, so a read spanning an exon junction is contiguous in transcriptome coordinates and its CIGAR never contains an
N.--splicedre-expresses each record in genome coordinates: the alignment is cut at the exon boundaries the annotation declares, the gaps reappear asN,@SQnames chromosomes, and records gainXS:Afor the transcript's genomic strand.It needs
--annotation <gtf|gff>for the exon structures and--genome <fasta>for the chromosome lengths@SQmust declare (a.faiis used when present). A transcript on the minus strand runs against the genomic axis, so its alignments are turned around whole: CIGAR, sequence, qualities, strand bits,TLENandMDtogether.M5is omitted, because the references are then chromosomes salmon never loaded and a guessed checksum is worse than none.Placements on a reference the annotation does not describe (a decoy, or a transcript absent from the GTF) are refused rather than half-projected, counted, and reported as a warning at the end of the run. With
--sampleUnaligneda fragment whose every placement was refused is written as unaligned instead of vanishing.Why the genome coordinates are the point
The correctness half of this work stands on its own, but it does not remove a second alignment: a transcriptome-coordinate BAM cannot be fed to the tools that need genome coordinates. Two concrete workflows in our pipeline:
RNA-seq QC. Picard
CollectRnaSeqMetrics, RSeQC (read_distribution.py,geneBody_coverage.py,inner_distance.py) and Qualimaprnaseqall take a genome-coordinate BAM plus an annotation. A salmon-quantified sample therefore gets aligned twice today: once by salmon for the quantification we actually want, and once by STAR purely to produce a file for QC. The second alignment is 100% of the QC cost and none of the analysis. With--splicedthe QC reads the same alignments the quantification used, which is also the more honest QC: it describes the file the numbers came from rather than a different aligner's view of the same library.Allele-specific expression. GATK
ASEReadCounterand phASER need genome coordinates and real CIGARs withNM/MD, the same second alignment for the same reason.I am one lab with one pipeline, not a queue of requests, so this is evidence rather than demand. If it is not enough, a decision criterion would help more than a merge: what would count? Two or three independent reports? A user issue naming a tool by name? Happy to keep this parked until then, and #1141 stands on its own either way.
Verification (unchanged from the original review)
Against samtools 1.24, on 500 genes over a 5 Mb simulated chromosome (three to eight exons, both strands) and on human chr22 (hg38) + GENCODE v44, 300k fragments, 1.85M records of which 685k span a junction:
NM/MDmatchsamtools calmdrecomputed against the genome, on every projected record.samtools quickcheckpasses, SAM round-trips throughsamtools view -b.SEQ/QUALorientation, strand bits,TLENsign andMDall agree after the flip.splice.rs.