Search in sources :

Example 1 with ClipStats

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;
}
Also used : VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine) Arrays(java.util.Arrays) VCFHeader(htsjdk.variant.vcf.VCFHeader) Collection(java.util.Collection) ClipStats(com.hartwig.hmftools.breakpointinspector.clipping.ClipStats) Set(java.util.Set) Collectors(java.util.stream.Collectors) Sets(com.google.common.collect.Sets) List(java.util.List) Stream(java.util.stream.Stream) Lists(com.google.common.collect.Lists) Pair(org.apache.commons.lang3.tuple.Pair) Collections(java.util.Collections) ClipStats(com.hartwig.hmftools.breakpointinspector.clipping.ClipStats)

Aggregations

Lists (com.google.common.collect.Lists)1 Sets (com.google.common.collect.Sets)1 ClipStats (com.hartwig.hmftools.breakpointinspector.clipping.ClipStats)1 VCFFilterHeaderLine (htsjdk.variant.vcf.VCFFilterHeaderLine)1 VCFHeader (htsjdk.variant.vcf.VCFHeader)1 Arrays (java.util.Arrays)1 Collection (java.util.Collection)1 Collections (java.util.Collections)1 List (java.util.List)1 Set (java.util.Set)1 Collectors (java.util.stream.Collectors)1 Stream (java.util.stream.Stream)1 Pair (org.apache.commons.lang3.tuple.Pair)1