Skip to content

Commit

Permalink
#213 Excluding local assembly anchoring from SC calculation; added ASC
Browse files Browse the repository at this point in the history
  • Loading branch information
d-cameron committed May 13, 2019
1 parent 144c48f commit 6ce4b34
Show file tree
Hide file tree
Showing 4 changed files with 31 additions and 4 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -277,10 +277,17 @@ 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()));
//supportingAS.stream().forEach(e -> processAnchor(allAnchoredBases, e.getSAMRecord()));
supportingRAS.stream().forEach(e -> processAnchor(allAnchoredBases, e.getSAMRecord()));
supportingCAS.stream().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()));
attribute(VcfInfoAttributes.ASSEMBLY_SUPPORT_CIGAR, makeCigar(assemblyAnchoredBases, nominalPosition).toString());

}

private void updateAssemblySupportTracking() {
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -170,8 +170,8 @@ public VariantContextDirectedEvidence applyConfidenceFilter(ProcessingContext pr
filters.add(VcfFilter.NO_ASSEMBLY.filter());
} else if (v.getBreakpointEvidenceCountLocalAssembly() == 0 || v.getBreakpointEvidenceCountRemoteAssembly() == 0) {
filters.add(VcfFilter.SINGLE_ASSEMBLY.filter());
} else if (v.getBreakpointEvidenceCountLocalAssembly() + v.getBreakpointEvidenceCountRemoteAssembly() > 0 &&
v.getBreakpointEvidenceCountReadPair() + v.getBreakpointEvidenceCountSoftClip() == 0) {
} else if (v.getBreakpointEvidenceCountLocalAssembly() + v.getBreakpointEvidenceCountRemoteAssembly() + v.getBreakpointEvidenceCountCompoundAssembly() > 0 &&
v.getBreakpointEvidenceCountReadPair() + v.getBreakpointEvidenceCountSoftClip() + v.getBreakpointEvidenceCountIndel() == 0) {
filters.add(VcfFilter.ASSEMBLY_ONLY.filter());
}
if (variant.getPhredScaledQual() < processContext.getVariantCallingParameters().breakpointLowQuality) {
Expand Down
3 changes: 2 additions & 1 deletion src/main/java/au/edu/wehi/idsv/vcf/VcfInfoAttributes.java
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,8 @@ public enum VcfInfoAttributes {
SUPPORT_INTERVAL ("SI", VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.Integer, "Support interval offsets from breakend position in which at least one supporting read/read pair/assembly is mapped."),
REMOTE_SUPPORT_INTERVAL ("RSI", VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.Integer, "Support interval offsets of partner breakend."),
INEXACT_HOMPOS ("IHOMPOS", 2, VCFHeaderLineType.Integer, "Position of inexact homology"),
SUPPORT_CIGAR ("SC", 1, VCFHeaderLineType.String, "CIGAR for displaying anchoring alignment of any contributing evidence and microhomologies."),
SUPPORT_CIGAR ("SC", 1, VCFHeaderLineType.String, "CIGAR for displaying anchoring alignment of any contributing evidence and microhomologies. Local assemblies are excluded due to https://github.com/PapenfussLab/gridss/issues/213"),
ASSEMBLY_SUPPORT_CIGAR ("ASC", 1, VCFHeaderLineType.String, "CIGAR encoding assembly contig anchoring alignments. Local assemblies are excluded due to https://github.com/PapenfussLab/gridss/issues/213."),
BREAKEND_ASSEMBLY_ID ("BEID", VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.String, "Identifiers of assemblies supporting the variant."),
BREAKEND_ASSEMBLY_ID_LOCAL_CONTIG_OFFSET ("BEIDL", VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.Integer, "Local chimeric alignment offset of corresponding BEID assembly."),
BREAKEND_ASSEMBLY_ID_REMOTE_CONTIG_OFFSET ("BEIDH", VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.Integer, "Remote chimeric alignment offset of corresponding BEID assembly."),
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -728,6 +728,25 @@ 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() {
ProcessingContext pc = getContext();
AssemblyEvidenceSource aes = AES(pc);
StructuralVariationCallBuilder builder = new StructuralVariationCallBuilder(pc, (VariantContextDirectedEvidence)new IdsvVariantContextBuilder(getContext()) {{
breakpoint(new BreakpointSummary(0, FWD, 100, 1, BWD, 200), "");
phredScore(10);
}}.make());
List<DirectedEvidence> support = Lists.<DirectedEvidence>newArrayList(
SR(Read(1, 200, "3S7M"), Read(0, 98, "3M7S"))
);
SAMRecord ass = AssemblyFactory.createAnchoredBreakend(pc, aes, new SequentialIdGenerator("asm"), BWD, support, fullSupport(support), 1, 200, 7, B("NNNNNNNNNN"), B("NNNNNNNNNN"));
SAMRecord beass = withMapq(44, Read(0, 98, "3M"))[0];
incorporateRealignment(AES(), ass, ImmutableList.of(beass));

builder.addEvidence(SingleReadEvidence.createEvidence(AES(), 1, beass).get(0));
VariantContextDirectedEvidence vc = builder.make();
assertEquals("2M1X", vc.getAttribute(VcfInfoAttributes.ASSEMBLY_SUPPORT_CIGAR.attribute()));
}
@Test
@Ignore("The assembler doesn't know which orientation it was assembling.")
public void spanning_assemblies_should_use_original_parent_assembly_direction_to_determine_local_remote_status() {
IndelEvidence r = IE(withMapq(40, Read(2, 90, "10M1I10M"))[0]);
Expand Down

0 comments on commit 6ce4b34

Please sign in to comment.