Java ReadClipper類代碼示例

import org.broadinstitute.hellbender.utils.clipping.ReadClipper; //導入依賴的package包/類
@Override
public void apply( GATKRead read, ReferenceContext ref, FeatureContext featureContext ) {
    if ( onlyDoRead == null || read.getName().equals(onlyDoRead) ) {
        if ( clippingRepresentation == ClippingRepresentation.HARDCLIP_BASES || clippingRepresentation == ClippingRepresentation.REVERT_SOFTCLIPPED_BASES )
            read = ReadClipper.revertSoftClippedBases(read);
        ReadClipperWithData clipper = new ReadClipperWithData(read, sequencesToClip);

        //
        // run all three clipping modules
        //
        clipBadQualityScores(clipper);
        clipCycles(clipper);
        clipSequences(clipper);
        accumulate(clipper);
    }
}
 

import org.broadinstitute.hellbender.utils.clipping.ReadClipper; //導入依賴的package包/類
@Test
public void testSimpleContexts() {
    final Random rnd = Utils.getRandomGenerator();
    final SAMFileHeader header = ArtificialReadUtils.createArtificialSamHeader();

    for(int i = 0; i < 10; i++) {
        final GATKRead read = ArtificialReadUtils.createRandomRead(header, 1000);
        read.setIsReverseStrand(rnd.nextBoolean());
        final GATKRead clippedRead = ReadClipper.clipLowQualEnds(read, RAC.LOW_QUAL_TAIL, ClippingRepresentation.WRITE_NS);
        final ReadCovariates readCovariates = new ReadCovariates(read.getLength(), 1, new CovariateKeyCache());
        covariate.recordValues(read, header, readCovariates, true);

        verifyCovariateArray(readCovariates.getMismatchesKeySet(), RAC.MISMATCHES_CONTEXT_SIZE, clippedRead, covariate, RAC.LOW_QUAL_TAIL);
        verifyCovariateArray(readCovariates.getInsertionsKeySet(), RAC.INDELS_CONTEXT_SIZE, clippedRead, covariate, RAC.LOW_QUAL_TAIL);
        verifyCovariateArray(readCovariates.getDeletionsKeySet(), RAC.INDELS_CONTEXT_SIZE, clippedRead, covariate, RAC.LOW_QUAL_TAIL);
    }
}
 

import org.broadinstitute.hellbender.utils.clipping.ReadClipper; //導入依賴的package包/類
/**
 * Correct a collection of reads based on stored k-mer counts
 * @param reads
 */
public final List<GATKRead> correctReads(final Collection<GATKRead> reads) {

    final List<GATKRead> correctedReads = new ArrayList<>(reads.size());
    if (DONT_CORRECT_IN_LONG_HOMOPOLYMERS && maxHomopolymerLengthInRegion > MAX_HOMOPOLYMER_THRESHOLD) {
        // just copy reads into output and exit
        correctedReads.addAll(reads);
    }
    else {
        computeKmerCorrectionMap();
        for (final GATKRead read: reads) {
            final GATKRead correctedRead = correctRead(read);
            if (trimLowQualityBases) {
                correctedReads.add(ReadClipper.hardClipLowQualEnds(correctedRead, minTailQuality));
            } else {
                correctedReads.add(correctedRead);
            }
        }
        if (debug) {
            logger.info("Number of corrected bases:" + readErrorCorrectionStats.numBasesCorrected);
            logger.info("Number of corrected reads:" + readErrorCorrectionStats.numReadsCorrected);
            logger.info("Number of skipped reads:" + readErrorCorrectionStats.numReadsUncorrected);
            logger.info("Number of solid kmers:" + readErrorCorrectionStats.numSolidKmers);
            logger.info("Number of corrected kmers:" + readErrorCorrectionStats.numCorrectedKmers);
            logger.info("Number of uncorrectable kmers:" + readErrorCorrectionStats.numUncorrectableKmers);
        }
    }
    return correctedReads;
}
 

import org.broadinstitute.hellbender.utils.clipping.ReadClipper; //導入依賴的package包/類
private static GATKRead handleTrimmed(final GATKRead read) {
    // store the start and the end, and remove the tags
    // the tags should be removed here because the completely trimmed flag is set
    // before removing and thus if the read was trimmed the length is not longer the same
    final int start = RTReadUtils.getTrimmingStartPoint(read);
    final int end = RTReadUtils.getTrimmingEndPoint(read);
    RTReadUtils.clearTrimmingPointTags(read);

    if (read.isUnmapped()) {
        final byte[] newBases = Arrays.copyOfRange(read.getBases(), start, end);
        final byte[] newQuals = Arrays.copyOfRange(read.getBaseQualities(), start, end);
        read.setBases(newBases);
        read.setBaseQualities(newQuals);
    } else {
        final int readLength = read.getLength();
        final ReadClipper clipper = new ReadClipper(read);
        // we have to clip the end first, because the read clipper works in a step-wise way
        if (end != readLength) {
            clipper.addOp(new ClippingOp(end, readLength));
        }
        if (start != 0) {
            // we should remove 1 for correctly clipping
            clipper.addOp(new ClippingOp(0, start - 1));
        }
        // this should modify in place here
        modifyWithClipped(read, clipper.clipRead(ClippingRepresentation.HARDCLIP_BASES));
        // return the same read
    }

    // it is modified in place
    return read;
}
 

import org.broadinstitute.hellbender.utils.clipping.ReadClipper; //導入依賴的package包/類
private GATKRead clipReadRightEnd(GATKRead read) {
    final byte[] quals = read.getBaseQualities();
    final int clipPoint = getRightClipPoint(quals);
    if (clipPoint != -1) {
        final ReadClipper readClipper = new ReadClipper(read);
        readClipper.addOp(new ClippingOp(clipPoint, read.getLength()));
        return readClipper.clipRead(ClippingRepresentation.HARDCLIP_BASES);
    } else {
        return read;
    }
}
 

import org.broadinstitute.hellbender.utils.clipping.ReadClipper; //導入依賴的package包/類
private GATKRead clipReadLeftEnd(GATKRead read) {
    final byte[] quals = read.getBaseQualities();
    final int clipPoint = getLeftClipPoint(quals);
    if (clipPoint != -1) {
        final ReadClipper readClipper = new ReadClipper(read);
        readClipper.addOp(new ClippingOp(0, clipPoint));
        return readClipper.clipRead(ClippingRepresentation.HARDCLIP_BASES);
    } else {
        return read;
    }
}
 

import org.broadinstitute.hellbender.utils.clipping.ReadClipper; //導入依賴的package包/類
@Override
public GATKRead apply(final GATKRead read) {
    if (read.getLength() < minClipLength) return read;
    final byte[] bases = read.getBases();

    //Be sure to find the best match in case one adapter ends with another adapter's sequence (within maxMismatches difference)
    int bestNumMismatches = maxMismatches;
    int bestAlignmentStart = bases.length;
    boolean foundAdapter = false;

    //Adapter list loop
    for (final String adapterSequence : adapterSequences) {
        //Start at the end of the read and walk backwards
        int alignmentLength = minClipLength;
        for (int alignmentStart = bases.length - minClipLength; alignmentStart >= 0; alignmentStart--) {
            int numMismatches = 0;
            //Check each base for mismatches
            for (int j = 0; j < alignmentLength; j++) {
                if (!SequenceUtil.isNoCall((byte) adapterSequence.charAt(j)) && bases[alignmentStart + j] != adapterSequence.charAt(j)) {
                    if (++numMismatches > maxMismatches) break;
                }
            }
            if (numMismatches < bestNumMismatches || (numMismatches == bestNumMismatches && alignmentStart < bestAlignmentStart)) {
                //We have a (better/earlier) match
                bestNumMismatches = numMismatches;
                bestAlignmentStart = alignmentStart;
                foundAdapter = true;
            }
            alignmentLength = alignmentLength < adapterSequence.length() ? alignmentLength + 1 : alignmentLength;
        }
    }
    if (foundAdapter) {
        //Hard clip from the beginning of the adapter to the end of the read
        final ReadClipper readClipper = new ReadClipper(read);
        readClipper.addOp(new ClippingOp(bestAlignmentStart, read.getLength()));
        return readClipper.clipRead(ClippingRepresentation.HARDCLIP_BASES);
    }
    return read;
}
 

import org.broadinstitute.hellbender.utils.clipping.ReadClipper; //導入依賴的package包/類
/**
 * Trim this region to no more than the span, producing a new assembly region with properly trimmed reads that
 * attempts to provide the best possible representation of this region covering the span.
 *
 * The challenge here is that span may (1) be larger than can be represented by this assembly region
 * + its original extension and (2) the extension must be symmetric on both sides.  This algorithm
 * therefore determines how best to represent span as a subset of the span of this
 * region with a padding value that captures as much of the span as possible.
 *
 * For example, suppose this active region is
 *
 * Active:    100-200 with extension of 50, so that the true span is 50-250
 * NewExtent: 150-225 saying that we'd ideally like to just have bases 150-225
 *
 * Here we represent the assembly region as a region from 150-200 with 25 bp of padding.
 *
 * The overall constraint is that the region can never exceed the original region, and
 * the extension is chosen to maximize overlap with the desired region
 *
 * @param span the new extend of the active region we want
 * @return a non-null, empty active region
 */
