Update to htsjdk use in GATK 3.3

biopet-framework/src/main/scala/nl/lumc/sasc/biopet/tools/ExtractAlignedFastq.scala

@@ -9,19 +9,20 @@ import java.io.File
 import scala.collection.mutable.{ Set => MSet }
 import scala.collection.JavaConverters._
 
-import htsjdk.samtools.SAMFileReader
 import htsjdk.samtools.QueryInterval
+import htsjdk.samtools.SamReaderFactory
+import htsjdk.samtools.ValidationStringency
 import htsjdk.samtools.fastq.{ BasicFastqWriter, FastqReader, FastqRecord }
-import htsjdk.tribble.Feature
-import htsjdk.tribble.BasicFeature
+import htsjdk.samtools.util.Interval
 
 import nl.lumc.sasc.biopet.core.ToolCommand
 
 object ExtractAlignedFastq extends ToolCommand {
 
   type FastqPair = (FastqRecord, FastqRecord)
+
   /**
-   * Function to create iterator over features given input interval string
+   * Function to create iterator over Interval given input interval string
    *
    * Valid interval strings are either of these:
    *    - chr5:10000-11000
@@ -38,7 +39,7 @@ object ExtractAlignedFastq extends ToolCommand {
    *
    * @param inStrings iterable yielding input interval string
    */
-  def makeFeatureFromString(inStrings: Iterable[String]): Iterator[Feature] = {
+  def makeIntervalFromString(inStrings: Iterable[String]): Iterator[Interval] = {
 
     // FIXME: can we combine these two patterns into one regex?
     // matches intervals with start and end coordinates
@@ -47,15 +48,16 @@ object ExtractAlignedFastq extends ToolCommand {
     val ptn2 = """([\w_-]+):([\d.,]+)""".r
     // make ints from coordinate strings
     // NOTE: while it is possible for coordinates to exceed Int.MaxValue, we are limited
-    // by the BasicFeature constructor only accepting ints
+    // by the Interval constructor only accepting ints
     def intFromCoord(s: String): Int = s.replaceAll(",", "").replaceAll("\\.", "").toInt
 
     inStrings.map(x => x match {
-        case ptn1(chr, start, end)  => new BasicFeature(chr, intFromCoord(start), intFromCoord(end))
-        case ptn2(chr, start)       => val startCoord = intFromCoord(start)
-                                       new BasicFeature(chr, startCoord, startCoord)
-        case _                      => throw new IllegalArgumentException("Invalid interval string: " + x)
-      })
+      case ptn1(chr, start, end) => new Interval(chr, intFromCoord(start), intFromCoord(end))
+      case ptn2(chr, start) =>
+        val startCoord = intFromCoord(start)
+        new Interval(chr, startCoord, startCoord)
+      case _ => throw new IllegalArgumentException("Invalid interval string: " + x)
+    })
       .toIterator
   }
 
@@ -69,34 +71,30 @@ object ExtractAlignedFastq extends ToolCommand {
    * @param commonSuffixLength length of suffix common to all read pairs
    * @return
    */
-  def makeMembershipFunction(iv: Iterator[Feature],
+  def makeMembershipFunction(iv: Iterator[Interval],
                              inAln: File,
                              minMapQ: Int = 0,
-                             commonSuffixLength: Int = 0
-                            ): (FastqPair => Boolean) = {
-
-    /** function to make interval queries for BAM files */
-    def makeQueryInterval(aln: SAMFileReader, feat: Feature): QueryInterval =
-      if (aln.getFileHeader.getSequenceIndex(feat.getChr) > -1)
-        aln.makeQueryInterval(feat.getChr, feat.getStart, feat.getEnd)
-      else if (feat.getChr.startsWith("chr")
-        && aln.getFileHeader.getSequenceIndex(feat.getChr.substring(3)) > -1)
-        aln.makeQueryInterval(feat.getChr.substring(3), feat.getStart, feat.getEnd)
-      else if (!feat.getChr.startsWith("chr")
-        && aln.getFileHeader.getSequenceIndex("chr" + feat.getChr) > -1)
-        aln.makeQueryInterval("chr" + feat.getChr, feat.getStart, feat.getEnd)
-      else
-        throw new IllegalStateException("Unexpected feature: " + feat.toString)
+                             commonSuffixLength: Int = 0): (FastqPair => Boolean) = {
 
-    val inAlnReader = new SAMFileReader(inAln)
+    val inAlnReader = SamReaderFactory
+      .make()
+      .validationStringency(ValidationStringency.LENIENT)
+      .open(inAln)
     require(inAlnReader.hasIndex)
 
+    def getSequenceIndex(name: String): Int = inAlnReader.getFileHeader.getSequenceIndex(name) match {
+      case x if x >= 0  =>
+        x
+      case otherwise    =>
+        throw new IllegalArgumentException("Chromosome " + name + " is not found in the alignment file")
+    }
+
     val queries: Array[QueryInterval] = iv.toList
-      // sort features
-      .sortBy(x => (x.getChr, x.getStart, x.getEnd))
-      // turn them into QueryInterval objects
-      .map(x => makeQueryInterval(inAlnReader, x))
-      // return as an array
+      // sort Interval
+      .sortBy(x => (x.getSequence, x.getStart, x.getEnd))
+      // transform to QueryInterval
+      .map(x => new QueryInterval(getSequenceIndex(x.getSequence), x.getStart, x.getEnd))
+      // cast to array
       .toArray
 
     lazy val selected: MSet[String] = inAlnReader
@@ -112,11 +110,11 @@ object ExtractAlignedFastq extends ToolCommand {
           logger.debug("Adding " + x.getReadName + " to set ...")
           acc += x.getReadName
         }
-       )
+      )
 
     (pair: FastqPair) => pair._2 match {
-      case null       => selected.contains(pair._1.getReadHeader)
-      case otherwise  =>
+      case null => selected.contains(pair._1.getReadHeader)
+      case otherwise =>
         require(commonSuffixLength < pair._1.getReadHeader.length)
         require(commonSuffixLength < pair._2.getReadHeader.length)
         selected.contains(pair._1.getReadHeader.dropRight(commonSuffixLength))
@@ -131,8 +129,8 @@ object ExtractAlignedFastq extends ToolCommand {
 
     val i1 = inputFastq1.iterator.asScala
     val i2 = inputFastq2 match {
-      case null       => Iterator.continually(null)
-      case otherwise  => otherwise.iterator.asScala
+      case null      => Iterator.continually(null)
+      case otherwise => otherwise.iterator.asScala
     }
     val o1 = outputFastq1
     val o2 = (inputFastq2, outputFastq2) match {
@@ -152,18 +150,18 @@ object ExtractAlignedFastq extends ToolCommand {
         case (rec1, rec2) =>
           o1.write(rec1)
           o2.write(rec2)
-    }
+      }
 
   }
 
-  case class Args (inputBam: File = null,
-                   intervals: List[String] = List.empty[String],
-                   inputFastq1: File = null,
-                   inputFastq2: File = null,
-                   outputFastq1: File = null,
-                   outputFastq2: File = null,
-                   minMapQ: Int = 0,
-                   commonSuffixLength: Int = 0) extends AbstractArgs
+  case class Args(inputBam: File = null,
+                  intervals: List[String] = List.empty[String],
+                  inputFastq1: File = null,
+                  inputFastq2: File = null,
+                  outputFastq1: File = null,
+                  outputFastq2: File = null,
+                  minMapQ: Int = 0,
+                  commonSuffixLength: Int = 0) extends AbstractArgs
 
   class OptParser extends AbstractOptParser {
 
@@ -172,36 +170,44 @@ object ExtractAlignedFastq extends ToolCommand {
         |$commandName - Select aligned FASTQ records
       """.stripMargin)
 
-    opt[File]('I', "input_file") required() valueName "<bam>" action { (x, c) =>
-      c.copy(inputBam = x) } validate {
+    opt[File]('I', "input_file") required () valueName "<bam>" action { (x, c) =>
+      c.copy(inputBam = x)
+    } validate {
       x => if (x.exists) success else failure("Input BAM file not found")
     } text "Input BAM file"
 
-    opt[String]('r', "interval") required() unbounded() valueName "<interval>" action { (x, c) =>
+    opt[String]('r', "interval") required () unbounded () valueName "<interval>" action { (x, c) =>
       // yes, we are appending and yes it's O(n) ~ preserving order is more important than speed here
-      c.copy(intervals = c.intervals :+ x) } text "Interval strings"
+      c.copy(intervals = c.intervals :+ x)
+    } text "Interval strings"
 
-    opt[File]('i', "in1") required() valueName "<fastq>" action { (x, c) =>
-      c.copy(inputFastq1 = x) } validate {
+    opt[File]('i', "in1") required () valueName "<fastq>" action { (x, c) =>
+      c.copy(inputFastq1 = x)
+    } validate {
       x => if (x.exists) success else failure("Input FASTQ file 1 not found")
     } text "Input FASTQ file 1"
 
-    opt[File]('j', "in2") optional() valueName "<fastq>" action { (x, c) =>
-      c.copy(inputFastq1 = x) } validate {
+    opt[File]('j', "in2") optional () valueName "<fastq>" action { (x, c) =>
+      c.copy(inputFastq1 = x)
+    } validate {
       x => if (x.exists) success else failure("Input FASTQ file 2 not found")
     } text "Input FASTQ file 2 (default: none)"
 
-    opt[File]('o', "out1") required() valueName "<fastq>" action { (x, c) =>
-      c.copy(outputFastq1 = x) } text "Output FASTQ file 1"
+    opt[File]('o', "out1") required () valueName "<fastq>" action { (x, c) =>
+      c.copy(outputFastq1 = x)
+    } text "Output FASTQ file 1"
 
-    opt[File]('p', "out2") optional() valueName "<fastq>" action { (x, c) =>
-      c.copy(outputFastq1 = x) } text "Output FASTQ file 2 (default: none)"
+    opt[File]('p', "out2") optional () valueName "<fastq>" action { (x, c) =>
+      c.copy(outputFastq1 = x)
+    } text "Output FASTQ file 2 (default: none)"
 
-    opt[Int]('Q', "min_mapq") optional() action { (x, c) =>
-      c.copy(minMapQ = x) } text "Minimum MAPQ of reads in target region to remove (default: 0)"
+    opt[Int]('Q', "min_mapq") optional () action { (x, c) =>
+      c.copy(minMapQ = x)
+    } text "Minimum MAPQ of reads in target region to remove (default: 0)"
 
-    opt[Int]('s', "read_suffix_length") optional() action { (x, c) =>
-      c.copy(commonSuffixLength = x) } text "Length of common suffix from each read pair (default: 0)"
+    opt[Int]('s', "read_suffix_length") optional () action { (x, c) =>
+      c.copy(commonSuffixLength = x)
+    } text "Length of common suffix from each read pair (default: 0)"
 
     note(
       """
@@ -229,7 +235,7 @@ object ExtractAlignedFastq extends ToolCommand {
       .getOrElse(sys.exit(1))
 
     val memFunc = makeMembershipFunction(
-      iv = makeFeatureFromString(commandArgs.intervals),
+      iv = makeIntervalFromString(commandArgs.intervals),
       inAln = commandArgs.inputBam,
       minMapQ = commandArgs.minMapQ,
       commonSuffixLength = commandArgs.commonSuffixLength)

biopet-framework/src/test/scala/nl/lumc/sasc/biopet/tools/ExtractAlignedFastqUnitTest.scala

@@ -14,8 +14,8 @@ import org.scalatest.mock.MockitoSugar
 import org.scalatest.testng.TestNGSuite
 import org.testng.annotations.{ DataProvider, Test }
 
+import htsjdk.samtools.util.Interval
 import htsjdk.samtools.fastq.{ BasicFastqWriter, FastqReader, FastqRecord }
-import htsjdk.tribble._
 
 class ExtractAlignedFastqUnitTest extends TestNGSuite with MockitoSugar with Matchers {
 
@@ -24,8 +24,8 @@ class ExtractAlignedFastqUnitTest extends TestNGSuite with MockitoSugar with Mat
   private def resourceFile(p: String): File =
     new File(Paths.get(getClass.getResource(p).toURI).toString)
 
-  private def makeFeature(chr: String, start: Int, end: Int): Feature =
-    new BasicFeature(chr, start ,end)
+  private def makeInterval(chr: String, start: Int, end: Int): Interval =
+    new Interval(chr, start, end)
 
   private def makeRecord(header: String): FastqRecord =
     new FastqRecord(header, "ATGC", "", "HIHI")
@@ -40,41 +40,52 @@ class ExtractAlignedFastqUnitTest extends TestNGSuite with MockitoSugar with Mat
     tName + " on " + f.getName + ", read " + rName + ": "
 
   @Test def testIntervalStartEnd() = {
-    val obs = makeFeatureFromString(List("chr5:1000-1100")).next()
-    val exp = new BasicFeature("chr5", 1000, 1100)
-    obs.getChr should === (exp.getChr)
+    val obs = makeIntervalFromString(List("chr5:1000-1100")).next()
+    val exp = new Interval("chr5", 1000, 1100)
+    obs.getSequence should === (exp.getSequence)
     obs.getStart should === (exp.getStart)
     obs.getEnd should === (exp.getEnd)
   }
 
   @Test def testIntervalStartEndComma() = {
-    val obs = makeFeatureFromString(List("chr5:1,000-1,100")).next()
-    val exp = new BasicFeature("chr5", 1000, 1100)
-    obs.getChr should === (exp.getChr)
+    val obs = makeIntervalFromString(List("chr5:1,000-1,100")).next()
+    val exp = new Interval("chr5", 1000, 1100)
+    obs.getSequence should === (exp.getSequence)
     obs.getStart should === (exp.getStart)
     obs.getEnd should === (exp.getEnd)
   }
 
   @Test def testIntervalStartEndDot() = {
-    val obs = makeFeatureFromString(List("chr5:1.000-1.100")).next()
-    val exp = new BasicFeature("chr5", 1000, 1100)
-    obs.getChr should === (exp.getChr)
+    val obs = makeIntervalFromString(List("chr5:1.000-1.100")).next()
+    val exp = new Interval("chr5", 1000, 1100)
+    obs.getSequence should === (exp.getSequence)
     obs.getStart should === (exp.getStart)
     obs.getEnd should === (exp.getEnd)
   }
 
   @Test def testIntervalStart() = {
-    val obs = makeFeatureFromString(List("chr5:1000")).next()
-    val exp = new BasicFeature("chr5", 1000, 1000)
-    obs.getChr should === (exp.getChr)
+    val obs = makeIntervalFromString(List("chr5:1000")).next()
+    val exp = new Interval("chr5", 1000, 1000)
+    obs.getSequence should === (exp.getSequence)
     obs.getStart should === (exp.getStart)
     obs.getEnd should === (exp.getEnd)
   }
 
-  @Test def testIntervalError() =
-    intercept[IllegalArgumentException] {
-      makeFeatureFromString(List("chr5:1000-")).next()
+  @Test def testIntervalError() = {
+    val thrown = intercept[IllegalArgumentException] {
+      makeIntervalFromString(List("chr5:1000-")).next()
+    }
+    thrown.getMessage should === ("Invalid interval string: chr5:1000-")
+  }
+
+  @Test def testMemFuncIntervalError() = {
+    val iv = Iterator(new Interval("chrP", 1, 1000))
+    val inAln = resourceFile("/single01.bam")
+    val thrown = intercept[IllegalArgumentException] {
+      makeMembershipFunction(iv, inAln)
     }
+    thrown.getMessage should === ("Chromosome chrP is not found in the alignment file")
+  }
 
   @DataProvider(name = "singleAlnProvider1", parallel = true)
   def singleAlnProvider1() = {
@@ -84,20 +95,20 @@ class ExtractAlignedFastqUnitTest extends TestNGSuite with MockitoSugar with Mat
 
     Array(
       Array("adjacent left",
-        makeFeature("chrQ", 30, 49), sBam01, sFastq1, sFastq1Default),
+        makeInterval("chrQ", 30, 49), sBam01, sFastq1, sFastq1Default),
       Array("adjacent right",
-        makeFeature("chrQ", 200, 210), sBam01, sFastq1, sFastq1Default),
+        makeInterval("chrQ", 200, 210), sBam01, sFastq1, sFastq1Default),
       Array("no overlap",
-        makeFeature("chrQ", 220, 230), sBam01, sFastq1, sFastq1Default),
+        makeInterval("chrQ", 220, 230), sBam01, sFastq1, sFastq1Default),
       Array("partial overlap",
-        makeFeature("chrQ", 430, 460), sBam01, sFastq1, sFastq1Default.updated("r04", true)),
+        makeInterval("chrQ", 430, 460), sBam01, sFastq1, sFastq1Default.updated("r04", true)),
       Array("enveloped",
-        makeFeature("chrQ", 693, 698), sBam01, sFastq1, sFastq1Default.updated("r03", true))
+        makeInterval("chrQ", 693, 698), sBam01, sFastq1, sFastq1Default.updated("r03", true))
     )
   }
 
   @Test(dataProvider = "singleAlnProvider1")
-  def testSingleBamDefault(name: String, feat: Feature, inAln: File,
+  def testSingleBamDefault(name: String, feat: Interval, inAln: File,
                            fastqMap: Map[String, FastqPair], resultMap: Map[String, Boolean]) = {
     require(resultMap.keySet == fastqMap.keySet)
     val memFunc = makeMembershipFunction(Iterator(feat), inAln)
@@ -116,16 +127,16 @@ class ExtractAlignedFastqUnitTest extends TestNGSuite with MockitoSugar with Mat
 
     Array(
       Array("less than minimum MAPQ",
-        makeFeature("chrQ", 830, 890), sBam02, 60, sFastq2, sFastq2Default),
+        makeInterval("chrQ", 830, 890), sBam02, 60, sFastq2, sFastq2Default),
       Array("greater than minimum MAPQ",
-        makeFeature("chrQ", 830, 890), sBam02, 20, sFastq2, sFastq2Default.updated("r07", true)),
+        makeInterval("chrQ", 830, 890), sBam02, 20, sFastq2, sFastq2Default.updated("r07", true)),
       Array("equal to minimum MAPQ",
-        makeFeature("chrQ", 260, 320), sBam02, 30, sFastq2, sFastq2Default.updated("r01", true))
+        makeInterval("chrQ", 260, 320), sBam02, 30, sFastq2, sFastq2Default.updated("r01", true))
     )
   }
 
   @Test(dataProvider = "singleAlnProvider2")
-  def testSingleBamMinMapQ(name: String, feat: Feature, inAln: File, minMapQ: Int,
+  def testSingleBamMinMapQ(name: String, feat: Interval, inAln: File, minMapQ: Int,
                            fastqMap: Map[String, FastqPair], resultMap: Map[String, Boolean]) = {
     require(resultMap.keySet == fastqMap.keySet)
     val memFunc = makeMembershipFunction(Iterator(feat), inAln, minMapQ)
@@ -148,22 +159,22 @@ class ExtractAlignedFastqUnitTest extends TestNGSuite with MockitoSugar with Mat
 
     Array(
       Array("adjacent left",
-        makeFeature("chrQ", 30, 49), pBam01, pFastq1, pFastq1Default),
+        makeInterval("chrQ", 30, 49), pBam01, pFastq1, pFastq1Default),
       Array("adjacent right",
-        makeFeature("chrQ", 200, 210), pBam01, pFastq1, pFastq1Default),
+        makeInterval("chrQ", 200, 210), pBam01, pFastq1, pFastq1Default),
       Array("no overlap",
-        makeFeature("chrQ", 220, 230), pBam01, pFastq1, pFastq1Default),
+        makeInterval("chrQ", 220, 230), pBam01, pFastq1, pFastq1Default),
       Array("partial overlap",
-        makeFeature("chrQ", 430, 460), pBam01, pFastq1, pFastq1Default.updated("r04", true)),
+        makeInterval("chrQ", 430, 460), pBam01, pFastq1, pFastq1Default.updated("r04", true)),
       Array("enveloped",
-        makeFeature("chrQ", 693, 698), pBam01, pFastq1, pFastq1Default.updated("r03", true)),
+        makeInterval("chrQ", 693, 698), pBam01, pFastq1, pFastq1Default.updated("r03", true)),
       Array("in intron",
-        makeFeature("chrQ", 900, 999), pBam01, pFastq1, pFastq1Default.updated("r05", true))
+        makeInterval("chrQ", 900, 999), pBam01, pFastq1, pFastq1Default.updated("r05", true))
     )
   }
 
   @Test(dataProvider = "pairAlnProvider1")
