Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -53,13 +53,12 @@ class ConsensusCallingIterator[ConsensusRead <: SimpleRead](sourceIterator: Iter
private var collectedStats: Boolean = false

protected val iter: Iterator[SamRecord] = {
// Only mapped reads may contribute to a consensus. An unmapped read has no cigar, and an empty cigar is a prefix
// of every cigar, so `filterToMostCommonAlignment` would match it against whichever alignment group it was tested
// against first rather than rejecting it as a minority alignment. Note that the mapped end of a half-mapped pair
// is still used; only the unmapped end is dropped.
// Templates with no mapped read at all are dropped here; the unmapped end of a half-mapped pair is passed through
// to the caller, which rejects it as `RejectionReason.Unmapped` so that it is counted and written to the rejects
// output rather than silently discarded.
val filteredIterator = sourceIterator
.filterNot (r => r.secondary || r.supplementary)
.filter(r => r.mapped)
.filter(r => r.mapped || (r.paired && r.mateMapped))

// Wrap our input iterator in a progress logging iterator if we have a progress logger
val progressIterator = progress match {
Expand Down
13 changes: 11 additions & 2 deletions src/main/scala/com/fulcrumgenomics/umi/UmiConsensusCaller.scala
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,7 @@ object UmiConsensusCaller {
case object HighDuplexDisagreement extends _RejectionReason("high_duplex_disagreement", "Too many errors between top/bottoms strands", false, false, true)
case object ClipOverlapFailed extends _RejectionReason("clip_overlap_failed", "See https://github.com/fulcrumgenomics/fgbio/issues/1090", false, false, true)
case object NotPrimaryFrPair extends _RejectionReason("not_primary_fr_pair", "Template did not have a single primary FR pair of reads", false, false, true)
case object Unmapped extends _RejectionReason("unmapped", "Read is unmapped", true, true, true)
}

/** Metric class for outputting consensus calling statistics. */
Expand Down Expand Up @@ -351,10 +352,18 @@ trait UmiConsensusCaller[ConsensusRead <: SimpleRead] {
*/
final def consensusReadsFromSamRecords(recs: Seq[SamRecord]): Seq[SamRecord] = {
this._totalReads += recs.size

// Unmapped reads cannot contribute to a consensus: they have no cigar, and an empty cigar is a prefix of every
// cigar, so `filterToMostCommonAlignment` would match one against whichever alignment group it was tested against
// first rather than rejecting it. Rejecting here rather than filtering upstream keeps them visible in the rejects
// output and in the rejection counts. The mapped end of a half-mapped pair is unaffected.
val (mapped, unmapped) = recs.partition(_.mapped)
if (unmapped.nonEmpty) rejectRecords(unmapped, RejectionReason.Unmapped)

// Ensure that mate cigar is set on all read pairs. This is needed when using
// `clipper.numBasesExtendingPastMate` subsequently
updateMateCigars(recs)
val result = consensusSamRecordsFromSamRecords(recs)
updateMateCigars(mapped)
val result = consensusSamRecordsFromSamRecords(mapped)
this._consensusReadsConstructed += result.size
result
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,8 @@ package com.fulcrumgenomics.umi

import com.fulcrumgenomics.bam.api.SamOrder
import com.fulcrumgenomics.testing.{SamBuilder, UnitSpec}
import com.fulcrumgenomics.umi.UmiConsensusCaller.{ConsensusKvMetric, RejectionReason}
import com.fulcrumgenomics.util.Metric
import com.fulcrumgenomics.umi.VanillaUmiConsensusCallerOptions._

/**
Expand Down Expand Up @@ -131,13 +133,15 @@ class CallMolecularConsensusReadsTest extends UnitSpec {
val rlen = 100
val output = newBam
val rejects = newBam
val stats = makeTempFile("call_molecular_consensus_reads_test.", ".txt")

new CallMolecularConsensusReads(
input = halfMappedFamily(rlen).toTempFile(),
output = output,
minReads = 1,
maxReads = Some(maxReads),
rejects = Some(rejects),
stats = Some(stats),
readGroupId = "ABC"
).execute()

Expand All @@ -156,8 +160,14 @@ class CallMolecularConsensusReadsTest extends UnitSpec {
r2[Int](ConsensusTags.PerRead.RawReadCount) shouldBe 2
r2.basesString shouldBe "C" * rlen

// The unmapped R2 must not have been used, whether it was rejected by the caller or filtered before reaching it.
readBamRecs(rejects).foreach { rec => rec.name shouldBe "halfmapped:3" }
// The unmapped R2 must not have been used, and must be accounted for rather than vanishing:
// it is written to the rejects BAM and counted under its own rejection reason.
val rejected = readBamRecs(rejects)
rejected.map(rec => (rec.name, rec.secondOfPair)) should contain theSameElementsAs Seq(("halfmapped:3", true))

val metrics = Metric.read[ConsensusKvMetric](stats).map(m => m.key -> m.value.toString).toMap
metrics(s"raw_reads_rejected_for_${RejectionReason.Unmapped.code}") shouldBe "1"
metrics("raw_reads_considered") shouldBe "6"
}
}
}
Loading