#213 fixes indel assemblies only being excluded from anchoring interv…

…al on one side of the assembly
PapenfussLab · Sep 9, 2019 · 6f89151 · 6f89151
1 parent 0f789d6
commit 6f89151
Show file tree

Hide file tree

Showing 2 changed files with 55 additions and 10 deletions.
diff --git a/src/main/java/au/edu/wehi/idsv/StructuralVariationCallBuilder.java b/src/main/java/au/edu/wehi/idsv/StructuralVariationCallBuilder.java
@@ -276,17 +276,27 @@ private void updateNominalCallPosition() {
 		supportingSR.stream().flatMap(l -> l.stream()).forEach(e -> processAnchor(allAnchoredBases, e.getSAMRecord()));
 		supportingIndel.stream().flatMap(l -> l.stream()).forEach(e -> processAnchor(allAnchoredBases, e.getSAMRecord()));
 		supportingDP.stream().flatMap(l -> l.stream()).forEach(e -> processAnchor(allAnchoredBases, e.getLocalledMappedRead()));
-		//supportingAS.stream().forEach(e -> processAnchor(allAnchoredBases, e.getSAMRecord()));
-		supportingRAS.stream().forEach(e -> processAnchor(allAnchoredBases, e.getSAMRecord()));
-		supportingCAS.stream().forEach(e -> processAnchor(allAnchoredBases, e.getSAMRecord()));
+		supportingAS.stream().filter(e -> !shouldfilterAssemblyFromSupportInterval(e)).forEach(e -> processAnchor(allAnchoredBases, e.getSAMRecord()));
+		supportingRAS.stream().filter(e -> !shouldfilterAssemblyFromSupportInterval(e)).forEach(e -> processAnchor(allAnchoredBases, e.getSAMRecord()));
+		supportingCAS.stream().filter(e -> !shouldfilterAssemblyFromSupportInterval(e)).forEach(e -> processAnchor(allAnchoredBases, e.getSAMRecord()));
 		attribute(VcfInfoAttributes.SUPPORT_CIGAR, makeCigar(allAnchoredBases, nominalPosition).toString());
 		RangeSet<Integer> assemblyAnchoredBases = TreeRangeSet.create();
-		// TODO: implement https://github.com/PapenfussLab/gridss/issues/213 so we can include local assemblies
-		//supportingAS.stream().forEach(e -> processAnchor(allAnchoredBases, e.getSAMRecord()));
-		supportingRAS.stream().forEach(e -> processAnchor(assemblyAnchoredBases, e.getSAMRecord()));
-		supportingCAS.stream().forEach(e -> processAnchor(assemblyAnchoredBases, e.getSAMRecord()));
+		supportingAS.stream().filter(e -> !shouldfilterAssemblyFromSupportInterval(e)).forEach(e -> processAnchor(allAnchoredBases, e.getSAMRecord()));
+		supportingRAS.stream().filter(e -> !shouldfilterAssemblyFromSupportInterval(e)).forEach(e -> processAnchor(assemblyAnchoredBases, e.getSAMRecord()));
+		supportingCAS.stream().filter(e -> !shouldfilterAssemblyFromSupportInterval(e)).forEach(e -> processAnchor(assemblyAnchoredBases, e.getSAMRecord()));
 		attribute(VcfInfoAttributes.ASSEMBLY_SUPPORT_CIGAR, makeCigar(assemblyAnchoredBases, nominalPosition).toString());
+	}
 
