Examples with Clipping com.hartwig.hmftools.breakpointinspector.clipping.Clipping used on opensource projects

Example 1 with Clipping

use of com.hartwig.hmftools.breakpointinspector.clipping.Clipping in project hmftools by hartwigmedical.

the class Analysis method processStructuralVariant.

StructuralVariantResult processStructuralVariant(final HMFVariantContext ctx) throws IOException {
    final QueryInterval[] intervals = QueryInterval.optimizeIntervals(new QueryInterval[] { new QueryInterval(ctx.MantaBP1.ReferenceIndex, Math.max(0, ctx.MantaBP1.Position + ctx.Uncertainty1.Start - range), ctx.MantaBP1.Position + ctx.Uncertainty1.End + range), new QueryInterval(ctx.MantaBP2.ReferenceIndex, Math.max(0, ctx.MantaBP2.Position + ctx.Uncertainty2.Start - range), ctx.MantaBP2.Position + ctx.Uncertainty2.End + range) });
    final File TEMP_REF_BAM = queryNameSortedBAM(refReader, intervals, "ref");
    final File TEMP_TUMOR_BAM = queryNameSortedBAM(tumorReader, intervals, "tumor");
    final SamReader SORTED_REF_READER = SamReaderFactory.makeDefault().open(TEMP_REF_BAM);
    final SamReader SORTED_TUMOR_READER = SamReaderFactory.makeDefault().open(TEMP_TUMOR_BAM);
    final BreakpointResult breakpoints = determineBreakpoints(ctx, SORTED_TUMOR_READER);
    final StructuralVariantResult result = new StructuralVariantResult();
    result.Breakpoints = breakpoints.Breakpoints;
    result.QueryIntervals = intervals;
    if (breakpoints.Error != BreakpointError.NONE) {
        result.Filters = Filter.getErrorFilter();
    } else {
        result.TumorStats = collectEvidence(ctx, SORTED_TUMOR_READER, result.Breakpoints);
        result.RefStats = collectEvidence(ctx, SORTED_REF_READER, result.Breakpoints);
        result.AlleleFrequency = AlleleFrequency.calculate(result.TumorStats);
        // load sample clipping
        SORTED_TUMOR_READER.forEach(r -> Clipping.getClips(r).forEach(c -> result.TumorStats.Sample_Clipping.add(c)));
        SORTED_REF_READER.forEach(r -> Clipping.getClips(r).forEach(c -> result.RefStats.Sample_Clipping.add(c)));
        result.Filters = Filter.getFilters(ctx, result.TumorStats, result.RefStats, result.Breakpoints, contamination);
        // adjust for homology
        final Location bp1 = result.Breakpoints.getLeft().add(ctx.OrientationBP1 > 0 ? 0 : -1);
        final Location bp2;
        if (!ctx.isInsert() && ctx.InsertSequence.isEmpty()) {
            bp2 = result.Breakpoints.getRight().add(-ctx.OrientationBP2 * ctx.HomologySequence.length()).add(ctx.OrientationBP2 > 0 ? 0 : -1);
        } else {
            bp2 = result.Breakpoints.getRight().add(ctx.OrientationBP2 > 0 ? 0 : -1);
        }
        result.Breakpoints = Pair.of(bp1, bp2);
    }
    result.FilterString = result.Filters.isEmpty() ? "PASS" : String.join(";", result.Filters);
    // clean up
    SORTED_REF_READER.close();
    SORTED_TUMOR_READER.close();
    if (!TEMP_REF_BAM.delete()) {
        LOGGER.error("couldn't delete {}", TEMP_REF_BAM);
    }
    if (!TEMP_TUMOR_BAM.delete()) {
        LOGGER.error("couldn't delete {}", TEMP_TUMOR_BAM);
    }
    return result;
}

Example 2 with Clipping

use of com.hartwig.hmftools.breakpointinspector.clipping.Clipping in project hmftools by hartwigmedical.

