use of com.hartwig.hmftools.breakpointinspector.clipping.ClipStats in project hmftools by hartwigmedical.
the class Filter method getFilters.
static Collection<String> getFilters(final HMFVariantContext ctx, final SampleStats tumorStats, final SampleStats refStats, final Pair<Location, Location> breakpoints, final float contamination) {
final int MIN_ANCHOR_LENGTH = 30;
final List<Filters> filters = Lists.newArrayList();
if (Stream.of(tumorStats.BP1_Stats, tumorStats.BP2_Stats).mapToInt(s -> s.PR_Only_Normal + s.PR_SR_Normal + s.PR_Only_Support + s.PR_SR_Support).anyMatch(i -> i < 10)) {
filters.add(Filters.MinDepth);
}
final int tumor_SR = Stream.of(tumorStats.BP1_Stats, tumorStats.BP2_Stats).mapToInt(Filter::supportSR).sum();
if (ctx.isInsert()) {
// no PR/SR checks
} else if (ctx.isShortVariant()) {
// short variant logic
final boolean bothSidesHaveSR = Stream.of(tumorStats.BP1_Stats, tumorStats.BP2_Stats).allMatch(s -> supportSR(s) > 0);
final boolean anchorLengthOkay = tumorStats.SR_Evidence.stream().anyMatch(p -> Stream.of(p.getLeft(), p.getRight()).anyMatch(r -> r.getAlignmentEnd() - r.getAlignmentStart() >= MIN_ANCHOR_LENGTH));
if (!bothSidesHaveSR) {
filters.add(Filters.SRSupportZero);
} else if (!anchorLengthOkay) {
filters.add(Filters.MinAnchorLength);
}
// must not have SR support in normal
final int ref_SR = Stream.of(refStats.BP1_Stats, refStats.BP2_Stats).mapToInt(Filter::supportSR).sum();
final int allowableNormalSupport = (int) (contamination * tumor_SR);
if (ref_SR > allowableNormalSupport) {
filters.add(Filters.SRNormalSupport);
}
} else {
// we only need to check BP1 as BP1 PR+PRSR == BP2 PR+PRSR
final int allowableNormalSupport = (int) (contamination * supportPR(tumorStats.BP1_Stats));
if (supportPR(refStats.BP1_Stats) > allowableNormalSupport) {
filters.add(Filters.PRNormalSupport);
}
final boolean anchorLengthOkay = tumorStats.PR_Evidence.stream().anyMatch(p -> Stream.of(p.getLeft(), p.getRight()).allMatch(r -> r.getAlignmentEnd() - r.getAlignmentStart() >= MIN_ANCHOR_LENGTH));
// only applicable for longer variants
final int tumor_PR = Stream.of(tumorStats.BP1_Stats, tumorStats.BP2_Stats).mapToInt(Filter::supportPR).sum();
if (tumor_PR == 0) {
filters.add(Filters.PRSupportZero);
} else if (!anchorLengthOkay) {
filters.add(Filters.MinAnchorLength);
}
}
// we must adjust from Manta breakpoint convention to our clipping position convention
final List<Location> adjusted_bp = Arrays.asList(breakpoints.getLeft().add(ctx.OrientationBP1), breakpoints.getRight().add(ctx.OrientationBP2));
final Set<String> concordant_reads = Sets.newHashSet();
for (final Location bp : adjusted_bp) {
for (final ClipStats t : tumorStats.Sample_Clipping.getSequencesAt(bp)) {
if (t.LongestClipSequence.length() < 5) {
continue;
}
final String tumorSeq = t.Left ? t.LongestClipSequence.substring(t.LongestClipSequence.length() - 5) : t.LongestClipSequence.substring(0, 5);
for (final ClipStats r : refStats.Sample_Clipping.getSequencesAt(bp)) {
if (t.Left != r.Left) {
continue;
} else if (r.LongestClipSequence.length() < 5) {
continue;
}
if (t.Left) {
if (tumorSeq.equals(r.LongestClipSequence.substring(r.LongestClipSequence.length() - 5))) {
concordant_reads.addAll(r.SupportingReads);
}
} else {
if (tumorSeq.equals(r.LongestClipSequence.substring(0, 5))) {
concordant_reads.addAll(r.SupportingReads);
}
}
}
}
}
if (concordant_reads.size() > (int) (contamination * tumor_SR)) {
filters.add(Filters.ClippingConcordance);
}
final Set<String> merged = Sets.newHashSet(ctx.Filter);
merged.addAll(filters.stream().map(Filters::toString).collect(Collectors.toList()));
return merged;
}
Aggregations