+	/**
+	 * https://github.com/PapenfussLab/gridss/issues/213
+	 * assemblies can assemble fragment that do not actually support the breakpoint
+	 * (e.g. nearby reads noise soft-clip). This means we can't use local assemblies
+	 * in anchor support interval calculation.
+	 * Note: this includes both sides of indel-spanning assemblies as well.
+	 */
+	private static boolean shouldfilterAssemblyFromSupportInterval(SingleReadEvidence e) {
+		// TODO: are both sides of indel assemblies
+		return !e.getSAMRecord().isSecondaryOrSupplementary();
 	}
 
 	private void updateAssemblySupportTracking() {

diff --git a/src/test/java/au/edu/wehi/idsv/StructuralVariationCallBuilderTest.java b/src/test/java/au/edu/wehi/idsv/StructuralVariationCallBuilderTest.java
@@ -728,7 +728,7 @@ public void anchor_cigar_should_use_local_coordinates() {
 		assertEquals("9M1X", vc.getAttribute(VcfInfoAttributes.SUPPORT_CIGAR.attribute()));
 	}
 	@Test
-	public void assembly_anchor_cigar_should_use_only_remote_assemblies() {
+	public void bug214_assembly_anchor_cigar_should_use_only_remote_assemblies() {
 		ProcessingContext pc = getContext();
 		AssemblyEvidenceSource aes = AES(pc);
 		StructuralVariationCallBuilder builder = new StructuralVariationCallBuilder(pc, (VariantContextDirectedEvidence)new IdsvVariantContextBuilder(getContext()) {{
@@ -742,9 +742,44 @@ public void assembly_anchor_cigar_should_use_only_remote_assemblies() {
 		SAMRecord beass = withMapq(44, Read(0, 98, "3M"))[0];
 		incorporateRealignment(AES(), ass, ImmutableList.of(beass));
 
+		// Remote stays
 		builder.addEvidence(SingleReadEvidence.createEvidence(AES(), 1, beass).get(0));
-		VariantContextDirectedEvidence vc = builder.make();
-		assertEquals("2M1X", vc.getAttribute(VcfInfoAttributes.ASSEMBLY_SUPPORT_CIGAR.attribute()));
+		assertEquals("2M1X", builder.make().getAttribute(VcfInfoAttributes.ASSEMBLY_SUPPORT_CIGAR.attribute()));
+
+		// Local gets filtered
+		builder = new StructuralVariationCallBuilder(pc, (VariantContextDirectedEvidence)new IdsvVariantContextBuilder(getContext()) {{
+			breakpoint(new BreakpointSummary(1, BWD, 200, 0, FWD, 100), "");
+			phredScore(10);
+		}}.make());
+		builder.addEvidence(SingleReadEvidence.createEvidence(AES(), 1, ass).get(0));
+		assertEquals("1X", builder.make().getAttribute(VcfInfoAttributes.ASSEMBLY_SUPPORT_CIGAR.attribute()));
+	}
+	@Test
+	public void bug214_assembly_anchor_cigar_should_not_use_indel_assembly() {
+		ProcessingContext pc = getContext();
+		AssemblyEvidenceSource aes = AES(pc);
+		StructuralVariationCallBuilder builder;
+		List<DirectedEvidence> support = Lists.<DirectedEvidence>newArrayList(
+				SR(Read(3, 98, "3M7S"), Read(3, 200, "7M"))
+		);
+		SAMRecord ass = AssemblyFactory.createAnchoredBreakend(pc, aes, new SequentialIdGenerator("asm"), FWD, support, fullSupport(support), 3, 100, 3, B("NNNNNNNNNN"), B("NNNNNNNNNN"));
+		ass.setCigarString("3M100D7M");
+		List<SingleReadEvidence> indelass = SingleReadEvidence.createEvidence(AES(), 1, ass);
+		assertEquals(2, indelass.size()); // both sides of indel assembly
+
+		builder = new StructuralVariationCallBuilder(pc, (VariantContextDirectedEvidence)new IdsvVariantContextBuilder(getContext()) {{
+			breakpoint(new BreakpointSummary(3, FWD, 100, 3, BWD, 200), "");
+			phredScore(10);
+		}}.make());
+		builder.addEvidence(indelass.get(0));
+		assertEquals("1X", builder.make().getAttribute(VcfInfoAttributes.ASSEMBLY_SUPPORT_CIGAR.attribute()));
+
+		builder = new StructuralVariationCallBuilder(pc, (VariantContextDirectedEvidence)new IdsvVariantContextBuilder(getContext()) {{
+			breakpoint(new BreakpointSummary(3, BWD, 200, 3, FWD, 100), "");
+			phredScore(10);
+		}}.make());
+		builder.addEvidence(indelass.get(1));
+		assertEquals("1X", builder.make().getAttribute(VcfInfoAttributes.ASSEMBLY_SUPPORT_CIGAR.attribute()));
 	}
 	@Test
 	@Ignore("The assembler doesn't know which orientation it was assembling.")