Fix for mdtag issues with insertions

adam-core/src/main/scala/org/bdgenomics/adam/rich/RichAlignmentRecord.scala

@@ -89,7 +89,7 @@ class RichAlignmentRecord(val record: AlignmentRecord) {
   }
 
   lazy val samtoolsCigar: Cigar = {
-    TextCigarCodec.decode(record.getCigar.toString)
+    TextCigarCodec.decode(record.getCigar)
   }
 
   // Returns the MdTag if the read is mapped, None otherwise

adam-core/src/main/scala/org/bdgenomics/adam/util/MdTag.scala

@@ -57,8 +57,8 @@ object MdTag {
     var cigarIdx = 0
     var mdTagStringOffset = 0
     var referencePos = referenceStart
-    var cigarReferencePosition = referenceStart
-
+    var cigarOperatorIndex = 0
+    var usedMatchingBases = 0
     if (mdTagInput == null || mdTagInput == "0") {
       new MdTag(referenceStart, List(), Map(), Map())
     } else {
@@ -68,74 +68,99 @@ object MdTag {
         val cigarElement = cigar.getCigarElement(cigarIdx)
         val nextDigit = digitPattern.findPrefixOf(mdTag.substring(mdTagStringOffset))
         (cigarElement.getOperator, nextDigit) match {
-          case (_, None) if mdTagStringOffset == 0 => throw new IllegalArgumentException(s"MdTag $mdTagInput does not start with a digit")
-          case (_, Some(matchingBases)) if matchingBases.toInt == 0 => {
+
+          case (_, None) if mdTagStringOffset == 0 =>
+            throw new IllegalArgumentException(s"MdTag $mdTagInput does not start with a digit")
+
+          case (_, Some(matchingBases)) if matchingBases.toInt == 0 =>
             mdTagStringOffset += 1
-          }
-          case (CigarOperator.MATCH_OR_MISMATCH | CigarOperator.M | CigarOperator.EQ, Some(matchingBases)) => {
-            val numMatchingBases = math.min(cigarElement.getLength, matchingBases.toInt)
+
+          case (CigarOperator.M | CigarOperator.EQ, Some(matchingBases)) =>
+            val numMatchingBases = math.min(cigarElement.getLength - cigarOperatorIndex, matchingBases.toInt - usedMatchingBases)
+
             if (numMatchingBases > 0) {
               matches ::= NumericRange(referencePos, referencePos + numMatchingBases, 1L)
+
+              // Move the reference position the length of the matching sequence
+              referencePos += numMatchingBases
+
+              // Move ahead in the current CIGAR operator
+              cigarOperatorIndex += numMatchingBases
+
+              // Move ahead in the current MdTag digit
+              usedMatchingBases += numMatchingBases
+            }
+
+            if (matchingBases.toInt == usedMatchingBases) {
+              mdTagStringOffset += matchingBases.length
+              usedMatchingBases = 0
             }
-            // Move the reference position the length of the matching sequence
-            referencePos += numMatchingBases.toInt
-            if (matchingBases.toInt <= cigarElement.getLength) mdTagStringOffset += matchingBases.size
 
             // If the M operator has been fully read move on to the next operator
-            if (referencePos >= cigarReferencePosition + cigarElement.getLength) {
+            if (cigarOperatorIndex == cigarElement.getLength) {
               cigarIdx += 1
-              cigarReferencePosition += cigarElement.getLength
+              cigarOperatorIndex = 0
             }
-          }
-          case (CigarOperator.MATCH_OR_MISMATCH | CigarOperator.M | CigarOperator.X, None) => {
+
+          case (CigarOperator.M | CigarOperator.X, None) =>
             basesPattern.findPrefixOf(mdTag.substring(mdTagStringOffset)) match {
               // Must have found a base or digit pattern in CigarOperator M
-              case None => throw new IllegalArgumentException(s"No match or mismatched bases found for ${cigar.toString} in MDTag $mdTag")
-              case Some(mismatchedBases) => {
+              case None =>
+                throw new IllegalArgumentException(
+                  s"No match or mismatched bases found for ${cigar.toString} in MDTag $mdTag"
+                )
+              case Some(mismatchedBases) =>
                 mismatchedBases.foreach {
                   base =>
                     mismatches += (referencePos -> base)
                     referencePos += 1
                 }
-                mdTagStringOffset += mismatchedBases.size
-              }
+                mdTagStringOffset += mismatchedBases.length
+                cigarOperatorIndex += mismatchedBases.length
             }
             // If the M operator has been fully read move on to the next operator
-            if (referencePos >= cigarReferencePosition + cigarElement.getLength) {
+            if (cigarOperatorIndex == cigarElement.getLength) {
               cigarIdx += 1
-              cigarReferencePosition += cigarElement.getLength
+              cigarOperatorIndex = 0
             }
-          }
-          case (CigarOperator.DELETION, None) => {
+
+          case (CigarOperator.DELETION, None) =>
             mdTag.charAt(mdTagStringOffset) match {
-              case '^' => {
+              case '^' =>
                 // Skip ahead of the deletion '^' character
                 mdTagStringOffset += 1
                 basesPattern.findPrefixOf(mdTag.substring(mdTagStringOffset)) match {
-                  case None => throw new IllegalArgumentException(s"No deleted bases found ${cigar.toString} in MDTag $mdTag")
-                  case Some(deletedBases) => {
+                  case None =>
+                    throw new IllegalArgumentException(s"No deleted bases found ${cigar.toString} in MDTag $mdTag")
+                  case Some(deletedBases) if deletedBases.length == cigarElement.getLength =>
                     deletedBases.foreach {
                       base =>
                         deletions += (referencePos -> base)
                         referencePos += 1
                     }
-                    mdTagStringOffset += deletedBases.size
-                  }
+                    mdTagStringOffset += deletedBases.length
+                  case Some(deletedBases) =>
+                    throw new IllegalArgumentException(
+                      s"Element ${cigarElement.getLength}${cigarElement.getOperator.toString} in cigar ${cigar.toString} contradicts number of bases listed in MDTag: ^${deletedBases}"
+                    )
                 }
-                cigarReferencePosition += cigarElement.getLength
                 cigarIdx += 1
-              }
-              case _ => throw new IllegalArgumentException(s"CIGAR ${cigar.toString} indicates deletion found but no deleted bases in MDTag $mdTagInput")
+                cigarOperatorIndex = 0
+
+              case _ =>
+                throw new IllegalArgumentException(
+                  s"CIGAR ${cigar.toString} indicates deletion found but no deleted bases in MDTag $mdTagInput"
+                )
             }
-          }
-          case _ if cigarElement.getOperator.consumesReferenceBases() => {
+
+          case (CigarOperator.N, _) =>
             referencePos += cigarElement.getLength
-            cigarReferencePosition += cigarElement.getLength
             cigarIdx += 1
-          }
-          case _ => {
+            cigarOperatorIndex = 0
+
+          case (CigarOperator.INSERTION | CigarOperator.H | CigarOperator.S, _) =>
             cigarIdx += 1
-          }
+
         }
       }
       new MdTag(referenceStart, matches, mismatches, deletions)
@@ -400,7 +425,7 @@ class MdTag(
    * @return True if this read has mismatches. We do not return true if the read has no mismatches but has deletions.
    */
   def hasMismatches: Boolean = {
-    !mismatches.isEmpty
+    mismatches.nonEmpty
   }
 
   /**
@@ -419,7 +444,7 @@ class MdTag(
    */
   def end(): Long = {
     val ends = matches.map(_.end - 1) ::: mismatches.keys.toList ::: deletions.keys.toList
-    ends.reduce(_ max _)
+    ends.max
   }
 
   /**
@@ -496,7 +521,7 @@ class MdTag(
    * @return MD string corresponding to [0-9]+(([A-Z]|\&#94;[A-Z]+)[0-9]+)
    * @see http://zenfractal.com/2013/06/19/playing-with-matches/
    */
-  override def toString(): String = {
+  override def toString: String = {
     if (matches.isEmpty && mismatches.isEmpty && deletions.isEmpty) {
       "0"
     } else {

adam-core/src/test/scala/org/bdgenomics/adam/util/MdTagSuite.scala

@@ -226,6 +226,138 @@ class MdTagSuite extends FunSuite {
     assert(tag.toString === "20")
   }
 
+  test("Get correct matches for mdtag with insertion") {
+    val tag = MdTag("10", 0L, TextCigarCodec.decode("5M3I5M"))
+    assert(tag.end === 9)
+
+    (0 until 9).foreach(locus => assert(tag.isMatch(locus)))
+
+    assert(tag.toString === "10")
+  }
+
+  test("Get correct matches for mdtag with mismatches and insertion") {
+    val tag = MdTag("2A7", 0L, TextCigarCodec.decode("5M3I5M"))
+    assert(tag.end === 9)
+
+    assert(tag.isMatch(0))
+    assert(tag.isMatch(1))
+    assert(tag.mismatches(2) === 'A')
+    (3 until 9).foreach(locus => assert(tag.isMatch(locus)))
+
+    assert(tag.toString === "2A7")
+  }
+
+  test("Get correct matches for mdtag with insertion between mismatches") {
+    val tag = MdTag("2A4A2", 0L, TextCigarCodec.decode("5M3I5M"))
+    assert(tag.end === 9)
+
+    assert(tag.isMatch(0))
+    assert(tag.isMatch(1))
+    assert(tag.mismatches(2L) === 'A')
+
+    assert(tag.mismatches(7L) === 'A')
+    assert(tag.isMatch(8))
+    assert(tag.isMatch(9))
+
+    assert(tag.toString === "2A4A2")
+  }
+
+  test("Get correct matches for mdtag with intron between mismatches") {
+    val tag = MdTag("2A4A2", 0L, TextCigarCodec.decode("5M3N5M"))
+    assert(tag.end === 12)
+
+    assert(tag.isMatch(0))
+    assert(tag.isMatch(1))
+    assert(tag.mismatches(2L) === 'A')
+    assert(tag.isMatch(3))
+    assert(tag.isMatch(4))
+
+    assert(!tag.isMatch(5))
+    assert(!tag.isMatch(6))
+    assert(!tag.isMatch(7))
+
+    assert(tag.isMatch(8))
+    assert(tag.isMatch(9))
+
+    assert(tag.mismatches(10L) === 'A')
+    assert(tag.isMatch(11))
+    assert(tag.isMatch(12))
+
+    assert(tag.toString === "2A4A2")
+  }
+
+  test("Get correct matches for mdtag with intron and deletion between mismatches") {
+    val tag = MdTag("2A4A0^AAA2", 0L, TextCigarCodec.decode("5M3N3M3D2M"))
+    assert(tag.end === 15)
+
+    assert(tag.isMatch(0))
+    assert(tag.isMatch(1))
+    assert(tag.mismatches(2L) === 'A')
+    assert(tag.isMatch(3))
+    assert(tag.isMatch(4))
+
+    assert(!tag.isMatch(5))
+    assert(!tag.isMatch(6))
+    assert(!tag.isMatch(7))
+
+    assert(tag.isMatch(8))
+    assert(tag.isMatch(9))
+
+    assert(tag.mismatches(10L) === 'A')
+
+    assert(tag.deletions(11L) === 'A')
+    assert(tag.deletions(12L) === 'A')
+    assert(tag.deletions(13L) === 'A')
+
+    assert(tag.toString === "2A4A0^AAA2")
+  }
+
+  test("Throw exception when number of deleted bases in mdtag disagrees with CIGAR") {
+    intercept[IllegalArgumentException] {
+      MdTag("2A4A0^AAA2", 0L, TextCigarCodec.decode("5M3N3M4D2M"))
+    }
+  }
+
+  test("Get correct matches for mdtag with mismatch, insertion and deletion") {
+    val tag = MdTag("2A3^AAA4", 0L, TextCigarCodec.decode("5M3I1M3D4M"))
+    assert(tag.end === 12)
+
+    assert(tag.isMatch(0))
+    assert(tag.isMatch(1))
+    assert(tag.mismatches(2L) === 'A')
+
+    (3 to 5).foreach(locus => assert(tag.isMatch(locus)))
+
+    assert(tag.deletedBase(6L) === Some('A'))
+    assert(tag.deletedBase(7L) === Some('A'))
+    assert(tag.deletedBase(8L) === Some('A'))
+
+    (9 to 12).foreach(locus => assert(tag.isMatch(locus)))
+
+    assert(tag.toString === "2A3^AAA4")
+  }
+
+  test("Get correct matches for mdtag with mismatches, insertion and deletion") {
+    val tag = MdTag("2A3^AAA2A1", 0L, TextCigarCodec.decode("5M3I1M3D4M"))
+    assert(tag.end === 12)
+
+    assert(tag.isMatch(0))
+    assert(tag.isMatch(1))
+    assert(tag.mismatches(2L) === 'A')
+
+    (3 to 5).foreach(locus => assert(tag.isMatch(locus)))
+
+    assert(tag.deletedBase(6L) === Some('A'))
+    assert(tag.deletedBase(7L) === Some('A'))
+    assert(tag.deletedBase(8L) === Some('A'))
+
+    (9 to 10).foreach(locus => assert(tag.isMatch(locus)))
+
+    assert(tag.mismatches(11L) === 'A')
+
+    assert(tag.toString === "2A3^AAA2A1")
+  }
+
   test("Get correct matches for MDTag with mismatches and deletions") {
 
     val tag1 = MdTag("40A5^TTT54", 0L, TextCigarCodec.decode("46M3D54M"))
@@ -237,6 +369,7 @@ class MdTagSuite extends FunSuite {
     assert(tag1.isMatch(41) === true)
 
     assert(tag1.mismatchedBase(40) === Some('A'))
+    assert(tag1.toString === "40A5^TTT54")
 
     val tag2 = MdTag("40A5^TTT0G53", 0L, TextCigarCodec.decode("46M3D54M"))
     assert(tag2.hasMismatches === true)
@@ -249,6 +382,22 @@ class MdTagSuite extends FunSuite {
     assert(tag2.mismatchedBase(40) === Some('A'))
     assert(tag2.mismatchedBase(49) === Some('G'))
     assert(tag2.isMatch(50) === true)
+    assert(tag2.toString === "40A5^TTT0G53")
+
+    val tag3 = MdTag("2^GA5^TC6", 0L, TextCigarCodec.decode("2M2D1M2I2M4I2M2D6M"))
+    (0 to 1).foreach(l => assert(tag3.isMatch(l)))
+    (2 to 3).foreach(l => assert(!tag3.isMatch(l)))
+    (4 to 8).foreach(l => assert(tag3.isMatch(l)))
+    (9 to 10).foreach(l => assert(!tag3.isMatch(l)))
+    (11 to 16).foreach(l => assert(tag3.isMatch(l)))
+
+    assert(tag3.deletedBase(2) == Some('G'))
+    assert(tag3.deletedBase(3) == Some('A'))
+
+    assert(tag3.deletedBase(9) == Some('T'))
+    assert(tag3.deletedBase(10) == Some('C'))
+
+    assert(tag3.toString === "2^GA5^TC6")
   }
 
   test("Get correct matches base from MDTag and CIGAR with N") {
@@ -453,4 +602,4 @@ class MdTagSuite extends FunSuite {
   test("handle '=' and 'X' operators") {
     testTag("ACCCAAGT", "ACCATAGA", "3=2X2=1X", 0, "3A0T2A0", 0, 7)
   }
-}
+}