diff --git a/src/main/scala/com/fulcrumgenomics/umi/ConsensusCallingIterator.scala b/src/main/scala/com/fulcrumgenomics/umi/ConsensusCallingIterator.scala index b5ab48065..5dba4cde7 100644 --- a/src/main/scala/com/fulcrumgenomics/umi/ConsensusCallingIterator.scala +++ b/src/main/scala/com/fulcrumgenomics/umi/ConsensusCallingIterator.scala @@ -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 { diff --git a/src/main/scala/com/fulcrumgenomics/umi/UmiConsensusCaller.scala b/src/main/scala/com/fulcrumgenomics/umi/UmiConsensusCaller.scala index 9d61fae27..8dc460866 100644 --- a/src/main/scala/com/fulcrumgenomics/umi/UmiConsensusCaller.scala +++ b/src/main/scala/com/fulcrumgenomics/umi/UmiConsensusCaller.scala @@ -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. */ @@ -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 } diff --git a/src/test/scala/com/fulcrumgenomics/umi/CallMolecularConsensusReadsTest.scala b/src/test/scala/com/fulcrumgenomics/umi/CallMolecularConsensusReadsTest.scala index 76620ebac..bac08c6af 100644 --- a/src/test/scala/com/fulcrumgenomics/umi/CallMolecularConsensusReadsTest.scala +++ b/src/test/scala/com/fulcrumgenomics/umi/CallMolecularConsensusReadsTest.scala @@ -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._ /** @@ -131,6 +133,7 @@ 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(), @@ -138,6 +141,7 @@ class CallMolecularConsensusReadsTest extends UnitSpec { minReads = 1, maxReads = Some(maxReads), rejects = Some(rejects), + stats = Some(stats), readGroupId = "ABC" ).execute() @@ -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" } } }