Java CigarElement类代码示例

import htsjdk.samtools.CigarElement; //导入依赖的package包/类
/**
 * Create a new pileup element
 *
 * @param read a non-null read to pileup
 * @param baseOffset the offset into the read's base / qual vector aligned to this position on the genome. If the
 *                   current cigar element is a deletion, offset should be the offset of the last M/=/X position.
 * @param currentElement a non-null CigarElement that indicates the cigar element aligning the read to the genome
 * @param currentCigarOffset the offset of currentElement in read.getCigar().getElement(currentCigarOffset) == currentElement)
 * @param offsetInCurrentCigar how far into the currentElement are we in our alignment to the genome?
 */
public PileupElement(final GATKSAMRecord read, final int baseOffset,
                     final CigarElement currentElement, final int currentCigarOffset,
                     final int offsetInCurrentCigar) {
    assert currentElement != null;

    this.read = read;
    this.offset = baseOffset;
    this.currentCigarElement = currentElement;
    this.currentCigarOffset = currentCigarOffset;
    this.offsetInCurrentCigar = offsetInCurrentCigar;

    // for performance regions these are assertions
    assert this.read != null;
    assert this.offset >= 0 && this.offset < this.read.getReadLength();
    assert this.currentCigarOffset >= 0;
    assert this.currentCigarOffset < read.getCigarLength();
    assert this.offsetInCurrentCigar >= 0;
    assert this.offsetInCurrentCigar < currentElement.getLength();
}
 

import htsjdk.samtools.CigarElement; //导入依赖的package包/类
/**
 * Helper function to get cigar elements between this and either the prev or next genomic position
 *
 * @param direction PREVIOUS if we want before, NEXT if we want after
 * @return a non-null list of cigar elements between this and the neighboring position in direction
 */
private LinkedList<CigarElement> getBetween(final Direction direction) {
    final int increment = direction == Direction.NEXT ? 1 : -1;
    LinkedList<CigarElement> elements = null;
    final int nCigarElements = read.getCigarLength();
    for ( int i = currentCigarOffset + increment; i >= 0 && i < nCigarElements; i += increment) {
        final CigarElement elt = read.getCigar().getCigarElement(i);
        if ( ON_GENOME_OPERATORS.contains(elt.getOperator()) )
            break;
        else {
            // optimization: don't allocate list if not necessary
            if ( elements == null )
                elements = new LinkedList<CigarElement>();

            if ( increment > 0 )
                // to keep the list in the right order, if we are incrementing positively add to the end
                elements.add(elt);
            else
                // counting down => add to front
                elements.addFirst(elt);
        }
    }

    // optimization: elements is null because nothing got added, just return the empty list
    return elements == null ? EMPTY_LINKED_LIST : elements;
}
 

import htsjdk.samtools.CigarElement; //导入依赖的package包/类
/**
 * Calculates the reference coordinate for the end of the read taking into account soft clips but not hard clips.
 *
 * Note: getUnclippedEnd() adds soft and hard clips, this function only adds soft clips.
 *
 * @return the unclipped end of the read taking soft clips (but not hard clips) into account
 */
public int getSoftEnd() {
    if ( softEnd == UNINITIALIZED ) {
        boolean foundAlignedBase = false;
        softEnd = getAlignmentEnd();
        final List<CigarElement> cigs = getCigar().getCigarElements();
        for (int i = cigs.size() - 1; i >= 0; --i) {
            final CigarElement cig = cigs.get(i);
            final CigarOperator op = cig.getOperator();

            if (op == CigarOperator.SOFT_CLIP) // assumes the soft clip that we found is at the end of the aligned read
                softEnd += cig.getLength();
            else if (op != CigarOperator.HARD_CLIP) {
                foundAlignedBase = true;
                break;
            }
        }
        if( !foundAlignedBase ) { // for example 64H14S, the soft end is actually the same as the alignment end
            softEnd = getAlignmentEnd();
        }
    }

    return softEnd;
}
 

import htsjdk.samtools.CigarElement; //导入依赖的package包/类
/**
 * Combines equal adjacent elements of a Cigar object
 *
 * @param rawCigar the cigar object
 * @return a combined cigar object
 */
