Implement contig-to-RefSeq translation, refactor SequenceDictionary with...

adam-cli/src/main/scala/org/bdgenomics/adam/cli/ADAM2Vcf.scala

@@ -19,10 +19,13 @@ package org.bdgenomics.adam.cli
 import org.bdgenomics.adam.avro.ADAMGenotype
 import org.bdgenomics.adam.rdd.ADAMContext._
 import org.bdgenomics.adam.rdd.variation.ADAMVariationContext._
-import org.kohsuke.args4j.Argument
+import org.kohsuke.args4j.{ Option => Args4jOption, Argument }
 import org.apache.spark.rdd.RDD
 import org.apache.spark.{ Logging, SparkContext }
 import org.apache.hadoop.mapreduce.Job
+import java.io.File
+import org.bdgenomics.adam.models.SequenceDictionary
+import scala.Option
 
 object ADAM2Vcf extends ADAMCommandCompanion {
 
@@ -35,17 +38,25 @@ object ADAM2Vcf extends ADAMCommandCompanion {
 }
 
 class ADAM2VcfArgs extends Args4jBase with ParquetArgs with SparkArgs {
+  @Args4jOption(required = false, name = "-dict", usage = "Reference dictionary")
+  var dictionaryFile: File = _
+
   @Argument(required = true, metaVar = "ADAM", usage = "The ADAM variant files to convert", index = 0)
   var adamFile: String = _
+
   @Argument(required = true, metaVar = "VCF", usage = "Location to write VCF data", index = 1)
   var outputPath: String = null
 }
 
-class ADAM2Vcf(val args: ADAM2VcfArgs) extends ADAMSparkCommand[ADAM2VcfArgs] with Logging {
+class ADAM2Vcf(val args: ADAM2VcfArgs) extends ADAMSparkCommand[ADAM2VcfArgs] with DictionaryCommand with Logging {
   val companion = ADAM2Vcf
 
   def run(sc: SparkContext, job: Job) {
+    var dictionary: Option[SequenceDictionary] = loadSequenceDictionary(args.dictionaryFile)
+    if (dictionary.isDefined)
+      log.info("Using contig translation")
+
     val adamGTs: RDD[ADAMGenotype] = sc.adamLoad(args.adamFile)
-    sc.adamVCFSave(args.outputPath, adamGTs.toADAMVariantContext)
+    sc.adamVCFSave(args.outputPath, adamGTs.toADAMVariantContext, dict = dictionary)
   }
 }

adam-cli/src/main/scala/org/bdgenomics/adam/cli/Bam2ADAM.scala

@@ -102,7 +102,7 @@ class Bam2ADAM(args: Bam2ADAMArgs) extends ADAMCommand {
     val samReader = new SAMFileReader(new File(args.bamFile), null, true)
     samReader.setValidationStringency(args.validationStringency)
 
-    val seqDict = SequenceDictionary.fromSAMReader(samReader)
+    val seqDict = SequenceDictionary(samReader)
     val rgDict = RecordGroupDictionary.fromSAMReader(samReader)
 
     println(seqDict)

adam-cli/src/main/scala/org/bdgenomics/adam/cli/DictionaryCommand.scala

@@ -0,0 +1,47 @@
+/*
+ * Copyright (c) 2014. Mount Sinai School of Medicine
+ *
+ * Licensed 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.adam.cli
+
+import java.io.{ FileOutputStream, File }
+import org.apache.commons.io.IOUtils
+import org.bdgenomics.adam.models.SequenceDictionary
+import net.sf.samtools.SAMFileReader
+
+trait DictionaryCommand {
+  private def getDictionaryFile(name: String): Option[File] = {
+    val stream = ClassLoader.getSystemClassLoader.getResourceAsStream("dictionaries/" + name)
+    if (stream == null)
+      return None
+    val file = File.createTempFile(name, ".dict")
+    file.deleteOnExit()
+    IOUtils.copy(stream, new FileOutputStream(file))
+    Some(file)
+  }
+
+  private def getDictionary(file: File) = Some(SequenceDictionary(SAMFileReader.getSequenceDictionary(file)))
+
+  def loadSequenceDictionary(file: File): Option[SequenceDictionary] = {
+    if (file != null) {
+      if (file.exists)
+        getDictionary(file)
+      else getDictionaryFile(file.getName) match {
+        case Some(file) => getDictionary(file)
+        case _          => None
+      }
+    } else None
+  }
+}

adam-cli/src/main/scala/org/bdgenomics/adam/cli/ListDict.scala

@@ -47,7 +47,7 @@ class ListDict(protected val args: ListDictArgs) extends ADAMSparkCommand[ListDi
 
     val dict = sc.adamDictionaryLoad[ADAMRecord](args.inputPath)
 
-    dict.recordsIn.sortBy(_.name.toString).foreach {
+    dict.records.foreach {
       rec: SequenceRecord =>
         println("%s\t%d".format(rec.name, rec.length))
     }

adam-cli/src/main/scala/org/bdgenomics/adam/cli/Vcf2ADAM.scala

@@ -16,13 +16,14 @@
 
 package org.bdgenomics.adam.cli
 
-import org.bdgenomics.adam.models.ADAMVariantContext
+import org.bdgenomics.adam.models.{ SequenceDictionary, ADAMVariantContext }
 import org.bdgenomics.adam.rdd.ADAMContext._
 import org.bdgenomics.adam.rdd.variation.ADAMVariationContext._
 import org.apache.hadoop.mapreduce.Job
-import org.apache.spark.SparkContext
+import org.apache.spark.{ Logging, SparkContext }
 import org.apache.spark.rdd.RDD
-import org.kohsuke.args4j.Argument
+import org.kohsuke.args4j.{ Option => Args4jOption, Argument }
+import java.io.File
 
 object Vcf2ADAM extends ADAMCommandCompanion {
   val commandName = "vcf2adam"
@@ -34,21 +35,31 @@ object Vcf2ADAM extends ADAMCommandCompanion {
 }
 
 class Vcf2ADAMArgs extends Args4jBase with ParquetArgs with SparkArgs {
-  @Argument(required = true, metaVar = "VCF",
-    usage = "The VCF file to convert", index = 0)
-  var vcfFile: String = _
-  @Argument(required = true, metaVar = "ADAM",
-    usage = "Location to write ADAM Variant data", index = 1)
+  @Args4jOption(required = false, name = "-dict", usage = "Reference dictionary")
+  var dictionaryFile: File = _
+
+  @Argument(required = true, metaVar = "VCF", usage = "The VCF file to convert", index = 0)
+  var vcfPath: String = _
+
+  @Argument(required = true, metaVar = "ADAM", usage = "Location to write ADAM Variant data", index = 1)
   var outputPath: String = null
 }
 
-class Vcf2ADAM(val args: Vcf2ADAMArgs) extends ADAMSparkCommand[Vcf2ADAMArgs] {
+class Vcf2ADAM(val args: Vcf2ADAMArgs) extends ADAMSparkCommand[Vcf2ADAMArgs] with DictionaryCommand with Logging {
   val companion = Vcf2ADAM
 
   def run(sc: SparkContext, job: Job) {
-    var adamVariants: RDD[ADAMVariantContext] = sc.adamVCFLoad(args.vcfFile)
-    adamVariants.flatMap(p => p.genotypes).adamSave(args.outputPath,
-      blockSize = args.blockSize, pageSize = args.pageSize,
+
+    var dictionary: Option[SequenceDictionary] = loadSequenceDictionary(args.dictionaryFile)
+    if (dictionary.isDefined)
+      log.info("Using contig translation")
+
+    var adamVariants: RDD[ADAMVariantContext] = sc.adamVCFLoad(args.vcfPath, dict = dictionary)
+
+    adamVariants.flatMap(p => p.genotypes).adamSave(
+      args.outputPath,
+      blockSize = args.blockSize,
+      pageSize = args.pageSize,
       compressCodec = args.compressionCodec,
       disableDictionaryEncoding = args.disableDictionary)
   }

adam-core/src/main/scala/org/bdgenomics/adam/converters/SAMRecordConverter.scala

@@ -17,9 +17,9 @@ package org.bdgenomics.adam.converters
 
 import net.sf.samtools.{ SAMReadGroupRecord, SAMRecord }
 
-import org.bdgenomics.adam.avro.{ ADAMContig, ADAMRecord }
+import org.bdgenomics.adam.avro.ADAMRecord
 import scala.collection.JavaConverters._
-import org.bdgenomics.adam.models.{ Attribute, RecordGroupDictionary, SequenceDictionary }
+import org.bdgenomics.adam.models.{ SequenceRecord, Attribute, RecordGroupDictionary, SequenceDictionary }
 import org.bdgenomics.adam.util.AttributeUtils
 
 class SAMRecordConverter extends Serializable {
@@ -35,11 +35,7 @@ class SAMRecordConverter extends Serializable {
     // This prevents looking up a -1 in the sequence dictionary
     val readReference: Int = samRecord.getReferenceIndex
     if (readReference != SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
-      builder
-        .setContig(ADAMContig.newBuilder
-          .setContigName(samRecord.getReferenceName)
-          .setContigLength(dict(samRecord.getReferenceName).length)
-          .setReferenceURL(dict(samRecord.getReferenceName).url).build)
+      builder.setContig(SequenceRecord.toADAMContig(dict(samRecord.getReferenceName).get))
 
       val start: Int = samRecord.getAlignmentStart
       if (start != 0) {
@@ -57,11 +53,7 @@ class SAMRecordConverter extends Serializable {
     val mateReference: Int = samRecord.getMateReferenceIndex
 
     if (mateReference != SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
-      builder
-        .setMateContig(ADAMContig.newBuilder
-          .setContigName(samRecord.getMateReferenceName)
-          .setContigLength(dict(samRecord.getMateReferenceName).length)
-          .setReferenceURL(dict(samRecord.getMateReferenceName).url).build)
+      builder.setMateContig(SequenceRecord.toADAMContig(dict(samRecord.getMateReferenceName).get))
 
       val mateStart = samRecord.getMateAlignmentStart
       if (mateStart > 0) {

adam-core/src/main/scala/org/bdgenomics/adam/converters/VariantAnnotationConverter.scala

@@ -16,7 +16,7 @@
 
 package org.bdgenomics.adam.converters
 
-import org.broadinstitute.variant.variantcontext.VariantContext
+import org.broadinstitute.variant.variantcontext.{ Genotype, VariantContext }
 import org.apache.avro.Schema
 import org.apache.avro.specific.SpecificRecord
 import org.broadinstitute.variant.vcf._
@@ -101,7 +101,9 @@ object VariantAnnotationConverter extends Serializable {
     AttrKey("gtFilters", VCFStandardHeaderLines.getFormatLine(VCFConstants.GENOTYPE_FILTER_KEY)),
     AttrKey("genotypeLikelihoods", VCFStandardHeaderLines.getFormatLine(VCFConstants.GENOTYPE_PL_KEY)),
     AttrKey("phaseQuality", attrAsInt _, new VCFFormatHeaderLine(VCFConstants.PHASE_QUALITY_KEY, 1, VCFHeaderLineType.Float, "Read-backed phasing quality")),
-    AttrKey("phaseSetId", attrAsInt _, new VCFFormatHeaderLine(VCFConstants.PHASE_SET_KEY, 1, VCFHeaderLineType.Integer, "Phase set")))
+    AttrKey("phaseSetId", attrAsInt _, new VCFFormatHeaderLine(VCFConstants.PHASE_SET_KEY, 1, VCFHeaderLineType.Integer, "Phase set")),
+    AttrKey("minReadDepth", attrAsInt _, new VCFFormatHeaderLine("MIN_DP", 1, VCFHeaderLineType.Integer, "Minimum DP observed within the GVCF block")),
+    AttrKey("strandBiasComponents", attrAsInt _, new VCFFormatHeaderLine("SB", 4, VCFHeaderLineType.Integer, "Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.")))
 
   lazy val infoHeaderLines: Seq[VCFCompoundHeaderLine] = INFO_KEYS.map(_.hdrLine)
   lazy val formatHeaderLines: Seq[VCFCompoundHeaderLine] = FORMAT_KEYS.map(_.hdrLine)
@@ -143,6 +145,17 @@ object VariantAnnotationConverter extends Serializable {
     fillRecord(VCF2VariantCallingAnnotations, vc, call)
   }
 
+  def convert(g: Genotype, genotype: ADAMGenotype): ADAMGenotype = {
+    for ((v, a) <- VariantAnnotationConverter.VCF2GenotypeAnnotations) {
+      // Add extended attributes if present
+      val attr = g.getExtendedAttribute(v)
+      if (attr != null && attr != VCFConstants.MISSING_VALUE_v4) {
+        genotype.put(a._1, a._2(attr))
+      }
+    }
+    genotype
+  }
+
   def mergeAnnotations(leftRecord: ADAMDatabaseVariantAnnotation, rightRecord: ADAMDatabaseVariantAnnotation): ADAMDatabaseVariantAnnotation = {
     val mergedAnnotation = ADAMDatabaseVariantAnnotation.newBuilder(leftRecord).build()
     val numFields = ADAMDatabaseVariantAnnotation.getClassSchema.getFields.size