Skip to content

Commit

Permalink
fixups
Browse files Browse the repository at this point in the history
  • Loading branch information
nh13 committed Mar 21, 2022
1 parent 95e9ec9 commit 5415287
Show file tree
Hide file tree
Showing 2 changed files with 36 additions and 6 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,7 @@ class CorrectOverlappingBases
threadDatum.templateMetric.total += 1
threadDatum.basesMetric.total += template.primaryReads.map(_.length).sum
// corrects
val stats = threadDatum.caller.correct(template)
val stats = threadDatum.caller.call(template)
val overlappingBases = stats.r1OverlappingBases + stats.r2OverlappingBases
val correctedBases = stats.r1CorrectedBases + stats.r2CorrectedBases
if (overlappingBases > 0) {
Expand Down Expand Up @@ -148,6 +148,14 @@ class CorrectOverlappingBases

}

/** Collects the the number of reads or bases that were examined, had overlap, and were corrected as part of
* the [[CorrectOverlappingBases]] tool.
*
* @param tpe template if the counts are per template, bases if counts are in units of bases.
* @param total the total number of templates (bases) examined
* @param overlapping the total number of templates (bases) that were overlapping
* @param corrected the total number of templates (bases) that were corrected.
*/
case class CorrectOverlappingBasesMetric
(
tpe: CountType,
Expand All @@ -166,6 +174,7 @@ case class CorrectOverlappingBasesMetric

sealed trait CountType extends EnumEntry

/** Enumeration for the type of counts in [[CorrectOverlappingBasesMetric]]. */
object CountType extends FgBioEnum[CountType] {
case object Templates extends CountType
case object Bases extends CountType
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -56,12 +56,26 @@ class OverlappingBasesConsensusCaller(onlyMaskDisagreements: Boolean = false,
private val r1Builders = Seq(r1BasesBuilder, r1QualsBuilder)
private val r2Builders = Seq(r2BasesBuilder, r2QualsBuilder)

def correct(template: Template): CorrectionStats = (template.r1, template.r2) match {
case (Some(r1), Some(r2)) if r1.mapped && r2.mapped && r1.matesOverlap.contains(true) => correct(r1, r2)
/** Consensus calls the overlapping bases if and only if the template is a paired end where both ends map with at
* least one base overlapping.
*
* @param template the template to potentially correct.
* @return summary statistics about how many bases were examined and modified
*/
def call(template: Template): CorrectionStats = (template.r1, template.r2) match {
case (Some(r1), Some(r2)) if r1.mapped && r2.mapped && r1.matesOverlap.contains(true) => call(r1, r2)
case _ => CorrectionStats(0, 0, 0, 0)
}

def correct(r1: SamRecord, r2: SamRecord): CorrectionStats = {
/** Consensus calls the overlapping bases if and only if both ends map with at least one base overlapping.
*
* @param r1 the first read in the pair
* @param r2 the second read in the pair
* @return summary statistics about how many bases were examined and modified
*/
def call(r1: SamRecord, r2: SamRecord): CorrectionStats = {
require(r1.mapped && r2.mapped && r1.paired && r2.paired && r1.name == r2.name && r1.matesOverlap.contains(true))

// Clear and resize the builders
r1Builders.foreach { builder =>
builder.clear()
Expand Down Expand Up @@ -91,7 +105,7 @@ class OverlappingBasesConsensusCaller(onlyMaskDisagreements: Boolean = false,
r2BasesBuilder.addAll(r2Bases.take(r2LastReadPos))
r2QualsBuilder.addAll(r2Quals.take(r2LastReadPos))

// Walk through the iterators by reference
// Walk through the iterators by reference position
while (r1Iter.hasNext && r2Iter.hasNext) {
val r1Head = r1Iter.head
val r2Head = r2Iter.head
Expand Down Expand Up @@ -172,7 +186,7 @@ class OverlappingBasesConsensusCaller(onlyMaskDisagreements: Boolean = false,
// Capture the raw consensus base prior to masking it to N, so that we can compute
// errors vs. the actually called base.
val (rawBase: Byte, rawQual: PhredScore) = {
if (base1 == base2 && maxQualOnAgreement) (base1, Math.max(qual1, qual2).toByte) // use the maximum base quality
if (base1 == base2 && maxQualOnAgreement) (base1, Math.max(qual1, qual2).toByte) // use the maximum base quality
else if (base1 == base2 && !maxQualOnAgreement) (base1, PhredScore.cap(qual1 + qual2)) // use the sum of base qualities
else if (onlyMaskDisagreements) (NoCall, NoCallQual) // disagreements are no-calls
else if (qual1 > qual2) (base1, PhredScore.cap(qual1 - qual2))
Expand All @@ -185,5 +199,12 @@ class OverlappingBasesConsensusCaller(onlyMaskDisagreements: Boolean = false,
}
}

/** Statistics for consensus calling overlapping bases in a read pair
*
* @param r1OverlappingBases the number of bases in R1 that overlap R2
* @param r2OverlappingBases the number of bases in R2 that overlap R1
* @param r1CorrectedBases the number of bases modified in R1
* @param r2CorrectedBases the number of bases modified in R2
*/
case class CorrectionStats(r1OverlappingBases: Int, r2OverlappingBases: Int, r1CorrectedBases: Int, r2CorrectedBases: Int)

0 comments on commit 5415287

Please sign in to comment.