Fastq record converter #1185
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -40,6 +40,96 @@ import scala.collection.JavaConversions._ | |
| */ | ||
| private[adam] class FastqRecordConverter extends Serializable with Logging { | ||
|
|
||
| /** | ||
| * Parse 4 lines at a time | ||
| * @see parseReadPairInFastq | ||
| * * | ||
| */ | ||
| private def parseReadInFastq(input: String, | ||
|
There was a problem hiding this comment. Is the use of |
||
| setFirstOfPair: Boolean = false, | ||
| setSecondOfPair: Boolean = false, | ||
| stringency: ValidationStringency = ValidationStringency.STRICT): (String, String, String) = { | ||
| val lines = input.split('\n') | ||
| require(lines.length == 4, | ||
| s"Input must have 4 lines (${lines.length.toString} found):\n${input}") | ||
|
|
||
| val readName = lines(0).drop(1) | ||
| if (readName.endsWith("/1") && setSecondOfPair) | ||
|
There was a problem hiding this comment. I've seen files in the wild that use There was a problem hiding this comment. I have seen both, as well. I am not aware if there is a specification that lists all possibilities. I am thinking of using regex to account for all of them gradually. There was a problem hiding this comment. There isn't a specification, only convention. See http://dx.doi.org/10.1093/nar/gkp1137 There was a problem hiding this comment. addressed by regex like |
||
| throw new Exception( | ||
| s"Found read name $readName ending in '/1' despite second-of-pair flag being set" | ||
| ) | ||
| else if (readName.endsWith("/2") && setFirstOfPair) | ||
| throw new Exception( | ||
| s"Found read name $readName ending in '/2' despite first-of-pair flag being set" | ||
| ) | ||
| val suffix = """(\/1$)|(\/2$)""".r | ||
| val readNameNoSuffix = suffix.replaceAllIn(readName, "") | ||
|
|
||
| val readSequence = lines(1) | ||
| val readQualitiesRaw = lines(3) | ||
|
|
||
| val readQualities = | ||
| if (stringency == ValidationStringency.STRICT) readQualitiesRaw | ||
| else { | ||
| if (readQualitiesRaw == "*") "B" * readSequence.length | ||
| else if (readQualitiesRaw.length < readSequence.length) readQualitiesRaw + ("B" * (readSequence.length - readQualitiesRaw.length)) | ||
| else if (readQualitiesRaw.length > readSequence.length) throw new NotImplementedError("Not implemented") | ||
|
There was a problem hiding this comment.
There was a problem hiding this comment. Also, what's the reason for padding There was a problem hiding this comment. IIRC, There was a problem hiding this comment. https://en.wikipedia.org/wiki/FASTQ_format
There was a problem hiding this comment. fixed |
||
| else readQualitiesRaw | ||
| } | ||
|
|
||
| if (stringency == ValidationStringency.STRICT) { | ||
| if (readQualitiesRaw == "*" && readSequence.length > 1) | ||
| throw new IllegalArgumentException(s"Fastq quality must be defined for\n $input") | ||
| } | ||
|
|
||
| require( | ||
| readSequence.length == readQualities.length, | ||
| s"The first read: ${readName}, has different sequence and qual length." | ||
| ) | ||
|
|
||
| (readNameNoSuffix, readSequence, readQualities) | ||
| } | ||
|
|
||
| private def parseReadPairInFastq(input: String): (String, String, String, String, String, String) = { | ||
| val lines = input.toString.split('\n') | ||
| require(lines.length == 8, | ||
|
There was a problem hiding this comment. I've mentioned this before; perhaps now is the time to fix it? FASTQ format allows for hard line wrapping, so there may be new line characters at any place in the record. There was a problem hiding this comment. Can you give an example for hard line wrapping? What does it look like? There was a problem hiding this comment. There are all kinds of correct and incorrect examples here https://github.com/biojava/biojava/tree/master/biojava-sequencing/src/test/resources/org/biojava/nbio/sequencing/io/fastq There was a problem hiding this comment. I think we should make this work for the simple case first (what we have currently implemented --> fastq record is 4 lines, interleaved read pair is 8 lines). In a follow on, we can make the arbitrary wrapping case work. In my experience, "simply" formatted files are much more common than arbitrarily formatted files. There was a problem hiding this comment. I tried to implement parsing for wrapped lines. Then I found that it would require exact match of sequence length and quality length. Otherwise, it's ambiguous to tell when the quality lines stop. This makes padding with |
||
| s"Record must have 8 lines (${lines.length.toString} found):\n${input}") | ||
|
|
||
| val (firstReadName, firstReadSequence, firstReadQualities) = | ||
| this.parseReadInFastq(lines.take(4).mkString("\n"), setFirstOfPair = true, setSecondOfPair = false) | ||
|
|
||
| val (secondReadName, secondReadSequence, secondReadQualities) = | ||
| this.parseReadInFastq(lines.drop(4).mkString("\n"), setFirstOfPair = false, setSecondOfPair = true) | ||
|
|
||
| ( | ||
| firstReadName, | ||
| firstReadSequence, | ||
| firstReadQualities, | ||
| secondReadName, | ||
| secondReadSequence, | ||
| secondReadQualities | ||
| ) | ||
| } | ||
|
|
||
| private def makeAlignmentRecord(readName: String, | ||
| sequence: String, | ||
| qual: String, | ||
| readInFragment: Int, | ||
| readPaired: Boolean = true, | ||
| recordGroupOpt: Option[String] = None): AlignmentRecord = { | ||
| val builder = AlignmentRecord.newBuilder | ||
| .setReadName(readName) | ||
| .setSequence(sequence) | ||
| .setQual(qual) | ||
| .setReadPaired(readPaired) | ||
| .setReadInFragment(readInFragment) | ||
|
|
||
| if (recordGroupOpt != None) | ||
|
There was a problem hiding this comment. With |
||
| recordGroupOpt.foreach(builder.setRecordGroupName) | ||
|
|
||
| builder.build | ||
| } | ||
|
|
||
| /** | ||
| * Converts a read pair in FASTQ format into two AlignmentRecords. | ||
| * | ||
|
|
@@ -59,57 +149,19 @@ private[adam] class FastqRecordConverter extends Serializable with Logging { | |
| * @see convertFragment | ||
| */ | ||
| def convertPair(element: (Void, Text)): Iterable[AlignmentRecord] = { | ||
| val ( | ||
| firstReadName, | ||
| firstReadSequence, | ||
| firstReadQualities, | ||
| secondReadName, | ||
| secondReadSequence, | ||
| secondReadQualities | ||
| ) = this.parseReadPairInFastq(element._2.toString) | ||
|
There was a problem hiding this comment. You don't need the |
||
|
|
||
| // build and return iterators | ||
| Iterable( | ||
| this.makeAlignmentRecord(firstReadName, firstReadSequence, firstReadQualities, 0), | ||
|
|
||
| this.makeAlignmentRecord(secondReadName, secondReadSequence, secondReadQualities, 1) | ||
| ) | ||
| } | ||
|
|
||
|
|
@@ -126,28 +178,15 @@ private[adam] class FastqRecordConverter extends Serializable with Logging { | |
| * @see convertPair | ||
| */ | ||
| def convertFragment(element: (Void, Text)): Fragment = { | ||
| val ( | ||
| firstReadName, | ||
| firstReadSequence, | ||
| firstReadQualities, | ||
| secondReadName, | ||
| secondReadSequence, | ||
| secondReadQualities | ||
| ) = this.parseReadPairInFastq(element._2.toString) | ||
|
|
||
|
|
||
| require( | ||
| firstReadName == secondReadName, | ||
| "Reads %s and %s in Fragment have different names.".format( | ||
|
|
@@ -156,17 +195,16 @@ private[adam] class FastqRecordConverter extends Serializable with Logging { | |
| ) | ||
| ) | ||
|
|
||
| val alignments = List( | ||
| this.makeAlignmentRecord(firstReadName, firstReadSequence, firstReadQualities, 0), | ||
|
|
||
| this.makeAlignmentRecord(secondReadName, secondReadSequence, secondReadQualities, 1) | ||
| ) | ||
|
|
||
| // build and return record | ||
| Fragment.newBuilder | ||
| .setReadName(firstReadName) | ||
| .setAlignments(alignments) | ||
| .build | ||
| } | ||
|
|
||
| /** | ||
|
|
@@ -193,77 +231,24 @@ private[adam] class FastqRecordConverter extends Serializable with Logging { | |
| setFirstOfPair: Boolean = false, | ||
| setSecondOfPair: Boolean = false, | ||
| stringency: ValidationStringency = ValidationStringency.STRICT): AlignmentRecord = { | ||
|
|
||
| val (readName, readSequence, readQualities) = | ||
| this.parseReadInFastq(element._2.toString, setFirstOfPair, setSecondOfPair, stringency) | ||
|
|
||
|
|
||
| // default to 0 | ||
| val readInFragment = | ||
| if (setSecondOfPair) 1 | ||
| else 0 | ||
|
|
||
| val readPaired = setFirstOfPair || setSecondOfPair | ||
|
|
||
| require( | ||
| (setFirstOfPair && setSecondOfPair) == false, | ||
| "setFirstOfPair and setSecondOfPair cannot be true at the same time" | ||
| ) | ||
|
|
||
| this.makeAlignmentRecord( | ||
|
There was a problem hiding this comment. I have removed all There was a problem hiding this comment.
|
||
| readName, readSequence, readQualities, | ||
| readInFragment, readPaired, recordGroupOpt) | ||
| } | ||
| } | ||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This doc comment doesn't match the method.
Perhaps something like
Return true if the read name suffix and flags match.