Search in sources :

Example 26 with CigarElement

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

the class CigarUtils method countRefBasesBasedOnCigar.

/**
     * Compute the number of reference bases between the start (inclusive) and end (exclusive) cigar elements.
     * Note: The method does NOT use CigarOperator.consumesReferenceBases, since it checks something different.
     * The idea is you remove some elements from the beginning of the cigar string,
     * and want to recalculate what if the new starting reference position,
     * you want to count all the elements that indicate existing bases in the reference
     * (everything but I and P).
     * For example original position = 10. cigar: 2M3I2D1M. If you remove the 2M the new starting position is 12.
     * If you remove the 2M3I it is still 12. If you remove 2M3I2D (not reasonable cigar), you will get position 14.
     */
@SuppressWarnings("fallthru")
public static int countRefBasesBasedOnCigar(final GATKRead read, final int cigarStartIndex, final int cigarEndIndex) {
    if (read == null) {
        throw new IllegalArgumentException("null read");
    }
    final List<CigarElement> elems = read.getCigarElements();
    if (cigarStartIndex < 0 || cigarEndIndex > elems.size() || cigarStartIndex > cigarEndIndex) {
        throw new IllegalArgumentException("invalid index:" + 0 + " -" + elems.size());
    }
    int result = 0;
    for (int i = cigarStartIndex; i < cigarEndIndex; i++) {
        final CigarElement cigarElement = elems.get(i);
        switch(cigarElement.getOperator()) {
            case M:
            case D:
            case N:
            case EQ:
            case X:
            case S:
            case H:
                result += cigarElement.getLength();
                break;
            case I:
            case //for these two, nothing happens.
            P:
                break;
            default:
                throw new GATKException("Unsupported cigar operator: " + cigarElement.getOperator());
        }
    }
    return result;
}
Also used : GATKException(org.broadinstitute.hellbender.exceptions.GATKException) CigarElement(htsjdk.samtools.CigarElement)

Example 27 with CigarElement

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

the class CigarUtils method trimReadToUnclippedBases.

/**
     * Removes all clipping operators from the cigar.
     */
public static Cigar trimReadToUnclippedBases(final Cigar cigar) {
    Utils.nonNull(cigar, "cigar is null");
    final List<CigarElement> elements = new ArrayList<>(cigar.numCigarElements());
    for (final CigarElement ce : cigar.getCigarElements()) {
        if (!isClipOperator(ce.getOperator())) {
            elements.add(ce);
        }
    }
    return new Cigar(elements);
}
Also used : Cigar(htsjdk.samtools.Cigar) CigarElement(htsjdk.samtools.CigarElement)

Example 28 with CigarElement

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

the class CigarUtils method startsWithDeletionIgnoringClips.

/**
     * Checks if cigar starts with a deletion (ignoring any clips at the beginning).
     */
private static boolean startsWithDeletionIgnoringClips(final List<CigarElement> elems) {
    final Iterator<CigarElement> iter = elems.iterator();
    boolean isClip = true;
    CigarOperator op = null;
    while (iter.hasNext() && isClip) {
        //consume clips at the beginning
        final CigarElement elem = iter.next();
        op = elem.getOperator();
        isClip = (op == CigarOperator.HARD_CLIP || op == CigarOperator.SOFT_CLIP);
    }
    //once all clips are consumed, is it a deletion or not?
    return op == CigarOperator.DELETION;
}
Also used : CigarOperator(htsjdk.samtools.CigarOperator) CigarElement(htsjdk.samtools.CigarElement)

Example 29 with CigarElement

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

the class CigarUtils method leftAlignCigarSequentially.

/**
     * Left align the given cigar sequentially. This is needed because AlignmentUtils doesn't accept cigars with more than one indel in them.
     * This is a target of future work to incorporate and generalize into AlignmentUtils for use by others.
     * @param cigar     the cigar to left align
     * @param refSeq    the reference byte array
     * @param readSeq   the read byte array
     * @param refIndex  0-based alignment start position on ref
     * @param readIndex 0-based alignment start position on read
     * @return          the left-aligned cigar
     */
