Reverse complement negative strand reads in fastq output

adam-cli/src/test/scala/org/bdgenomics/adam/cli/ADAM2FastqSuite.scala

@@ -0,0 +1,82 @@
+/**
+ * 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.cli
+
+import com.google.common.io.Files
+import java.io.File
+import org.apache.spark.rdd.RDD
+import org.bdgenomics.adam.rdd.ADAMContext._
+import org.bdgenomics.adam.util.ADAMFunSuite
+import org.bdgenomics.formats.avro.AlignmentRecord
+
+class Adam2FastqSuite extends ADAMFunSuite {
+
+  sparkTest("convert SAM to paired FASTQ") {
+    val readsFilepath = ClassLoader.getSystemClassLoader.getResource("bqsr1.sam").getFile
+
+    // The following fastq files were generated by Picard's SamToFastq
+
+    // First generate mapped on SAM file (sorted by readName to make the comparison easier)
+
+    // samtools view -H adam-core/src/test/resources/bqsr1.sam > adam-core/src/test/resources/bqsr1-readnamesorted.sam
+    // samtools view -F 4 adam-core/src/test/resources/bqsr1.sam | sort -k1 >> adam-core/src/test/resources/bqsr1-readnamesorted.sam
+
+    // java -jar picard.jar
+    // SamToFastq 
+    // I=adam-core/src/test/resources/bqsr1-readnamesorted.sam 
+    // FASTQ=adam-core/src/test/resources/bqsr1-r1.fq 
+    // SECOND_END_FASTQ=adam-core/src/test/resources/bqsr1-r2.fq 
+    // VALIDATION_STRINGENCY=SILENT 
+
+    // VALIDATION_STRINGENCY=SILENT is necessary since they are unpaired reads and this matches the ADAM default
+
+    val fastq1Path = ClassLoader.getSystemClassLoader.getResource("bqsr1-r1.fq").getFile
+    val fastq2Path = ClassLoader.getSystemClassLoader.getResource("bqsr1-r2.fq").getFile
+
+    val outputDir = Files.createTempDir()
+    val outputFastqR1File = outputDir.getAbsolutePath + "/bqsr1-r1.fq"
+    val outputFastqR2File = outputDir.getAbsolutePath + "/bqsr1-r2.fq"
+
+    // Only looking at mapped reads
+    // This is because Picard and ADAM disagree on setting negative strand on unmapped reads
+    // Picard allows unmapped reads to set the negative strand flag and therefore reverse-complemented on output
+    val reads: RDD[AlignmentRecord] =
+      sc
+        .loadAlignments(readsFilepath)
+        .filter(r => r.getReadMapped != null && r.getReadMapped)
+
+    reads.adamSaveAsFastq(outputFastqR1File, Some(outputFastqR2File), sort = true)
+
+    val goldR1Reads =
+      scala.io.Source.fromFile(new File(fastq1Path)).getLines().toSeq
+
+    val goldR2Reads =
+      scala.io.Source.fromFile(new File(fastq2Path)).getLines().toSeq
+
+    val outputR1Reads = scala.io.Source.fromFile(new File(outputFastqR1File + "/part-00000")).getLines().toSeq
+    val outputR2Reads = scala.io.Source.fromFile(new File(outputFastqR2File + "/part-00000")).getLines().toSeq
+
+    assert(outputR1Reads.length === goldR1Reads.length)
+    assert(outputR2Reads.length === goldR2Reads.length)
+
+    outputR1Reads.zip(goldR1Reads).foreach(kv => assert(kv._1 === kv._2))
+    outputR2Reads.zip(goldR2Reads).foreach(kv => assert(kv._1 === kv._2))
+
+  }
+
+}

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

@@ -18,11 +18,10 @@
 package org.bdgenomics.adam.converters
 
 import htsjdk.samtools.{ SAMFileHeader, SAMRecord }
-import org.bdgenomics.adam.models.{ RecordGroupDictionary, SAMFileHeaderWritable, SequenceDictionary }
-import org.bdgenomics.adam.rdd.ADAMContext._
+import org.bdgenomics.adam.instrumentation.Timers._
+import org.bdgenomics.adam.models._
 import org.bdgenomics.adam.rich.RichAlignmentRecord
 import org.bdgenomics.formats.avro.AlignmentRecord
-import org.bdgenomics.adam.instrumentation.Timers._
 
 class AlignmentRecordConverter extends Serializable {
 
@@ -56,8 +55,14 @@ class AlignmentRecordConverter extends Serializable {
     "@%s%s\n%s\n+\n%s".format(
       adamRecord.getReadName,
       readNameSuffix,
-      adamRecord.getSequence,
-      adamRecord.getQual
+      if (adamRecord.getReadMapped && adamRecord.getReadNegativeStrand)
+        Alphabet.dna.reverseComplement(adamRecord.getSequence)
+      else
+        adamRecord.getSequence,
+      if (adamRecord.getReadMapped && adamRecord.getReadNegativeStrand)
+        adamRecord.getQual.reverse
+      else
+        adamRecord.getQual
     )
   }
 

adam-core/src/test/scala/org/bdgenomics/adam/converters/AlignmentRecordConverterSuite.scala

@@ -17,6 +17,9 @@
  */
 package org.bdgenomics.adam.converters
 
