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

Example 6 with CigarElement

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

the class Civar method toCigar.

/**
     * Transforms a civar into the equivalent Cigar.
     * @return never {@code null}.
     */
public Cigar toCigar(final int templateLength) {
    int minSeqLen = minimumTemplateSequenceSize();
    Utils.validateArg(expands() || templateLength == minSeqLen, "the sequence provided does not match this Civar size " + templateLength + " != " + minSeqLen);
    Utils.validateArg(templateLength >= minSeqLen, "the sequence provided is too small for this Civar " + templateLength + " < " + minSeqLen);
    int starCount = starCount();
    int paddingTotal = templateLength - minSeqLen;
    int starPadding = starCount == 0 ? 0 : paddingTotal / starCount;
    int excessPadding = starCount == 0 ? paddingTotal : paddingTotal % starCount;
    // We first get the equivalent cigar elements for the elements in the Civar.
    final List<CigarElement> cigarElements = new LinkedList<>();
    for (final Element e : elements()) {
        final int size = e.size(starPadding, excessPadding);
        excessPadding -= e.excessPaddingUsed(excessPadding);
        switch(e.operator()) {
            case EMBEDDED:
                cigarElements.addAll(e.embedded.toCigar(size).getCigarElements());
                break;
            case MATCH:
            case TRANSITION:
            case COMPLEMENT:
            case TRANSVERSION:
                cigarElements.add(new CigarElement(size, CigarOperator.M));
                break;
            case INSERTION:
                cigarElements.add(new CigarElement(size, CigarOperator.I));
                break;
            case DELETION:
                cigarElements.add(new CigarElement(size, CigarOperator.D));
                break;
            default:
        }
    }
    // No we look for consequitive elements with the same operator and we merge them.
    final ListIterator<CigarElement> it = cigarElements.listIterator();
    while (it.hasNext()) {
        final CigarElement thisElement = it.next();
        if (!it.hasNext())
            continue;
        final CigarElement nextElement = it.next();
        if (thisElement.getOperator() == nextElement.getOperator()) {
            final int nextLength = nextElement.getLength();
            it.remove();
            it.previous();
            it.set(new CigarElement(thisElement.getLength() + nextLength, thisElement.getOperator()));
        } else
            it.previous();
    }
    return new Cigar(cigarElements);
}

Example 7 with CigarElement

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

the class PileupElement method getBetween.

/**
     * Helper function to get cigar elements between this and either the prev or next genomic position.
     * Note that this method stops accumulating elements at the first "on-genome" operator it encounters.
     *
     * @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 List<CigarElement> getBetween(final Direction direction) {
    final int increment = direction.getIncrement();
    final List<CigarElement> elements = new ArrayList<>();
    final List<CigarElement> cigarElements = read.getCigarElements();
    final int nCigarElements = cigarElements.size();
    for (int i = currentCigarOffset + increment; i >= 0 && i < nCigarElements; i += increment) {
        final CigarElement elt = cigarElements.get(i);
        if (ON_GENOME_OPERATORS.contains(elt.getOperator())) {
            break;
        } else {
            elements.add(elt);
        }
    }
    if (increment < 0) {
        Collections.reverse(elements);
    }
    return elements;
}

Example 8 with CigarElement

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

the class Haplotype method getConsolidatedPaddedCigar.

/**
     * Get the haplotype cigar extended by padSize M at the tail, consolidated into a clean cigar
     *
     * @param padSize how many additional Ms should be appended to the end of this cigar.  Must be >= 0
     * @return a newly allocated Cigar that consolidate(getCigar + padSize + M)
     */
public Cigar getConsolidatedPaddedCigar(final int padSize) {
    Utils.validateArg(padSize >= 0, () -> "padSize must be >= 0 but got " + padSize);
    final Cigar extendedHaplotypeCigar = new Cigar(getCigar().getCigarElements());
    if (padSize > 0) {
        extendedHaplotypeCigar.add(new CigarElement(padSize, CigarOperator.M));
    }
    return AlignmentUtils.consolidateCigar(extendedHaplotypeCigar);
}

Example 9 with CigarElement

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

the class AlignmentUtils method leftAlignSingleIndel.

/**
     * See the documentation for AlignmentUtils.leftAlignIndel() for more details.
     *
     * This flavor of the left alignment works if and only if the alignment has one - and only one - indel.
     * An exception is thrown if there are no indels or more than 1 indel in the alignment.
     *
     * @param cigar     structure of the original alignment -- cannot be null
     * @param refSeq    reference sequence the read is aligned to
     * @param readSeq   read sequence
     * @param refIndex  0-based alignment start position on ref
     * @param readIndex 0-based alignment start position on read
     * @param cleanupCigar if true, we'll cleanup the resulting cigar element, removing 0 length elements and deletions from the first cigar position
     * @return a non-null cigar, in which the single indel is guaranteed to be placed at the leftmost possible position across a repeat (if any)
     */