public AssemblyRegion trim(final SimpleInterval span, final SimpleInterval extendedSpan) {
    Utils.nonNull(span, "Active region extent cannot be null");
    Utils.nonNull(extendedSpan, "Active region extended span cannot be null");
    Utils.validateArg(extendedSpan.contains(span), "The requested extended span must fully contain the requested span");

    final SimpleInterval subActive = getSpan().intersect(span);
    final int requiredOnRight = Math.max(extendedSpan.getEnd() - subActive.getEnd(), 0);
    final int requiredOnLeft = Math.max(subActive.getStart() - extendedSpan.getStart(), 0);
    final int requiredExtension = Math.min(Math.max(requiredOnLeft, requiredOnRight), getExtension());

    final AssemblyRegion result = new AssemblyRegion( subActive, Collections.<ActivityProfileState>emptyList(), isActive, requiredExtension, header );

    final List<GATKRead> myReads = getReads();
    final SimpleInterval resultExtendedLoc = result.getExtendedSpan();
    final int resultExtendedLocStart = resultExtendedLoc.getStart();
    final int resultExtendedLocStop = resultExtendedLoc.getEnd();

    final List<GATKRead> trimmedReads = new ArrayList<>(myReads.size());
    for( final GATKRead read : myReads ) {
        final GATKRead clippedRead = ReadClipper.hardClipToRegion(read, resultExtendedLocStart, resultExtendedLocStop);
        if( result.readOverlapsRegion(clippedRead) && clippedRead.getLength() > 0 ) {
            trimmedReads.add(clippedRead);
        }
    }
    result.clearReads();

    trimmedReads.sort(new ReadCoordinateComparator(header));
    result.addAll(trimmedReads);
    return result;
}
 

import org.broadinstitute.hellbender.utils.clipping.ReadClipper; //導入依賴的package包/類
private ReadTransformer makeReadTransform() {
    ReadTransformer f0 = BaseRecalibrationEngine::consolidateCigar;

    ReadTransformer f = f0.andThen(this::setDefaultBaseQualities)
            .andThen(this::resetOriginalBaseQualities)
            .andThen(ReadClipper::hardClipAdaptorSequence)
            .andThen(ReadClipper::hardClipSoftClippedBases);

    return f;
}
 

import org.broadinstitute.hellbender.utils.clipping.ReadClipper; //導入依賴的package包/類
/**
 * Given a read, clips low quality ends (by overwriting with N) and returns the underlying bases, after
 * reverse-complementing for negative-strand reads.
 * @param read the read
 * @param lowQTail every base quality lower than or equal to this in the tail of the read will be replaced with N.
 * @return bases of the read.
 */
@VisibleForTesting
static byte[] getStrandedClippedBytes(final GATKRead read, final byte lowQTail) {

    // Write N's over the low quality tail of the reads to avoid adding them into the context
    final GATKRead clippedRead = ReadClipper.clipLowQualEnds(read, lowQTail, ClippingRepresentation.WRITE_NS);

    final byte[] bases = clippedRead.getBases();
    if (read.isReverseStrand()) {
        return BaseUtils.simpleReverseComplement(bases);
    } else {
        return bases;
    }
}
 

import org.broadinstitute.hellbender.utils.clipping.ReadClipper; //導入依賴的package包/類
/**
 * Pull out an individual split position for a read
 *
 * @param read               the read being split
 * @param cigarStartIndex    the index of the first cigar element to keep
 * @param cigarEndIndex      the index of the last cigar element to keep
 * @param forSplitPositions  the manager for keeping track of split positions; can be null
 * @return a non-null read representing the section of the original read being split out
 */
private static GATKRead splitReadBasedOnCigar(final GATKRead read, final int cigarStartIndex, final int cigarEndIndex, final OverhangFixingManager forSplitPositions) {
    int cigarFirstIndex = cigarStartIndex;
    int cigarSecondIndex = cigarEndIndex;

    //in case a section of the read ends or starts with D (for example the first section in 1M1D1N1M is 1M1D), we should trim this cigar element
    // it can be 'if', but it was kept as 'while' to make sure the code can work with Cigar strings that were not "cleaned"
    while(read.getCigarElement(cigarFirstIndex).getOperator().equals(CigarOperator.D)) {
        cigarFirstIndex++;
    }
    while(read.getCigarElement(cigarSecondIndex-1).getOperator().equals(CigarOperator.D)) {
        cigarSecondIndex--;
    }
    if(cigarFirstIndex > cigarSecondIndex) {
        throw new IllegalArgumentException("Cannot split this read (might be an empty section between Ns, for example 1N1D1N): " + read.getCigar().toString());
    }

    // we keep only the section of the read that is aligned to the reference between startRefIndex and stopRefIndex (inclusive).
    // the other sections of the read are clipped:
    final int startRefIndex = read.getUnclippedStart() + CigarUtils.countRefBasesBasedOnCigar(read, 0, cigarFirstIndex); //goes through the prefix of the cigar (up to cigarStartIndex) and move the reference index.
    final int stopRefIndex = startRefIndex + CigarUtils.countRefBasesBasedOnCigar(read,cigarFirstIndex,cigarSecondIndex)-1; //goes through a consecutive non-N section of the cigar (up to cigarEndIndex) and move the reference index.

    if ( forSplitPositions != null ) {
        final String contig = read.getContig();
        final int splitStart = startRefIndex + CigarUtils.countRefBasesBasedOnCigar(read,cigarFirstIndex,cigarEndIndex);  //we use cigarEndIndex instead of cigarSecondIndex so we won't take into account the D's at the end.
        final int splitEnd = splitStart + read.getCigarElement(cigarEndIndex).getLength() - 1;
        forSplitPositions.addSplicePosition(contig, splitStart, splitEnd);
    }

    return ReadClipper.softClipToRegionIncludingClippedBases(read,startRefIndex,stopRefIndex);
}
 

