Commit 612f9c14 authored by bow's avatar bow
Browse files

Improve mate query to better-disambiguate multiple-mapped reads

parent e96290dc
...@@ -194,16 +194,35 @@ object WipeReads extends ToolCommand { ...@@ -194,16 +194,35 @@ object WipeReads extends ToolCommand {
* @param rec SAMRecord whose mate will be queried * @param rec SAMRecord whose mate will be queried
* @return SAMRecord wrapped in Option * @return SAMRecord wrapped in Option
*/ */
def monadicMateQuery(inBAM: SAMFileReader, rec: SAMRecord): Option[SAMRecord] = def monadicMateQuery(inBAM: SAMFileReader, rec: SAMRecord): Option[SAMRecord] = {
// catching unpaired read here since queryMate will raise an exception if not // adapted from htsjdk's SAMFileReader.queryMate to better deal with multiple-mapped reads
if (!rec.getReadPairedFlag) { if (!rec.getReadPairedFlag)
None return None
} else {
inBAM.queryMate(rec) match { if (rec.getFirstOfPairFlag == rec.getSecondOfPairFlag)
case null => None throw new IllegalArgumentException("SAMRecord must either be first or second of pair, but not both")
case qres @ _ => Some(qres)
val it =
if (rec.getMateReferenceIndex == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX)
inBAM.queryUnmapped()
else
inBAM.queryAlignmentStart(rec.getMateReferenceName, rec.getMateAlignmentStart)
val qres =
try {
it.asScala.toList
.filter(x => x.getReadPairedFlag
&& rec.getAlignmentStart == x.getMateAlignmentStart
&& rec.getMateAlignmentStart == x.getAlignmentStart)
} finally {
it.close()
} }
}
if (qres.length == 1)
Some(qres.head)
else
None
}
/** function to make IntervalTree from our RawInterval objects /** function to make IntervalTree from our RawInterval objects
* *
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment