hold --mgSplit second-pass MAPQ to the whole-genome value - #1988
Closed
glennhickey wants to merge 1 commit into
Closed
hold --mgSplit second-pass MAPQ to the whole-genome value#1988glennhickey wants to merge 1 commit into
glennhickey wants to merge 1 commit into
Conversation
Under --mgSplit the second pass maps each chromosome's contigs against a graph
holding only that chromosome, so a repeat that is ambiguous genome-wide has
nothing left to be ambiguous with and comes back uniquely mapped. MAPQ is
inflated exactly where it is least trustworthy.
Adding sequence to a graph only lowers MAPQ and removing it only raises it, so
each pass is an upper bound on the MAPQ the sequence would earn against a
whole-genome pangenome: the first has fewer competitors in the graph and more in
the genome, the second the reverse. The minimum of the two is the tightest bound
available, not merely the cautious reading.
make_mapq_track projects the whole-genome PAF down to (query, interval, mapq) and
apply_mapq_floor rewrites each second-pass MAPQ as min(second pass, first pass),
keeping the original in om:i:.
Three things are deliberate:
raw PAF the track is built before filter_paf, which deletes everything
under minMAPQ -- the strongest evidence the floor has
no bin test the bound is what the first pass said about the query interval,
wherever it put the sequence. A repeat with one copy here and one
on another chromosome is rightly MAPQ 0 whichever of the two the
first pass happened to print, so conditioning on the reference
contig would miss half the cases
no clamp a whole-genome MAPQ of 0 writes 0, and filter_paf deletes the line
when cactus-align reads it. That is not a new filter: it is the
one a whole-genome run already applies, and that --mgSplit loses by
hiding the competing chromosome from the mapper. minMAPQ is read
only to report how many lines that comes to
The reference is exempt. Its alignment reaches the PAF from minigraph and from
rgfa2paf off the rGFA tags, and only the first carries a whole-genome MAPQ, so
the track holds the minigraph copy (0 for a self-mapping) while the bin holds the
rgfa2paf copy (60); crossing them would zero the reference and filter it out.
The floor runs before the collapse PAF is merged in, whose minimap2
self-alignment MAPQs are a different quantity. Without --mgSplit no track is
built and nothing changes, and the step-by-step entry points have no first-pass
track so they skip it too.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0151awXbza1QLwmsYBMW2yqk
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
1 participant
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.
Under --mgSplit the second pass maps each chromosome's contigs against a graph holding only that chromosome, so a repeat that is ambiguous genome-wide has nothing left to be ambiguous with and comes back uniquely mapped. MAPQ is inflated exactly where it is least trustworthy.
This PR remembers the 1st-pass (genome-wide) mapq, and applies it (via
min()) to 2nd-pass mapqs.