import org.broadinstitute.hellbender.utils.clipping.ReadClipper; //導入依賴的package包/類
@Override
void fixSplit(SplitRead read, Splice splice) {
    if ( read.read.getStart() == 10000) {
        final GATKRead clippedRead = ReadClipper.softClipByReadCoordinates(read.read, 0, 2);
        read.setRead(clippedRead);
    }
}
 

import org.broadinstitute.hellbender.utils.clipping.ReadClipper; //導入依賴的package包/類
public static void finalizeRegion(final AssemblyRegion region,
                                  final boolean errorCorrectReads,
                                  final boolean dontUseSoftClippedBases,
                                  final byte minTailQuality,
                                  final SAMFileHeader readsHeader,
                                  final SampleList samplesList) {
    if ( region.isFinalized() ) {
        return;
    }

    // Loop through the reads hard clipping the adaptor and low quality tails
    final List<GATKRead> readsToUse = new ArrayList<>(region.getReads().size());
    for( final GATKRead myRead : region.getReads() ) {
        final byte minTailQualityToUse = errorCorrectReads ? HaplotypeCallerEngine.MIN_TAIL_QUALITY_WITH_ERROR_CORRECTION : minTailQuality;
        GATKRead clippedRead = ReadClipper.hardClipLowQualEnds(myRead, minTailQualityToUse);

        // remove soft clips if we cannot reliably clip off adapter sequence or if the user doesn't want to use soft clips at all
        // otherwie revert soft clips so that we see the alignment start and end assuming the soft clips are all matches
        // TODO -- WARNING -- still possibility that unclipping the soft clips will introduce bases that aren't
        // TODO -- truly in the extended region, as the unclipped bases might actually include a deletion
        // TODO -- w.r.t. the reference.  What really needs to happen is that kmers that occur before the
        // TODO -- reference haplotype start must be removed
        clippedRead = dontUseSoftClippedBases || ! ReadUtils.hasWellDefinedFragmentSize(clippedRead) ?
                ReadClipper.hardClipSoftClippedBases(clippedRead) : ReadClipper.revertSoftClippedBases(clippedRead);

        clippedRead = clippedRead.isUnmapped() ? clippedRead : ReadClipper.hardClipAdaptorSequence(clippedRead);
        if ( ! clippedRead.isEmpty() && clippedRead.getCigar().getReadLength() > 0 ) {
            clippedRead = ReadClipper.hardClipToRegion( clippedRead, region.getExtendedSpan().getStart(), region.getExtendedSpan().getEnd() );
            if ( region.readOverlapsRegion(clippedRead) && clippedRead.getLength() > 0 ) {
                readsToUse.add(clippedRead);
            }
        }
    }

    // TODO -- Performance optimization: we partition the reads by sample 4 times right now; let's unify that code.
    // final List<GATKRead> downsampledReads = DownsamplingUtils.levelCoverageByPosition(ReadUtils.sortReadsByCoordinate(readsToUse), maxReadsInRegionPerSample, minReadsPerAlignmentStart);
    Collections.sort(readsToUse, new ReadCoordinateComparator(readsHeader)); // TODO: sort may be unnecessary here

    // handle overlapping read pairs from the same fragment
    cleanOverlappingReadPairs(readsToUse, samplesList, readsHeader);

    region.clearReads();
    region.addAll(readsToUse);
    region.setFinalized(true);
}
 

import org.broadinstitute.hellbender.utils.clipping.ReadClipper; //導入依賴的package包/類
private GATKRead maskRead(final GATKRead read, final int maskStart, final int maskEnd) {
    final ReadClipper readClipper = new ReadClipper(read);
    readClipper.addOp(new ClippingOp(maskStart, maskEnd - 1));
    return readClipper.clipRead(ClippingRepresentation.WRITE_NS_Q0S);
}