Corrected scoring of (bwa) split read edge case in which the split al…

…ignments overlap
PapenfussLab · Dec 13, 2016 · e74d307 · e74d307
1 parent 9d83b63
commit e74d307
Show file tree

Hide file tree

Showing 2 changed files with 28 additions and 2 deletions.
diff --git a/src/main/java/au/edu/wehi/idsv/SplitReadEvidence.java b/src/main/java/au/edu/wehi/idsv/SplitReadEvidence.java
@@ -135,8 +135,12 @@ public float getBreakpointQual() {
 		}
 		int softClipLength = getBreakendSequence().length;
 		if (getSAMRecord().getSupplementaryAlignmentFlag()) {
-			// TODO: better model that can handle multiple realignments and realignment CIGARs
-			softClipLength = getAnchorSequence().length + getUntemplatedSequence().length();
+			ChimericAlignment caThis = new ChimericAlignment(getSAMRecord());
+			// The first record should be the primary
+			ChimericAlignment caPrimary = ChimericAlignment.getChimericAlignments(getSAMRecord()).get(0);
+			// before 
+			BreakendDirection primaryDirectionTowardThis = caThis.getFirstAlignedBaseReadOffset() < caPrimary.getFirstAlignedBaseReadOffset() ^ caPrimary.isNegativeStrand ? BreakendDirection.Backward : BreakendDirection.Forward;
+			softClipLength = SAMRecordUtil.getSoftClipLength(caPrimary.cigar.getCigarElements(), primaryDirectionTowardThis);
 		}
 		return (float)getEvidenceSource().getContext().getConfig().getScoring().getModel().scoreSplitRead(getEvidenceSource().getMetrics(),
 				softClipLength,

diff --git a/src/test/java/au/edu/wehi/idsv/SplitReadEvidenceTest.java b/src/test/java/au/edu/wehi/idsv/SplitReadEvidenceTest.java
@@ -17,6 +17,7 @@
 
 import au.edu.wehi.idsv.picard.InMemoryReferenceSequenceFile;
 import au.edu.wehi.idsv.picard.SynchronousReferenceLookupAdapter;
+import au.edu.wehi.idsv.sam.ChimericAlignment;
 import au.edu.wehi.idsv.sam.SAMRecordUtil;
 import htsjdk.samtools.SAMFileHeader;
 import htsjdk.samtools.SAMRecord;
@@ -259,4 +260,25 @@ public void indel_mismapping_false_positive_assembly_should_not_throw_homology_e
 		}
 		folder.delete();
 	}
+	@Test
+	public void score_should_be_symmetrical_even_if_anchors_overlap() {
+		SAMRecord primary = Read(0, 100, "6M4S");
+		primary.setReadBases(B("NNNNNNNNNN"));
+		primary.setBaseQualities(B("1234567890"));
+		primary.setMappingQuality(40);
+
+		SAMRecord supp = Read(0, 200, "3S7M");
+		supp.setSupplementaryAlignmentFlag(true);
+		supp.setReadBases(B("NNNNNNNNNN"));
+		supp.setBaseQualities(B("1234567890"));
+		supp.setMappingQuality(20);
+
+		primary.setAttribute("SA", new ChimericAlignment(supp).toString());
+		supp.setAttribute("SA", new ChimericAlignment(primary).toString());
+
+		StubSAMEvidenceSource ses = new StubSAMEvidenceSource(getContext(), null, 0, 0, 100);
+		SplitReadEvidence ep = SplitReadEvidence.create(ses, primary).get(0);
+		SplitReadEvidence es = SplitReadEvidence.create(ses, supp).get(0);
+		assertEquals(ep.getBreakpointQual(), es.getBreakpointQual(), 0);
+	}
 }