-  def testPairBamDefault(name: String, feat: Feature, inAln: File,
+  def testPairBamDefault(name: String, feat: Interval, inAln: File,
                          fastqMap: Map[String, FastqPair], resultMap: Map[String, Boolean]) = {
     require(resultMap.keySet == fastqMap.keySet)
     val memFunc = makeMembershipFunction(Iterator(feat), inAln, commonSuffixLength = 2)
@@ -175,13 +186,13 @@ class ExtractAlignedFastqUnitTest extends TestNGSuite with MockitoSugar with Mat
   }
 
   @Test def testWriteSingleBamDefault() = {
-      val memFunc = (recs: FastqPair) => Set("r01", "r03").contains(recs._1.getReadHeader)
-      val in1 = new FastqReader(resourceFile("/single01.fq"))
-      val mo1 = mock[BasicFastqWriter]
-      selectFastqReads(memFunc, in1, mo1)
-      verify(mo1, times(2)).write(anyObject.asInstanceOf[FastqRecord])
-      verify(mo1).write(new FastqRecord("r01", "A", "", "H"))
-      verify(mo1).write(new FastqRecord("r03", "G", "", "H"))
+    val memFunc = (recs: FastqPair) => Set("r01", "r03").contains(recs._1.getReadHeader)
+    val in1 = new FastqReader(resourceFile("/single01.fq"))
+    val mo1 = mock[BasicFastqWriter]
+    selectFastqReads(memFunc, in1, mo1)
+    verify(mo1, times(2)).write(anyObject.asInstanceOf[FastqRecord])
+    verify(mo1).write(new FastqRecord("r01", "A", "", "H"))
+    verify(mo1).write(new FastqRecord("r03", "G", "", "H"))
   }
 
   @Test def testWritePairBamDefault() = {