Misc. INDEL realigner bugfixes

adam-core/src/main/scala/org/bdgenomics/adam/algorithms/consensus/ConsensusGenerator.scala

@@ -83,10 +83,13 @@ object ConsensusGenerator {
    * Provides a generator to extract consensuses from a known set of INDELs.
    *
    * @param rdd The previously called INDEL variants.
+   * @param flankSize The number of bases to flank each known INDEL by. Default
+   *   is 0 bases.
    * @return A consensus generator that looks at previously called INDELs.
    */
-  def fromKnownIndels(rdd: VariantRDD): ConsensusGenerator = {
-    new ConsensusGeneratorFromKnowns(rdd.rdd)
+  def fromKnownIndels(rdd: VariantRDD,
+                      flankSize: Int = 0): ConsensusGenerator = {
+    new ConsensusGeneratorFromKnowns(rdd.rdd, flankSize)
   }
 
   /**

adam-core/src/main/scala/org/bdgenomics/adam/algorithms/consensus/ConsensusGeneratorFromKnowns.scala

@@ -24,6 +24,7 @@ import org.bdgenomics.adam.rdd.ADAMContext._
 import org.bdgenomics.adam.rdd.read.realignment.IndelRealignmentTarget
 import org.bdgenomics.adam.rich.RichAlignmentRecord
 import org.bdgenomics.formats.avro.Variant
+import scala.math.max
 import scala.transient
 
 /**
@@ -36,7 +37,8 @@ import scala.transient
  * @param file Path to file containing variants.
  * @param sc Spark context to use.
  */
-private[adam] class ConsensusGeneratorFromKnowns(rdd: RDD[Variant]) extends ConsensusGenerator {
+private[adam] class ConsensusGeneratorFromKnowns(rdd: RDD[Variant],
+                                                 val flankSize: Int) extends ConsensusGenerator {
 
   private val indelTable = rdd.context.broadcast(IndelTable(rdd))
 
@@ -50,7 +52,10 @@ private[adam] class ConsensusGeneratorFromKnowns(rdd: RDD[Variant]) extends Cons
 
     Some(rdd.filter(v => v.getReferenceAllele.length != v.getAlternateAllele.length)
       .map(v => ReferenceRegion(v))
-      .map(r => new IndelRealignmentTarget(Some(r), r)))
+      .map(r => {
+        new IndelRealignmentTarget(Some(r),
+          ReferenceRegion(r.referenceName, max(0L, r.start - flankSize), r.end + flankSize))
+      }))
   }
 
   /**

adam-core/src/main/scala/org/bdgenomics/adam/rdd/read/realignment/RealignIndels.scala

@@ -284,6 +284,7 @@ private[read] class RealignIndels(
 
         // loop over all consensuses and evaluate
         consensus.foreach(c => {
+
           // generate a reference sequence from the consensus
           val consensusSequence = c.insertIntoReference(reference, refRegion)
 
@@ -294,7 +295,7 @@ private[read] class RealignIndels(
 
             // if the read's mismatch quality improves over the original alignment, save
             // its alignment in the consensus sequence, else store -1
-            if (qual < originalQual) {
+            if (qual <= originalQual) {
               (r, (qual, pos))
             } else {
               (r, (originalQual, -1))
@@ -359,23 +360,41 @@ private[read] class RealignIndels(
                 // compensate the end
                 builder.setEnd(refStart + remapping + r.getSequence.length + endPenalty)
 
+                val startLength = bestConsensus.index.start - (refStart + remapping)
+                val adjustedIdElement = if (endLength < 0) {
+                  new CigarElement(idElement.getLength + endLength.toInt, idElement.getOperator)
+                } else if (startLength < 0) {
+                  new CigarElement(idElement.getLength + startLength.toInt, idElement.getOperator)
+                } else {
+                  idElement
+                }
+
                 val cigarElements = List[CigarElement](
-                  new CigarElement((bestConsensus.index.start - (refStart + remapping)).toInt, CigarOperator.M),
-                  idElement,
+                  new CigarElement(startLength.toInt, CigarOperator.M),
+                  adjustedIdElement,
                   new CigarElement(endLength.toInt, CigarOperator.M)
-                )
+                ).filter(_.getLength > 0)
 
                 new Cigar(cigarElements)
               }
 
               // update mdtag and cigar
-              builder.setMismatchingPositions(MdTag.moveAlignment(r, newCigar, reference.drop(remapping), refStart + remapping).toString())
-              builder.setOldPosition(r.getStart())
-              builder.setOldCigar(r.getCigar())
-              builder.setCigar(newCigar.toString)
-              new RichAlignmentRecord(builder.build())
-            } else
+              try {
+                builder.setMismatchingPositions(MdTag.moveAlignment(r, newCigar, reference.drop(remapping), refStart + remapping).toString())
+                builder.setOldPosition(r.getStart())
+                builder.setOldCigar(r.getCigar())
+                builder.setCigar(newCigar.toString)
+                new RichAlignmentRecord(builder.build())
+              } catch {
+                case t: Throwable => {
+                  log.warn("Caught %s when trying to move alignment to %s for %s.".format(
+                    t, newCigar, r))
+                  new RichAlignmentRecord(r)
+                }
+              }
+            } else {
               new RichAlignmentRecord(builder.build())
+            }
           })
 
           log.info("On " + refRegion + ", realigned " + realignedReadCount + " reads to " +