Skip to content

Commit 0fa5563

Browse files
rob-pclaude
andcommitted
fix(chain): count decoy hits in fillMemCollection (decoy/transcript equal footing)
findOptChain early-returns `if (maxHits == 0)` and builds no chains. maxHits came from fillMemCollection, which returned `maxNonDecoyHits` — the max MEM-hit count over *non-decoy* references only. A mate that maps ONLY to a decoy (e.g. the intronic mate of a genomic / pre-mRNA fragment, which maps uniquely to the genome and to no transcript) therefore returned 0 and was given no chains at all. It could never form the concordant decoy pair its (exonic) mate participates in, so the genomic fragment leaked through as a spurious single-mate transcript orphan instead of being attributed to the decoy. Count hits over ALL references (transcript and decoy) so the two are chained on equal footing — matching the salmon Rust implementation, whose chaining/consensus considers decoy and transcript candidates together. At the default FILTER_AFTER_CHAINING policy the only behavioral change is that decoy-only mates are now chained (and can pair); transcript-only reads are unaffected. Validated on SRR1039508 + GRCh38 decoy index: a fragment with one exonic mate (maps to transcript+genome) and one intronic mate (genome only) now forms the concordant genome pair and is correctly counted as decoy. On the 3M-read subset C++ goes from 94.019% mapped / 42,044 decoy to 93.742% / 154,800 — matching the Rust port's 93.749% / 154,284 (was a 0.27-pt over-report of transcript mappings). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01B7JMur5DmDpECddErpi2JS
1 parent 41a9245 commit 0fa5563

1 file changed

Lines changed: 14 additions & 3 deletions

File tree

src/MemChainer.cpp

Lines changed: 14 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -71,7 +71,17 @@ size_t MemClusterer::fillMemCollection(std::vector<std::pair<int, pufferfish::ut
7171
return 0;
7272
}
7373

74-
size_t maxNonDecoyHits{0};
74+
// Max MEM-hit count over ALL references (transcript AND decoy), so transcript
75+
// and decoy targets are chained on equal footing. Previously this counted only
76+
// non-decoy references, so a mate mapping *only* to a decoy (e.g. an intronic
77+
// read of a genomic / pre-mRNA fragment) returned 0 here, and findOptChain's
78+
// `if (maxHits == 0) return false` then built it no chains at all — so it could
79+
// never form the concordant decoy pair its mate participates in, and the
80+
// genomic fragment leaked through as a spurious transcript orphan. Counting
81+
// decoy hits equally lets such a mate be chained and paired (matching the Rust
82+
// implementation, whose chaining/consensus considers decoy and transcript
83+
// candidates together).
84+
size_t maxHits{0};
7585
size_t totSize{0};
7686
for (auto &hit : core::range<decltype(hits.begin())>(hits.begin(), hits.end())) {
7787
auto &refs = hit.second.refRange;
@@ -104,12 +114,13 @@ size_t MemClusterer::fillMemCollection(std::vector<std::pair<int, pufferfish::ut
104114
auto& refHits = trMemMap[encodeRefKey(tid, refPosOri.isFW)];
105115
refHits.emplace_back(memItr, refPosOri.pos, refPosOri.isFW);
106116
auto nh = refHits.size();
107-
maxNonDecoyHits = (tid < firstDecoyIndex) ? std::max(nh, maxNonDecoyHits) : maxNonDecoyHits;
117+
maxHits = std::max(nh, maxHits);
108118
//}
109119
}
110120
}
111121
}
112-
return maxNonDecoyHits;
122+
(void)firstDecoyIndex; // no longer used to gate hit counting (see maxHits above)
123+
return maxHits;
113124
}
114125

115126
bool MemClusterer::findOptChain(std::vector<std::pair<int, pufferfish::util::ProjectedHits>> &hits,

0 commit comments

Comments
 (0)