Supports access to indexed fa and fasta files

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

@@ -0,0 +1,87 @@
+/**
+ * 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.adam.util
+
+import htsjdk.samtools.ValidationStringency
+import htsjdk.samtools.reference.{ FastaSequenceIndex, IndexedFastaSequenceFile }
+import java.net.URI
+import java.nio.file.Paths
+import org.apache.hadoop.fs.Path
+import org.apache.spark.SparkContext
+import org.bdgenomics.adam.models.{ SequenceDictionary, ReferenceRegion }
+import org.bdgenomics.utils.misc.Logging
+
+/**
+ * Loads and extracts sequences directly from indexed fasta or fa files. filePath requires fai index in the
+ * same directory with same naming convention.
+ *
+ * @param filePath path to fasta or fa index
+ */
+case class IndexedFastaFile(sc: SparkContext,
+                            filePath: String,
+                            stringency: ValidationStringency = ValidationStringency.STRICT)
+    extends ReferenceFile with Logging {
+
+  // Generate IndexedFastaSequenceFile from path and fai index
+  private val ref: IndexedFastaSequenceFile = {
+
+    // get absolute path and scheme to create URI
+    val path = new Path(filePath).toString
+    val scheme = new Path(path).getFileSystem(sc.hadoopConfiguration).getScheme
+    val pathWithScheme = s"${scheme}://${path}"
+
+    val uri = new URI(pathWithScheme)
+
+    val file = Paths.get(uri).toFile
+    val uriIdx = new URI(pathWithScheme + ".fai")
+    val pathIdx = Paths.get(uriIdx)
+
+    val idx = new FastaSequenceIndex(pathIdx.toFile)
+    new IndexedFastaSequenceFile(file, idx)
+  }
+
+  // Get sequence dictionary. If sequence dictionary is not defined,
+  // generate sequence dictionary from file
+  val sequences =
+    try {
+      SequenceDictionary(ref.getSequenceDictionary)
+    } catch {
+      case e: Throwable => {
+        if (stringency == ValidationStringency.STRICT) {
+          throw e
+        } else {
+          if (stringency == ValidationStringency.LENIENT) {
+            log.warn("Caught exception %s when loading FASTA sequence dictionary. Using empty dictionary instead.".format(e))
+          }
+          SequenceDictionary.empty
+        }
+      }
+    }
+
+  /**
+   * Extracts base sequence from FastaSequenceIndex
+   * @param region The desired ReferenceRegion to extract.
+   * @return The reference sequence at the desired locus.
+   */
+  def extract(region: ReferenceRegion): String = {
+    ref.getSubsequenceAt(region.referenceName, region.start, region.end)
+      .getBaseString
+  }
+}
+

adam-core/src/test/resources/HLA_DQB1_05_01_01_02.dict

@@ -0,0 +1,2 @@
+@HD	VN:1.5
+@SQ	SN:HLA-DQB1*05:01:01:02	LN:7090	M5:0f304adf7acf3bd4b7c54c1394c85a4b	UR:file:/Users/akmorrow/ADAM/adam/adam-core/src/test/resources/HLA_DQB1_05_01_01_02.fa

adam-core/src/test/resources/HLA_DQB1_05_01_01_02.fa.fai

@@ -0,0 +1 @@
+HLA-DQB1*05:01:01:02	7090	39	72	73

adam-core/src/test/scala/org/bdgenomics/adam/util/IndexedFastaFileSuite.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.adam.util
+
+import htsjdk.samtools.{ ValidationStringency, SAMException }
+import org.bdgenomics.adam.models.ReferenceRegion
+
+class IndexedFastaFileSuite extends ADAMFunSuite {
+
+  val filePath = resourceUrl("HLA_DQB1_05_01_01_02.fa").getPath
+
+  sparkTest("correctly generates sequence dictionary from .dict file") {
+    val indexedFasta = IndexedFastaFile(sc, filePath)
+    assert(indexedFasta.sequences.records.length == 1)
+  }
+
+  sparkTest("correctly gets sequence") {
+    val indexedFasta = IndexedFastaFile(sc, filePath)
+    val region = ReferenceRegion("HLA-DQB1*05:01:01:02", 1, 50)
+    val sequence = indexedFasta.extract(region)
+    assert(sequence == "TTCTAAGACCTTTGCTCTTCTCCCCAGGACTTAAGGCTCTTCAGCGTGTC")
+  }
+
+  sparkTest("fails when fai index is not provided") {
+    val pathWithoutIndex = resourceUrl("hs38DH_chr1_10.fa").getPath
+
+    try {
+      IndexedFastaFile(sc, pathWithoutIndex)
+      assert(false)
+    } catch {
+      case s: SAMException => assert(true)
+      case e: Throwable    => assert(false)
+    }
+  }
+
+  sparkTest("passes when dict is not provided and ValidationStringency = LENIENT") {
+    val pathWithoutDict = resourceUrl("artificial.fa").getPath
+
+    try {
+      IndexedFastaFile(sc, pathWithoutDict, ValidationStringency.LENIENT)
+      assert(true)
+    } catch {
+      case s: SAMException => assert(false)
+      case e: Exception    => assert(false)
+    }
+  }
+}
+