Skip to content

eval: placement comparison is strand-asymmetric under soft-clipping #23

Description

@nh13

holodeck eval compares the leftmost mapped position to the leftmost true position (eval.rs:136, read.rs:157). Strand is encoded in the read name and parsed, but never used in the comparison. Since the leftmost coordinate is the read's 5' end on the forward strand and its 3' end on the reverse strand, soft-clipping is scored differently depending on which strand a read came from:

clip forward reverse
5' (cycle 1) penalized ignored
3' (last cycle) ignored penalized

Error rates rise toward the 3' end, so aligners clip 3' ends more often, and the default --wiggle 5 turns any clip longer than five bases into a mismap. The reported mismapping rate is consequently a function of strand and of the error model as much as of the aligner.

Proposal

Add --placement to choose which coordinates are compared, and report clipping fidelity as columns rather than folding it into the pass/fail decision.

  • cycle1 (new default) compares the 5' end of the read on both strands. It is strand-symmetric, anchored on where the molecule began being read, and unaffected by 3' trimming. A strand mismatch counts as mismapped.
  • leftmost is today's behaviour. Retained because it matches wgsim_eval and most published accuracy tables, but documented as strand-asymmetric.
  • overlap is correct when the true and mapped spans reciprocally overlap by at least --min-overlap (default 0.1, matching paftools.js mapeval). It answers "the right locus at all?" and ignores clipping at both ends.

Two columns are added, computed only over reads that pass the placement test, since deltas from mismapped and unmapped reads are meaningless:

  • pct_both_ends, the percentage with both |Δ5'| and |Δ3'| within --wiggle.
  • mean_end_delta, the mean of |Δ5'| + |Δ3'| in bases.

Both are defined over both ends under every mode so the numbers stay comparable across --placement values. Under cycle1 the 5' term is near zero by construction, so they reduce to a clean 3' clipping signal without hard-coding which end is which.

Truth spans

cycle1 and overlap need the true reference end, which the encoded read name does not carry. It can be approximated as min(read length, fragment length), exact unless the read spans a haplotype indel. The approximation affects only reverse-strand reads under cycle1, as a forward read's 5' anchor is the encoded position itself, so an indel longer than --wiggle inside a read yields spurious reverse-strand mismaps.

--truth is presently a stub (eval.rs:61). Implementing it to take exact spans from the golden BAM removes the approximation; until then eval should warn when cycle1 runs without it.

Considered and not proposed

  • both-ends (Σ|Δ|) as a mode. It scores a read placed 10 kb away and a read placed correctly but clipped 20 bases identically, which is the conflation the two columns exist to remove.
  • midpoint and contig-only. Neither answers a question that cycle1 and overlap do not already bracket.
  • Oracle scoring, in which an alignment is correct if it is among all optimally-scoring positions. This is the right answer for repeats where the true position is genuinely ambiguous and the aligner is not wrong to pick an equally good copy, but it requires alignment inside holodeck. Worth a note in the docs regardless, as it explains residual mismaps at high MAPQ.
  • Pair-level scoring. Aligners use mate rescue and per-read counting double-counts a pair that fails together, but that changes the unit of aggregation rather than the placement anchor, so it belongs in its own issue.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions