Examples with CigarElement htsjdk.samtools.CigarElement used on opensource projects

Example 11 with CigarElement

use of htsjdk.samtools.CigarElement in project gatk by broadinstitute.

the class AlignmentUtils method isInsideDeletion.

/**
     * Is the offset inside a deletion?
     *
     * @param cigar         the read's CIGAR -- cannot be null
     * @param offset        the offset into the CIGAR
     * @return true if the offset is inside a deletion, false otherwise
     */
public static boolean isInsideDeletion(final Cigar cigar, final int offset) {
    Utils.nonNull(cigar);
    if (offset < 0)
        return false;
    // pos counts read bases
    int pos = 0;
    int prevPos = 0;
    for (final CigarElement ce : cigar.getCigarElements()) {
        switch(ce.getOperator()) {
            case I:
            case S:
            case D:
            case M:
            case EQ:
            case X:
                prevPos = pos;
                pos += ce.getLength();
                break;
            case H:
            case P:
            case N:
                break;
            default:
                throw new GATKException("Unsupported cigar operator: " + ce.getOperator());
        }
        // Is the offset inside a deletion?
        if (prevPos < offset && pos >= offset && ce.getOperator() == CigarOperator.D) {
            return true;
        }
    }
    return false;
}

Example 12 with CigarElement

use of htsjdk.samtools.CigarElement in project gatk by broadinstitute.

the class AlignmentUtils method getNumHardClippedBases.

/**
     * Count the number of bases hard clipped from read
     *
     * If read's cigar is null, return 0
     *
     * @param r a non-null read
     * @return a positive integer
     */
public static int getNumHardClippedBases(final GATKRead r) {
    if (r == null)
        throw new IllegalArgumentException("Read cannot be null");
    int n = 0;
    final Cigar cigar = r.getCigar();
    if (cigar == null)
        return 0;
    for (final CigarElement e : cigar.getCigarElements()) if (e.getOperator() == CigarOperator.H)
        n += e.getLength();
    return n;
}

Example 13 with CigarElement

use of htsjdk.samtools.CigarElement in project gatk by broadinstitute.

the class AlignmentUtils method moveCigarLeft.

/**
     * 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<>(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);
}

Example 14 with CigarElement

use of htsjdk.samtools.CigarElement in project gatk by broadinstitute.

the class AlignmentUtils method getNumAlignmentBlocks.

/**
     * Returns number of alignment blocks (continuous stretches of aligned bases) in the specified alignment.
     * This method follows closely the SAMRecord::getAlignmentBlocks() implemented in samtools library, but
     * it only counts blocks without actually allocating and filling the list of blocks themselves. Hence, this method is
     * a much more efficient alternative to r.getAlignmentBlocks.size() in the situations when this number is all that is needed.
     * Formally, this method simply returns the number of M elements in the cigar.
     *
     * @param r alignment
     * @return number of continuous alignment blocks (i.e. 'M' elements of the cigar; all indel and clipping elements are ignored).
     */
public static int getNumAlignmentBlocks(final GATKRead r) {
    Utils.nonNull(r);
    final Cigar cigar = r.getCigar();
    if (cigar == null)
        return 0;
    int n = 0;
    for (final CigarElement e : cigar.getCigarElements()) {
        if (ALIGNED_TO_GENOME_OPERATORS.contains(e.getOperator()))
            n++;
    }
    return n;
}

Example 15 with CigarElement

use of htsjdk.samtools.CigarElement in project gatk by broadinstitute.

the class AlignmentUtils method calcNumHighQualitySoftClips.

/**
     * Calculate the number of bases that are soft clipped in read with quality score greater than threshold
     *
     * Handles the case where the cigar is null (i.e., the read is unmapped), returning 0
     *
     * @param read a non-null read.
     * @param qualThreshold consider bases with quals > this value as high quality.  Must be >= 0
     * @return positive integer
     */
public static int calcNumHighQualitySoftClips(final GATKRead read, final byte qualThreshold) {
    if (read == null)
        throw new IllegalArgumentException("Read cannot be null");
    if (qualThreshold < 0)
        throw new IllegalArgumentException("Expected qualThreshold to be a positive byte but saw " + qualThreshold);
    if (// the read is unmapped
    read.getCigar() == null)
        return 0;
    final byte[] qual = read.getBaseQualities();
    int numHQSoftClips = 0;
    int alignPos = 0;
    for (final CigarElement ce : read.getCigarElements()) {
        final int elementLength = ce.getLength();
        switch(ce.getOperator()) {
            case S:
                for (int jjj = 0; jjj < elementLength; jjj++) {
                    if (qual[alignPos++] > qualThreshold) {
                        numHQSoftClips++;
                    }
                }
                break;
            case M:
            case I:
            case EQ:
            case X:
                alignPos += elementLength;
                break;
            case H:
            case P:
            case D:
            case N:
                break;
            default:
                throw new IllegalStateException("Unsupported cigar operator: " + ce.getOperator());
        }
    }
    return numHQSoftClips;
}