public static Cigar combineAdjacentCigarElements(Cigar rawCigar) {
    Cigar combinedCigar = new Cigar();
    CigarElement lastElement = null;
    int lastElementLength = 0;
    for (CigarElement cigarElement : rawCigar.getCigarElements()) {
        if (lastElement != null &&
                ((lastElement.getOperator() == cigarElement.getOperator()) ||
                        (lastElement.getOperator() == CigarOperator.I && cigarElement.getOperator() == CigarOperator.D) ||
                        (lastElement.getOperator() == CigarOperator.D && cigarElement.getOperator() == CigarOperator.I)))
            lastElementLength += cigarElement.getLength();
        else
        {
            if (lastElement != null)
                combinedCigar.add(new CigarElement(lastElementLength, lastElement.getOperator()));

            lastElement = cigarElement;
            lastElementLength = cigarElement.getLength();
        }
    }
    if (lastElement != null)
        combinedCigar.add(new CigarElement(lastElementLength, lastElement.getOperator()));

    return combinedCigar;
}
 

import htsjdk.samtools.CigarElement; //导入依赖的package包/类
public static final int countRefBasesBasedOnCigar(final GATKSAMRecord read, final int cigarStartIndex, final int cigarEndIndex){
    int result = 0;
    for(int i = cigarStartIndex; i<cigarEndIndex;i++){
        final CigarElement cigarElement = read.getCigar().getCigarElement(i);
        switch (cigarElement.getOperator()) {
            case M:
            case S:
            case D:
            case N:
            case H:
                result += cigarElement.getLength();
                break;
            case I:
                break;
            default:
                throw new ReviewedGATKException("Unsupported cigar operator: " + cigarElement.getOperator());
        }
    }
    return result;
}
 

import htsjdk.samtools.CigarElement; //导入依赖的package包/类
/**
     * Make sure that the SW didn't fail in some terrible way, and throw exception if it did
     */
    private static boolean isSWFailure(final SmithWaterman alignment) {
        // check that the alignment starts at the first base, which it should given the padding
        if ( alignment.getAlignmentStart2wrt1() > 0 ) {
            return true;
//          throw new IllegalStateException("SW failure ref " + paddedRef + " vs. " + paddedPath + " should always start at 0, but got " + alignment.getAlignmentStart2wrt1() + " with cigar " + alignment.getCigar());
        }

        // check that we aren't getting any S operators (which would be very bad downstream)
        for ( final CigarElement ce : alignment.getCigar().getCigarElements() ) {
            if ( ce.getOperator() == CigarOperator.S )
                return true;
            // soft clips at the end of the alignment are really insertions
//                throw new IllegalStateException("SW failure ref " + paddedRef + " vs. " + paddedPath + " should never contain S operators but got cigar " + alignment.getCigar());
        }

        return false;
    }
 

import htsjdk.samtools.CigarElement; //导入依赖的package包/类
/**
 * Move the indel in a given cigar string one base to the left
 *
 * @param cigar          original cigar
 * @param indexOfIndel   the index of the indel cigar element
 * @return non-null cigar with indel moved one base to the left
 */
private static Cigar moveCigarLeft(Cigar cigar, int indexOfIndel) {
    // get the first few elements
    ArrayList<CigarElement> elements = new ArrayList<CigarElement>(cigar.numCigarElements());
    for (int i = 0; i < indexOfIndel - 1; i++)
        elements.add(cigar.getCigarElement(i));

    // get the indel element and move it left one base
    CigarElement ce = cigar.getCigarElement(indexOfIndel - 1);
    elements.add(new CigarElement(Math.max(ce.getLength() - 1, 0), ce.getOperator()));
    elements.add(cigar.getCigarElement(indexOfIndel));
    if (indexOfIndel + 1 < cigar.numCigarElements()) {
        ce = cigar.getCigarElement(indexOfIndel + 1);
        elements.add(new CigarElement(ce.getLength() + 1, ce.getOperator()));
    } else {
        elements.add(new CigarElement(1, CigarOperator.M));
    }

    // get the last few elements
    for (int i = indexOfIndel + 2; i < cigar.numCigarElements(); i++)
        elements.add(cigar.getCigarElement(i));
    return new Cigar(elements);
}
 