the class Analysis method determineBreakpointsImprecise.

private BreakpointResult determineBreakpointsImprecise(final HMFVariantContext ctx, final SamReader reader) {
    final Pair<Integer, Integer> ctxOrientation = Pair.of(ctx.OrientationBP1, ctx.OrientationBP2);
    final PairedReads interesting = new PairedReads();
    final PairedReads clipped_proper = new PairedReads();
    final PairedReads secondary_pairs = new PairedReads();
    final List<SAMRecord> currentReads = Lists.newArrayList();
    final SAMRecordIterator iterator = reader.iterator();
    while (iterator.hasNext() || !currentReads.isEmpty()) {
        final SAMRecord record = iterator.hasNext() ? iterator.next() : null;
        if (record != null) {
            if (currentReads.isEmpty() || record.getReadName().equals(currentReads.get(0).getReadName())) {
                currentReads.add(record);
                continue;
            }
        }
        currentReads.sort(comparator);
        final PairedReads pairs = pairs(currentReads);
        currentReads.clear();
        if (record != null) {
            currentReads.add(record);
        }
        for (final Pair<SAMRecord, SAMRecord> pair : pairs) {
            final boolean correctOrientation = orientation(pair).equals(ctxOrientation);
            final boolean correctChromosome = Location.fromSAMRecord(pair.getLeft()).sameChromosomeAs(ctx.MantaBP1) && Location.fromSAMRecord(pair.getRight()).sameChromosomeAs(ctx.MantaBP2);
            final boolean hasExpectedClipping = clippedOnCorrectSide(pair.getLeft(), ctx.OrientationBP1) || clippedOnCorrectSide(pair.getRight(), ctx.OrientationBP2);
            final boolean sameChromosome = pair.getLeft().getReferenceIndex().equals(pair.getRight().getReferenceIndex());
            final boolean potentialSROnly = sameChromosome && Stream.of(ctx.OrientationBP1, ctx.OrientationBP2).anyMatch(orientation -> clippedOnCorrectSide(orientation > 0 ? pair.getRight() : pair.getLeft(), orientation));
            final boolean secondary = stream(pair).anyMatch(SAMRecord::isSecondaryOrSupplementary);
            final boolean proper = stream(pair).anyMatch(SAMRecord::getProperPairFlag) && !secondary;
            LOGGER.trace("determineBreakpoints {} {}->{} {} {}->{} correctOrientation({}) correctChromosome({}) hasExpectedClipping({}) proper({}) potentialSROnly({}) secondary({})", pair.getLeft().getReadName(), pair.getLeft().getAlignmentStart(), pair.getLeft().getMateAlignmentStart(), pair.getRight().getReadName(), pair.getRight().getAlignmentStart(), pair.getRight().getMateAlignmentStart(), correctOrientation, correctChromosome, hasExpectedClipping, proper, potentialSROnly, secondary);
            if (secondary && potentialSROnly) {
                secondary_pairs.add(pair);
            } else if ((!proper || hasExpectedClipping) && correctChromosome && correctOrientation) {
                interesting.add(pair);
            } else if (proper && potentialSROnly) {
                clipped_proper.add(pair);
            }
        }
    }
    iterator.close();
    // load clipping info
    Clipping bp1_clipping = new Clipping();
    Clipping bp2_clipping = new Clipping();
    for (final Pair<SAMRecord, SAMRecord> pair : interesting) {
        if (ctx.OrientationBP1 > 0) {
            bp1_clipping.add(Clipping.getRightClip(pair.getLeft()));
        } else {
            bp1_clipping.add(Clipping.getLeftClip(pair.getLeft()));
        }
        if (ctx.OrientationBP2 > 0) {
            bp2_clipping.add(Clipping.getRightClip(pair.getRight()));
        } else {
            bp2_clipping.add(Clipping.getLeftClip(pair.getRight()));
        }
    }
    for (final Pair<SAMRecord, SAMRecord> pair : clipped_proper) {
        if (stream(pair).allMatch(r -> Location.fromSAMRecord(r).sameChromosomeAs(ctx.MantaBP1))) {
            if (ctx.OrientationBP1 > 0) {
                bp1_clipping.add(Clipping.getRightClip(pair.getRight()));
            } else {
                bp1_clipping.add(Clipping.getLeftClip(pair.getLeft()));
            }
        }
        if (stream(pair).allMatch(r -> Location.fromSAMRecord(r).sameChromosomeAs(ctx.MantaBP2))) {
            if (ctx.OrientationBP2 > 0) {
                bp2_clipping.add(Clipping.getRightClip(pair.getRight()));
            } else {
                bp2_clipping.add(Clipping.getLeftClip(pair.getLeft()));
            }
        }
    }
    for (final Pair<SAMRecord, SAMRecord> pair : secondary_pairs) {
        if (stream(pair).allMatch(r -> Location.fromSAMRecord(r).sameChromosomeAs(ctx.MantaBP1))) {
            if (ctx.OrientationBP1 > 0) {
                bp1_clipping.add(Clipping.getRightClip(pair.getRight()));
            } else {
                bp1_clipping.add(Clipping.getLeftClip(pair.getLeft()));
            }
        }
        if (stream(pair).allMatch(r -> Location.fromSAMRecord(r).sameChromosomeAs(ctx.MantaBP2))) {
            if (ctx.OrientationBP2 > 0) {
                bp2_clipping.add(Clipping.getRightClip(pair.getRight()));
            } else {
                bp2_clipping.add(Clipping.getLeftClip(pair.getLeft()));
            }
        }
    }
    // determinate candidates based on clipping info
    final List<Location> bp1_candidates = bp1_clipping.getSequences().stream().map(c -> c.Alignment).filter(c -> withinRange(c, ctx.MantaBP1, ctx.Uncertainty1)).collect(Collectors.toList());
    if (bp1_candidates.isEmpty()) {
        final Location candidate = interesting.stream().map(Pair::getLeft).map(r -> Location.fromSAMRecord(r, ctx.OrientationBP1 < 0).add(ctx.OrientationBP1 > 0 ? 1 : -1)).filter(l -> withinRange(l, ctx.MantaBP1, ctx.Uncertainty1)).max((a, b) -> ctx.OrientationBP1 > 0 ? a.compareTo(b) : b.compareTo(a)).orElse(null);
        if (candidate != null) {
            bp1_candidates.add(candidate);
        }
    }
    final List<Location> bp2_candidates = bp2_clipping.getSequences().stream().map(c -> c.Alignment).filter(c -> withinRange(c, ctx.MantaBP2, ctx.Uncertainty2)).collect(Collectors.toList());
    if (bp2_candidates.isEmpty()) {
        final Location candidate = interesting.stream().map(Pair::getRight).map(r -> Location.fromSAMRecord(r, ctx.OrientationBP2 < 0).add(ctx.OrientationBP2 > 0 ? 1 : -1)).filter(l -> withinRange(l, ctx.MantaBP2, ctx.Uncertainty2)).max((a, b) -> ctx.OrientationBP2 > 0 ? a.compareTo(b) : b.compareTo(a)).orElse(null);
        if (candidate != null) {
            bp2_candidates.add(candidate);
        }
    }
    // NOTE: we include homology on both sides here and take it out later
    LOGGER.trace("bp1_candidates={} bp2_candidates={}", bp1_candidates, bp2_candidates);
    final Location breakpoint1 = bp1_candidates.isEmpty() ? null : bp1_candidates.get(0).add(-ctx.OrientationBP1);
    final Location breakpoint2 = bp2_candidates.isEmpty() ? null : bp2_candidates.get(0).add(-ctx.OrientationBP2);
    return BreakpointResult.from(Pair.of(breakpoint1, breakpoint2));
}