Pass assembly variants to realigner as knowns

avocado-cli/src/main/scala/org/bdgenomics/avocado/cli/Reassemble.scala

@@ -20,8 +20,7 @@ package org.bdgenomics.avocado.cli
 import org.apache.spark.SparkContext
 import org.bdgenomics.adam.rdd.ADAMContext._
 import org.bdgenomics.adam.rdd.ADAMSaveAnyArgs
-import org.bdgenomics.adam.rdd.read.{ AlignedReadRDD, AlignmentRecordRDD, MDTagging }
-import org.bdgenomics.avocado.realigner.Realigner
+import org.bdgenomics.avocado.realigner.{ ConsensusRealigner, Realigner }
 import org.bdgenomics.utils.cli._
 import org.kohsuke.args4j.{ Argument, Option => Args4jOption }
 
@@ -58,6 +57,10 @@ class ReassembleArgs extends Args4jBase with ADAMSaveAnyArgs with ParquetArgs {
     name = "-defer_merging",
     usage = "Defers merging single file output")
   var deferMerging: Boolean = false
+  @Args4jOption(required = false,
+    name = "-use_consensus_realigner",
+    usage = "If provided, uses consensus-mode realigner.")
+  var useConsensusRealigner: Boolean = false
 
   // required by ADAMSaveAnyArgs
   var sortFastqOutput: Boolean = false