public static Cigar leftAlignSingleIndel(Cigar cigar, final byte[] refSeq, final byte[] readSeq, final int refIndex, final int readIndex, final boolean cleanupCigar) {
    ensureLeftAlignmentHasGoodArguments(cigar, refSeq, readSeq, refIndex, readIndex);
    int indexOfIndel = -1;
    for (int i = 0; i < cigar.numCigarElements(); i++) {
        CigarElement ce = cigar.getCigarElement(i);
        if (ce.getOperator() == CigarOperator.D || ce.getOperator() == CigarOperator.I) {
            // if there is more than 1 indel, exception out
            if (indexOfIndel != -1)
                throw new IllegalArgumentException("attempting to left align a CIGAR that has more than 1 indel in its alignment");
            indexOfIndel = i;
        }
    }
    // if there is no indel, exception out
    if (indexOfIndel == -1)
        throw new IllegalArgumentException("attempting to left align a CIGAR that has no indels in its alignment");
    // if the alignment starts with an insertion (so that there is no place on the read to move that insertion further left), we are done
    if (indexOfIndel == 0)
        return cigar;
    final int indelLength = cigar.getCigarElement(indexOfIndel).getLength();
    byte[] altString = createIndelString(cigar, indexOfIndel, refSeq, readSeq, refIndex, readIndex);
    if (altString == null)
        return cigar;
    Cigar newCigar = cigar;
    for (int i = 0; i < indelLength; i++) {
        newCigar = moveCigarLeft(newCigar, indexOfIndel);
        byte[] newAltString = createIndelString(newCigar, indexOfIndel, refSeq, readSeq, refIndex, readIndex);
        // check to make sure we haven't run off the end of the read
        boolean reachedEndOfRead = cigarHasZeroSizeElement(newCigar);
        if (Arrays.equals(altString, newAltString)) {
            cigar = newCigar;
            i = -1;
            if (reachedEndOfRead)
                cigar = cleanupCigar ? cleanUpCigar(cigar) : cigar;
        }
        if (reachedEndOfRead)
            break;
    }
    return cigar;
}

Example 10 with CigarElement

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

the class AlignmentUtils method calcAlignmentByteArrayOffset.

/**
     * Calculate the index into the read's bases of the beginning of the encompassing cigar element for a given cigar and offset
     *
     * @param cigar            the read's CIGAR -- cannot be null
     * @param offset           the offset to use for the calculation or -1 if in the middle of a deletion
     * @param isDeletion       are we in the middle of a deletion?
     * @param alignmentStart   the alignment start of the read
     * @param refLocus         the reference position of the offset
     * @return a non-negative int index
     */
public static int calcAlignmentByteArrayOffset(final Cigar cigar, final int offset, final boolean isDeletion, final int alignmentStart, final int refLocus) {
    if (cigar == null)
        throw new IllegalArgumentException("attempting to find the alignment position from a CIGAR that is null");
    if (offset < -1)
        throw new IllegalArgumentException("attempting to find the alignment position with an offset that is negative (and not -1)");
    if (alignmentStart < 0)
        throw new IllegalArgumentException("attempting to find the alignment position from an alignment start that is negative");
    if (refLocus < 0)
        throw new IllegalArgumentException("attempting to find the alignment position from a reference position that is negative");
    if (offset >= cigar.getReadLength())
        throw new IllegalArgumentException("attempting to find the alignment position of an offset than is larger than the read length");
    int pileupOffset = offset;
    // Reassign the offset if we are in the middle of a deletion because of the modified representation of the read bases
    if (isDeletion) {
        pileupOffset = refLocus - alignmentStart;
        final CigarElement ce = cigar.getCigarElement(0);
        if (ce.getOperator() == CigarOperator.S) {
            pileupOffset += ce.getLength();
        }
    }
    int pos = 0;
    int alignmentPos = 0;
    for (int iii = 0; iii < cigar.numCigarElements(); iii++) {
        final CigarElement ce = cigar.getCigarElement(iii);
        final int elementLength = ce.getLength();
        switch(ce.getOperator()) {
            case I:
            case // TODO -- I don't think that soft clips should be treated the same as inserted bases here. Investigation needed.
            S:
                pos += elementLength;
                if (pos >= pileupOffset) {
                    return alignmentPos;
                }
                break;
            case D:
                if (!isDeletion) {
                    alignmentPos += elementLength;
                } else {
                    if (pos + elementLength - 1 >= pileupOffset) {
                        return alignmentPos + (pileupOffset - pos);
                    } else {
                        pos += elementLength;
                        alignmentPos += elementLength;
                    }
                }
                break;
            case M:
            case EQ:
            case X:
                if (pos + elementLength - 1 >= pileupOffset) {
                    return alignmentPos + (pileupOffset - pos);
                } else {
                    pos += elementLength;
                    alignmentPos += elementLength;
                }
                break;
            case H:
            case P:
            case N:
                break;
            default:
                throw new GATKException("Unsupported cigar operator: " + ce.getOperator());
        }
    }
    return alignmentPos;
}