import htsjdk.samtools.CigarElement; //导入依赖的package包/类
/**
 * Hard clips any leading insertions in the read. Only looks at the beginning of the read, not the end.
 *
 * @return a new read without leading insertions
 */
private GATKSAMRecord hardClipLeadingInsertions() {
    if (read.isEmpty())
        return read;

    for(CigarElement cigarElement : read.getCigar().getCigarElements()) {
        if (cigarElement.getOperator() != CigarOperator.HARD_CLIP && cigarElement.getOperator() != CigarOperator.SOFT_CLIP &&
                cigarElement.getOperator() != CigarOperator.INSERTION)
            break;

        else if (cigarElement.getOperator() == CigarOperator.INSERTION)
            this.addOp(new ClippingOp(0, cigarElement.getLength() - 1));

    }
    return clipRead(ClippingRepresentation.HARDCLIP_BASES);
}
 

import htsjdk.samtools.CigarElement; //导入依赖的package包/类
/**
 * Compute the offset of the first "real" position in the cigar on the genome
 *
 * This is defined as a first position after a run of Hs followed by a run of Ss
 *
 * @param cigar A non-null cigar
 * @return the offset (from 0) of the first on-genome base
 */
private int calcHardSoftOffset(final Cigar cigar) {
    final List<CigarElement> elements = cigar.getCigarElements();

    int size = 0;
    int i = 0;
    while ( i < elements.size() && elements.get(i).getOperator() == CigarOperator.HARD_CLIP ) {
        size += elements.get(i).getLength();
        i++;
    }
    while ( i < elements.size() && elements.get(i).getOperator() == CigarOperator.SOFT_CLIP ) {
        size += elements.get(i).getLength();
        i++;
    }

    return size;
}
 

import htsjdk.samtools.CigarElement; //导入依赖的package包/类
public static int getReadCoordinateForReferenceCoordinate(final int alignmentStart, final Cigar cigar, final int refCoord, final ClippingTail tail, final boolean allowGoalNotReached) {
    Pair<Integer, Boolean> result = getReadCoordinateForReferenceCoordinate(alignmentStart, cigar, refCoord, allowGoalNotReached);
    int readCoord = result.getFirst();

    // Corner case one: clipping the right tail and falls on deletion, move to the next
    // read coordinate. It is not a problem for the left tail because the default answer
    // from getReadCoordinateForReferenceCoordinate is to give the previous read coordinate.
    if (result.getSecond() && tail == ClippingTail.RIGHT_TAIL)
        readCoord++;

    // clipping the left tail and first base is insertion, go to the next read coordinate
    // with the same reference coordinate. Advance to the next cigar element, or to the
    // end of the read if there is no next element.
    final CigarElement firstElementIsInsertion = readStartsWithInsertion(cigar);
    if (readCoord == 0 && tail == ClippingTail.LEFT_TAIL && firstElementIsInsertion != null)
        readCoord = Math.min(firstElementIsInsertion.getLength(), cigar.getReadLength() - 1);

    return readCoord;
}
 

import htsjdk.samtools.CigarElement; //导入依赖的package包/类
/**
 * Calculates the reference coordinate for a read coordinate
 *
 * @param read   the read
 * @param offset the base in the read (coordinate in the read)
 * @return the reference coordinate correspondent to this base
 */
public static long getReferenceCoordinateForReadCoordinate(GATKSAMRecord read, int offset) {
    if (offset > read.getReadLength())
        throw new ReviewedGATKException(String.format(OFFSET_OUT_OF_BOUNDS_EXCEPTION, offset, read.getReadLength()));

    long location = read.getAlignmentStart();
    Iterator<CigarElement> cigarElementIterator = read.getCigar().getCigarElements().iterator();
    while (offset > 0 && cigarElementIterator.hasNext()) {
        CigarElement cigarElement = cigarElementIterator.next();
        long move = 0;
        if (cigarElement.getOperator().consumesReferenceBases())
            move = (long) Math.min(cigarElement.getLength(), offset);
        location += move;
        offset -= move;
    }
    if (offset > 0 && !cigarElementIterator.hasNext())
        throw new ReviewedGATKException(OFFSET_NOT_ZERO_EXCEPTION);

    return location;
}
 