@@ -79,7 +82,11 @@ class Reassemble(
     val reads = sc.loadAlignments(args.inputPath)
 
     // realign the reads
-    val reassembledReads = Realigner.realign(reads, args.kmerLength)
+    val reassembledReads = if (args.useConsensusRealigner) {
+      ConsensusRealigner.realign(reads, args.kmerLength)
+    } else {
+      Realigner.realign(reads, args.kmerLength)
+    }
 
     // save the reads
     reassembledReads.save(args)

avocado-core/src/main/scala/org/bdgenomics/avocado/genotyping/BiallelicGenotyper.scala

@@ -222,10 +222,26 @@ private[avocado] object BiallelicGenotyper extends Serializable with Logging {
             } else if (isInsertion(variant)) {
               val insAllele = variant.getAlternateAllele.tail
               val insObserved = observed.filter(_._2 == insAllele)
+              println("for %s observed %s in %s".format(insAllele, observed.mkString(","), read))
               if (observed.size == 2 &&
                 insObserved.size == 1) {
                 Some((variant, insObserved.head._3.duplicate(Some(false))))
-              } else if (observed.forall(_._3.isRef)) {
+              } else if (!observed.forall(_._3.isRef)) {
+                val nonRef = observed.filter(!_._3.isRef)
+                if (nonRef.size == 1) {
+                  val (_, allele, nonRefObs) = nonRef.head
+                  if (allele.length == insAllele.length) {
+                    val matchingBases = allele.zip(insAllele).count(p => p._1 == p._2)
+                    Some((variant, nonRefObs.scale(matchingBases,
+                      insAllele.length,
+                      Some(false))))
+                  } else {
+                    Some((variant, observed.head._3.nullOut))
+                  }
+                } else {
+                  Some((variant, observed.head._3.nullOut))
+                }
+              } else if (read.getEnd != variant.getEnd) {
                 Some((variant, observed.head._3.invert))
               } else {
                 Some((variant, observed.head._3.nullOut))

avocado-core/src/main/scala/org/bdgenomics/avocado/genotyping/DiscoverVariants.scala

@@ -201,7 +201,11 @@ object DiscoverVariants extends Serializable with Logging {
               kv
             }
             case Insertion(length) => {
-              val insQuals = qual.substring(idx - 1, idx + length).map(_.toInt - 33).sum / length
+              val insQuals = if (length > 0) {
+                qual.substring(idx - 1, idx + length).map(_.toInt - 33).sum / length
+              } else {
+                0
+              }
               val newVar = if (insQuals >= phredThreshold) {
                 Variant.newBuilder
                   .setContigName(contigName)

avocado-core/src/main/scala/org/bdgenomics/avocado/genotyping/Observer.scala

@@ -131,9 +131,13 @@ private[genotyping] object Observer extends Serializable {
 
         // get bases and quals
         val bases = readSequence.substring(oldReadIdx, readIdx)
-        val qual = readQualities.substring(oldReadIdx, readIdx)
-          .map(_.toInt - 33)
-          .sum / length
+        val qual = if (length > 0) {
+          readQualities.substring(oldReadIdx, readIdx)
+            .map(_.toInt - 33)
+            .sum / length
+        } else {
+          0
+        }
 
         // the key is the (site, allele, sampleId)
         // insertions associate to the site to their left, hence the -1

avocado-core/src/main/scala/org/bdgenomics/avocado/models/Observation.scala

@@ -91,6 +91,24 @@ case class Observation(alleleForwardStrand: Int,
       isRef = setRef.getOrElse(isRef))
   }
 
+  /**
+   * @return Makes a copy where underlying arrays are not shared.
+   */
+  def scale(num: Int,
+            denum: Int,
+            setRef: Option[Boolean] = None): Observation = {
+    val scaleBy = num.toDouble / denum.toDouble
+    Observation(alleleForwardStrand,
+      otherForwardStrand,
+      squareMapQ,
+      alleleLogLikelihoods.map(v => v * scaleBy),
+      otherLogLikelihoods.map(v => v * scaleBy),
+      alleleCoverage,
+      otherCoverage,
+      totalCoverage = totalCoverage,
+      isRef = setRef.getOrElse(isRef))
+  }
+
   /**
    * @return Returns this observation, but with allele/other swapped.
    *

avocado-core/src/main/scala/org/bdgenomics/avocado/realigner/ConsensusRealigner.scala

@@ -0,0 +1,64 @@
+/**
+ * Licensed to Big Data Genomics (BDG) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The BDG licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.bdgenomics.avocado.realigner
+
+import org.bdgenomics.adam.algorithms.consensus.ConsensusGenerator
+import org.bdgenomics.adam.rdd.read.AlignmentRecordRDD
+import org.bdgenomics.adam.rdd.variant.VariantRDD
+import org.bdgenomics.avocado.genotyping.DiscoverVariants
+import org.bdgenomics.formats.avro.Variant
+
+object ConsensusRealigner {
+
+  /**
+   * Realigns a set of reads against the reference genome.
+   *
+   * @param reads Reads to realign.
+   * @param kmerLength The length k of the k-mers.
+   * @return Returns the realigned reads.
+   */
+  def realign(reads: AlignmentRecordRDD,
+              kmerLength: Int): AlignmentRecordRDD = {
+
+    // use the realigner to make a first pass over the reads
+    val realignedReads = Realigner.realign(reads, kmerLength)
+
+    // from here, we'll discover any potential variants
+    val variants = filterIndels(DiscoverVariants(realignedReads))
+
+    // we'll pass these discovered variants to ADAM's indel realigner
+    reads.realignIndels(ConsensusGenerator.fromKnownIndels(variants))
+  }
+
+  /**
+   * @param v The variant to filter.
+   * @return Returns false if the variant is a SNV/MNV.
+   */
+  private[realigner] def discardSnvs(v: Variant): Boolean = {
+    v.getReferenceAllele.length != v.getAlternateAllele.length
+  }
+
+  /**
+   * @param rdd The RDD of variants to filter.
+   * @return Returns a new RDD of variants where the SNVs & MNVs have been
+   *   removed and only INDEL variants remain.
+   */
+  private[realigner] def filterIndels(rdd: VariantRDD): VariantRDD = {
+    rdd.transform(r => r.filter(discardSnvs))
+  }
+}

avocado-core/src/main/scala/org/bdgenomics/avocado/util/PrefilterReads.scala

@@ -18,7 +18,7 @@
 package org.bdgenomics.avocado.util
 
 import org.bdgenomics.adam.models.SequenceDictionary
-import org.bdgenomics.adam.rdd.read.{ AlignedReadRDD, AlignmentRecordRDD }
+import org.bdgenomics.adam.rdd.read.AlignmentRecordRDD
 import org.bdgenomics.formats.avro.AlignmentRecord
 
 private[avocado] trait PrefilterReadsArgs extends Serializable {
@@ -85,7 +85,7 @@ private[avocado] object PrefilterReads extends Serializable {
     val readRdd = rdd.rdd.filter(readFn)
 
     // construct a new aligned read rdd
-    AlignedReadRDD(readRdd, sequences, rdd.recordGroups)
+    AlignmentRecordRDD(readRdd, sequences, rdd.recordGroups)
   }
 
   /**

avocado-core/src/test/scala/org/bdgenomics/avocado/genotyping/BiallelicGenotyperSuite.scala

@@ -20,6 +20,7 @@ package org.bdgenomics.avocado.genotyping
 import org.bdgenomics.adam.rdd.ADAMContext._
 import org.bdgenomics.avocado.AvocadoFunSuite
 import org.bdgenomics.avocado.models.Observation
+import org.bdgenomics.avocado.realigner.ConsensusRealigner
 import org.bdgenomics.avocado.util.{
   HardFilterGenotypes,
   TestHardFilterGenotypesArgs
@@ -437,14 +438,16 @@ class BiallelicGenotyperSuite extends AvocadoFunSuite {
     })
   }
 
-  ignore("call hom alt C->CCCCT insertion at 1/866511") {
+  sparkTest("call hom alt C->CCCCT insertion at 1/866511") {
     val readPath = resourceUrl("NA12878.chr1.866511.sam")
     val reads = sc.loadAlignments(readPath.toString)
       .transform(rdd => {
         rdd.filter(_.getMapq > 0)
       })
 
-    val gts = BiallelicGenotyper.discoverAndCall(reads,
+    val realignedReads = ConsensusRealigner.realign(reads, 20)
+
+    val gts = BiallelicGenotyper.discoverAndCall(realignedReads,
       2,
       optPhredThreshold = Some(30)).transform(rdd => {
         rdd.filter(gt => gt.getAlleles.forall(_ == GenotypeAllele.ALT))

avocado-core/src/test/scala/org/bdgenomics/avocado/genotyping/DiscoverVariantsSuite.scala

@@ -22,7 +22,7 @@ import org.bdgenomics.adam.models.{
   SequenceDictionary,
   SequenceRecord
 }
-import org.bdgenomics.adam.rdd.read.AlignedReadRDD
+import org.bdgenomics.adam.rdd.read.AlignmentRecordRDD
 import org.bdgenomics.avocado.AvocadoFunSuite
 import org.bdgenomics.formats.avro.{ AlignmentRecord, Variant }
 
@@ -233,7 +233,7 @@ class DiscoverVariantsSuite extends AvocadoFunSuite {
       snpReadMCigar, snpReadEqCigar,
       insertRead,
       deleteRead))
-    val readRdd = AlignedReadRDD(rdd,
+    val readRdd = AlignmentRecordRDD(rdd,
       SequenceDictionary(
         SequenceRecord("1", 50L),
         SequenceRecord("2", 40L),

avocado-core/src/test/scala/org/bdgenomics/avocado/realigner/ConsensusRealignerSuite.scala

@@ -0,0 +1,85 @@
+/**
+ * Licensed to Big Data Genomics (BDG) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The BDG licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.bdgenomics.avocado.realigner
+
+import org.bdgenomics.adam.models.SequenceDictionary
+import org.bdgenomics.adam.rdd.read.AlignmentRecordRDD
+import org.bdgenomics.adam.rdd.variant.VariantRDD
+import org.bdgenomics.formats.avro.Variant
+
+class ConsensusRealignerSuite extends SparkRealignerSuite {
+
+  val allowLegacyCigars = true
+
+  def realign(rdd: AlignmentRecordRDD,
+              kmerLength: Int): AlignmentRecordRDD = {
+    ConsensusRealigner.realign(rdd, kmerLength)
+  }
+
+  val snv = Variant.newBuilder
+    .setStart(0L)
+    .setReferenceAllele("A")
+    .setAlternateAllele("G")
+    .build
+  val mnv = Variant.newBuilder
+    .setStart(1L)
+    .setReferenceAllele("AC")
+    .setAlternateAllele("GT")
+    .build
+  val ins = Variant.newBuilder
+    .setStart(2L)
+    .setReferenceAllele("A")
+    .setAlternateAllele("ACC")
+    .build
+  val del = Variant.newBuilder
+    .setStart(3L)
+    .setReferenceAllele("GCT")
+    .setAlternateAllele("G")
+    .build
+
+  test("filter out a snv") {
+    assert(!ConsensusRealigner.discardSnvs(snv))
+  }
+
+  test("filter out a mnv") {
+    assert(!ConsensusRealigner.discardSnvs(mnv))
+  }
+
+  test("keep an insert") {
+    assert(ConsensusRealigner.discardSnvs(ins))
+  }
+
+  test("keep a deletion") {
+    assert(ConsensusRealigner.discardSnvs(del))
+  }
+
+  sparkTest("filter snv/mnv variants out") {
+    val vRdd = VariantRDD(sc.parallelize(Seq(snv, mnv, ins, del)),
+      SequenceDictionary.empty)
+    val filteredVariants = ConsensusRealigner.filterIndels(vRdd)
+      .rdd
+      .collect
+
+    val variantIds = filteredVariants.map(_.getStart)
+      .toSet
+
+    assert(variantIds.size === 2)
+    assert(variantIds(2L))
+    assert(variantIds(3L))
+  }
+}