+import java.io.File
+
+import htsjdk.samtools.{ SamReaderFactory, SAMRecord }
 import org.bdgenomics.formats.avro.{ AlignmentRecord, Contig }
 import org.bdgenomics.adam.models.{
   RecordGroupDictionary,
@@ -27,12 +30,14 @@ import org.bdgenomics.adam.models.{
 }
 import org.scalatest.FunSuite
 
+import scala.collection.JavaConversions._
+
 class AlignmentRecordConverterSuite extends FunSuite {
 
   // allocate converters
   val adamRecordConverter = new AlignmentRecordConverter
 
-  def make_read(start: Long, cigar: String, mdtag: String, length: Int, id: Int = 0, nullQuality: Boolean = false): AlignmentRecord = {
+  def makeRead(start: Long, cigar: String, mdtag: String, length: Int, id: Int = 0, nullQuality: Boolean = false): AlignmentRecord = {
     val sequence: String = "A" * length
     val builder = AlignmentRecord.newBuilder()
       .setReadName("read" + id.toString)
@@ -54,7 +59,7 @@ class AlignmentRecordConverterSuite extends FunSuite {
   }
 
   test("testing the fields in a converted ADAM Read") {
-    val adamRead = make_read(3L, "2M3D2M", "2^AAA2", 4)
+    val adamRead = makeRead(3L, "2M3D2M", "2^AAA2", 4)
 
     // add reference details
     adamRead.setRecordGroupName("record_group")
@@ -102,7 +107,7 @@ class AlignmentRecordConverterSuite extends FunSuite {
   }
 
   test("converting a read with null quality is OK") {
-    val adamRead = make_read(3L, "2M3D2M", "2^AAA2", 4, nullQuality = true)
+    val adamRead = makeRead(3L, "2M3D2M", "2^AAA2", 4, nullQuality = true)
 
     // add reference details
     adamRead.setRecordGroupName("record_group")
@@ -165,5 +170,80 @@ class AlignmentRecordConverterSuite extends FunSuite {
     assert(fastq(2) === "+")
     assert(fastq(3) === ".+**.+;:**.")
   }
+
+  def getSAMRecordFromReadName(readName: String): (AlignmentRecord, AlignmentRecord) = {
+    val samToADAMConverter = new SAMRecordConverter
+    val SAMTestFile = new File(getClass.getClassLoader.getResource("bqsr1.sam").getFile)
+    val newSAMReader = SamReaderFactory.makeDefault().open(SAMTestFile)
+
+    // Obtain SAMRecord
+    val newSAMRecord = newSAMReader.iterator().dropWhile(r => r.getReadName != readName)
+    val newSequenceRecord = new SequenceRecord("22", 51304566)
+    val newSequenceDictionary = SequenceDictionary(newSequenceRecord)
+    val firstRecord = samToADAMConverter.convert(newSAMRecord.next(), newSequenceDictionary, new RecordGroupDictionary(Seq()))
+    val secondRecord = samToADAMConverter.convert(newSAMRecord.next(), newSequenceDictionary, new RecordGroupDictionary(Seq()))
+    (firstRecord, secondRecord)
+  }
+
+  test("reverse complement reads when converting to fastq") {
+
+    // SRR062634.10022079      83      22      16082719        0       5S95M   =       16082635        -179    
+    // AAGTAGCTGGGACTACACGCACGCACCACCATGCCTGGCTAATTTTTGTATTTTTAGTAGAGATGAGGTTTCACCATATTGGCCAGGCTGGTTTTGAATT    
+    // #####EB5BB<840&:2?>A?-AC8=,5@AABCB?CEDBDC@6BB,CA0CB,B-DEDEDEDEA:D?DE5EBEC?E?5?D:AEEEDEDDEEE=BEEBDD-?    
+    // RG:Z:SRR062634  XC:i:95 XT:A:R  NM:i:2  SM:i:0  AM:i:0  X0:i:3  X1:i:0  XM:i:2  XO:i:0  XG:i:0  MD:Z:15G0T78    
+    // XA:Z:GL000244.1,+31092,100M,2;14,+19760216,100M,2;
+
+    val (firstRecord, secondRecord) = getSAMRecordFromReadName("SRR062634.10022079")
+
+    val firstRecordFastq = adamRecordConverter.convertToFastq(firstRecord, maybeAddSuffix = true)
+      .toString
+      .split('\n')
+
+    assert(firstRecordFastq(0) === "@SRR062634.10022079/2")
+    assert(firstRecordFastq(1) === "CTGGAGTGCAGTGGCATGATTTCAGCTCACTGTCGTCTCTGCCTCCCTGACTCAAGTGATTCTCCTGCCTCAGCCTCCCACGTCGCTCGGACTCCACGCC")
+    assert(firstRecordFastq(2) === "+")
+    assert(firstRecordFastq(3) === "A:=D5D5E?D?DDD:.@@@@=?EE=DADDB@D=DD??ED=:CCCC?D:E=EEB=-C>C=@=EEEEB5EC-?A>=C-C?DC+34+4A>-?5:=/-A=@>>:")
+
+    val secondRecordFastq = adamRecordConverter.convertToFastq(secondRecord, maybeAddSuffix = true)
+      .toString
+      .split('\n')
+
+    assert(secondRecordFastq(0) === "@SRR062634.10022079/1")
+    assert(secondRecordFastq(1) === "AATTCAAAACCAGCCTGGCCAATATGGTGAAACCTCATCTCTACTAAAAATACAAAAATTAGCCAGGCATGGTGGTGCGTGCGTGTAGTCCCAGCTACTT")
+    assert(secondRecordFastq(2) === "+")
+    assert(secondRecordFastq(3) === "?-DDBEEB=EEEDDEDEEEA:D?5?E?CEBE5ED?D:AEDEDEDED-B,BC0AC,BB6@CDBDEC?BCBAA@5,=8CA-?A>?2:&048<BB5BE#####")
+
+  }
+
+  test("converting to fastq with unmapped reads") {
+    //SRR062634.10448889      117     22      16079761        0       *       =       16079761        0       
+    // TTTCTTTCTTTTATATATATATACACACACACACACACACACACACATATATGTATATATACACGTATATGTATGTATATATGTATATATACACGTATAT    
+    // @DF>C;FDC=EGEGGEFDGEFDD?DFDEEGFGFGGGDGGGGGGGEGGGGFGGGFGGGGGGFGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG    
+    // RG:Z:SRR062634
+
+    val (secondRecord, firstRecord) = getSAMRecordFromReadName("SRR062634.10448889")
+
+    val firstRecordFastq = adamRecordConverter.convertToFastq(firstRecord, maybeAddSuffix = true)
+      .toString
+      .split('\n')
+
+    assert(firstRecord.getReadMapped)
+    assert(firstRecord.getReadNegativeStrand)
+    assert(firstRecordFastq(0) === "@SRR062634.10448889/2")
+    assert(firstRecordFastq(1) === "ACCTGTCTCAGCCTCCCAAAGTGCTGCGATTACAGTCATGAGCCACCGCACTTGGCTGGGTTTTCGTTTTCTTTCTTTTATATATATATACACACACACA")
+    assert(firstRecordFastq(2) === "+")
+    assert(firstRecordFastq(3) === "GGGGGGGGGGGGGGGGGGGGGEGGGGGGGGGGGGGGGGGGGGGGGGGFGEGEEDGGFDF?AEEEBDADEEDEEE;DFC@'B:B=B=B=BADCBCBCA=DA")
+
+    val secondRecordFastq = adamRecordConverter.convertToFastq(secondRecord, maybeAddSuffix = true)
+      .toString
+      .split('\n')
+
+    assert(!secondRecord.getReadMapped)
+    assert(!secondRecord.getReadNegativeStrand)
+    assert(secondRecordFastq(0) === "@SRR062634.10448889/1")
+    assert(secondRecordFastq(1) === secondRecord.getSequence)
+    assert(secondRecordFastq(2) === "+")
+    assert(secondRecordFastq(3) === secondRecord.getQual)
+  }
 }