import htsjdk.samtools.CigarElement; //导入依赖的package包/类
/**
 * Creates a CIGAR for the case where the prefix/suffix match combination encompasses an entire sequence
 *
 * @param reference            the reference sequence
 * @param alternate            the alternate sequence
 * @param matchingPrefix       the prefix match size
 * @param matchingSuffix       the suffix match size
 * @return non-null CIGAR
 */
private Cigar makeCigarForStrictPrefixAndSuffix(final byte[] reference, final byte[] alternate, final int matchingPrefix, final int matchingSuffix) {

    final List<CigarElement> result = new ArrayList<CigarElement>();

    // edge case: no D or I element
    if ( reference.length == alternate.length ) {
        result.add(makeElement(State.MATCH, matchingPrefix + matchingSuffix));
    } else {
        // add the first M element
        if ( matchingPrefix > 0 )
            result.add(makeElement(State.MATCH, matchingPrefix));

        // add the D or I element
        if ( alternate.length > reference.length )
            result.add(makeElement(State.INSERTION, alternate.length - reference.length));
        else // if ( reference.length > alternate.length )
            result.add(makeElement(State.DELETION, reference.length - alternate.length));

        // add the last M element
        if ( matchingSuffix > 0 )
            result.add(makeElement(State.MATCH, matchingSuffix));
    }

    return new Cigar(result);
}
 

import htsjdk.samtools.CigarElement; //导入依赖的package包/类
protected static Cigar reclipCigar(Cigar cigar, SAMRecord read) {
    ArrayList<CigarElement> elements = new ArrayList<CigarElement>();

    int i = 0;
    int n = read.getCigar().numCigarElements();
    while (i < n && isClipOperator(read.getCigar().getCigarElement(i).getOperator()))
        elements.add(read.getCigar().getCigarElement(i++));

    elements.addAll(cigar.getCigarElements());

    i++;
    while (i < n && !isClipOperator(read.getCigar().getCigarElement(i).getOperator()))
        i++;

    while (i < n && isClipOperator(read.getCigar().getCigarElement(i).getOperator()))
        elements.add(read.getCigar().getCigarElement(i++));

    return new Cigar(elements);
}
 

import htsjdk.samtools.CigarElement; //导入依赖的package包/类
private int getNumClippedBasesAtStart(final SAMRecord read) {
    // compute total number of clipped bases (soft or hard clipped)
    // check for hard clips (never consider these bases):
    final Cigar c = read.getCigar();
    final CigarElement first = c.getCigarElement(0);

    int numStartClippedBases = 0;
    if (first.getOperator() == CigarOperator.H) {
        numStartClippedBases = first.getLength();
    }
    final byte[] unclippedReadBases = read.getReadBases();
    final byte[] unclippedReadQuals = read.getBaseQualities();

    // Do a stricter base clipping than provided by CIGAR string, since this one may be too conservative,
    // and may leave a string of Q2 bases still hanging off the reads.
    for (int i = numStartClippedBases; i < unclippedReadBases.length; i++) {
        if (unclippedReadQuals[i] < PairHMMIndelErrorModel.BASE_QUAL_THRESHOLD)
            numStartClippedBases++;
        else
            break;

    }

    return numStartClippedBases;
}
 

import htsjdk.samtools.CigarElement; //导入依赖的package包/类
private int getNumClippedBasesAtEnd(final GATKSAMRecord read) {
    // compute total number of clipped bases (soft or hard clipped)
    // check for hard clips (never consider these bases):
    final Cigar c = read.getCigar();
    CigarElement last = c.getCigarElement(c.numCigarElements() - 1);

    int numEndClippedBases = 0;
    if (last.getOperator() == CigarOperator.H) {
        numEndClippedBases = last.getLength();
    }
    final byte[] unclippedReadBases = read.getReadBases();
    final byte[] unclippedReadQuals = read.getBaseQualities();

    // Do a stricter base clipping than provided by CIGAR string, since this one may be too conservative,
    // and may leave a string of Q2 bases still hanging off the reads.
    for (int i = unclippedReadBases.length - numEndClippedBases - 1; i >= 0; i--) {
        if (unclippedReadQuals[i] < PairHMMIndelErrorModel.BASE_QUAL_THRESHOLD)
            numEndClippedBases++;
        else
            break;
    }

    return numEndClippedBases;
}