public static Cigar leftAlignCigarSequentially(final Cigar cigar, final byte[] refSeq, final byte[] readSeq, int refIndex, int readIndex) {
    Utils.nonNull(cigar, "cigar null");
    Utils.nonNull(refSeq, "refSeq null");
    Utils.nonNull(readSeq, "readSeq null");
    final Cigar cigarToReturn = new Cigar();
    Cigar cigarToAlign = new Cigar();
    for (int i = 0; i < cigar.numCigarElements(); i++) {
        final CigarElement ce = cigar.getCigarElement(i);
        if (ce.getOperator() == CigarOperator.D || ce.getOperator() == CigarOperator.I) {
            cigarToAlign.add(ce);
            final Cigar leftAligned = AlignmentUtils.leftAlignSingleIndel(cigarToAlign, refSeq, readSeq, refIndex, readIndex, false);
            for (final CigarElement toAdd : leftAligned.getCigarElements()) {
                cigarToReturn.add(toAdd);
            }
            refIndex += cigarToAlign.getReferenceLength();
            readIndex += cigarToAlign.getReadLength();
            cigarToAlign = new Cigar();
        } else {
            cigarToAlign.add(ce);
        }
    }
    if (!cigarToAlign.isEmpty()) {
        for (final CigarElement toAdd : cigarToAlign.getCigarElements()) {
            cigarToReturn.add(toAdd);
        }
    }
    final Cigar result = AlignmentUtils.consolidateCigar(cigarToReturn);
    Utils.validate(result.getReferenceLength() == cigar.getReferenceLength(), () -> "leftAlignCigarSequentially failed to produce a valid CIGAR.  Reference lengths differ.  Initial cigar " + cigar + " left aligned into " + result);
    return result;
}
Also used : Cigar(htsjdk.samtools.Cigar) CigarElement(htsjdk.samtools.CigarElement)

Example 30 with CigarElement

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

the class CigarUtils method hasConsecutiveIndels.

/**
     * Checks if cigar has consecutive I/D elements.
     */
private static boolean hasConsecutiveIndels(final List<CigarElement> elems) {
    boolean prevIndel = false;
    for (final CigarElement elem : elems) {
        final CigarOperator op = elem.getOperator();
        final boolean isIndel = (op == CigarOperator.INSERTION || op == CigarOperator.DELETION);
        if (prevIndel && isIndel) {
            return true;
        }
        prevIndel = isIndel;
    }
    return false;
}
Also used : CigarOperator(htsjdk.samtools.CigarOperator) CigarElement(htsjdk.samtools.CigarElement)

Aggregations

CigarElement (htsjdk.samtools.CigarElement)164 Cigar (htsjdk.samtools.Cigar)97 ArrayList (java.util.ArrayList)50 SAMRecord (htsjdk.samtools.SAMRecord)49 CigarOperator (htsjdk.samtools.CigarOperator)34 SAMRecordIterator (htsjdk.samtools.SAMRecordIterator)32 SamReader (htsjdk.samtools.SamReader)31 SAMFileHeader (htsjdk.samtools.SAMFileHeader)24 SAMSequenceDictionaryProgress (com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress)22 Test (org.testng.annotations.Test)19 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)17 BaseTest (org.broadinstitute.hellbender.utils.test.BaseTest)17 File (java.io.File)16 SAMFileWriter (htsjdk.samtools.SAMFileWriter)15 SAMReadGroupRecord (htsjdk.samtools.SAMReadGroupRecord)14 Interval (htsjdk.samtools.util.Interval)14 IOException (java.io.IOException)14 VisibleForTesting (com.google.common.annotations.VisibleForTesting)13 List (java.util.List)13 GATKException (org.broadinstitute.hellbender.exceptions.GATKException)13