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

Example 31 with CigarElement

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

the class LIBS_position method stepForwardOnGenome.

/**
     * Steps forward on the genome.  Returns false when done reading the read, true otherwise.
     */
@SuppressWarnings("fallthrough")
public boolean stepForwardOnGenome() {
    if (currentOperatorIndex == numOperators)
        return false;
    CigarElement curElement = read.getCigar().getCigarElement(currentOperatorIndex);
    if (currentPositionOnOperator >= curElement.getLength()) {
        if (++currentOperatorIndex == numOperators)
            return false;
        curElement = read.getCigar().getCigarElement(currentOperatorIndex);
        currentPositionOnOperator = 0;
    }
    switch(curElement.getOperator()) {
        // insertion w.r.t. the reference
        case I:
        //                    break;
        case // soft clip
        S:
            currentReadOffset += curElement.getLength();
        // hard clip
        case H:
        case // padding
        P:
            currentOperatorIndex++;
            return stepForwardOnGenome();
        // deletion w.r.t. the reference
        case D:
        case // reference skip (looks and gets processed just like a "deletion", just different logical meaning)
        N:
            currentPositionOnOperator++;
            currentGenomeOffset++;
            break;
        case M:
        case EQ:
        case X:
            sawMop = true;
            currentReadOffset++;
            currentPositionOnOperator++;
            currentGenomeOffset++;
            break;
        default:
            throw new IllegalStateException("No support for cigar op: " + curElement.getOperator());
    }
    final boolean isFirstOp = currentOperatorIndex == 0;
    final boolean isLastOp = currentOperatorIndex == numOperators - 1;
    final boolean isFirstBaseOfOp = currentPositionOnOperator == 1;
    final boolean isLastBaseOfOp = currentPositionOnOperator == curElement.getLength();
    isBeforeDeletionStart = isBeforeOp(read.getCigar(), currentOperatorIndex, CigarOperator.D, isLastOp, isLastBaseOfOp);
    isBeforeDeletedBase = isBeforeDeletionStart || (!isLastBaseOfOp && curElement.getOperator() == CigarOperator.D);
    isAfterDeletionEnd = isAfterOp(read.getCigar(), currentOperatorIndex, CigarOperator.D, isFirstOp, isFirstBaseOfOp);
    isAfterDeletedBase = isAfterDeletionEnd || (!isFirstBaseOfOp && curElement.getOperator() == CigarOperator.D);
    isBeforeInsertion = isBeforeOp(read.getCigar(), currentOperatorIndex, CigarOperator.I, isLastOp, isLastBaseOfOp) || (!sawMop && curElement.getOperator() == CigarOperator.I);
    isAfterInsertion = isAfterOp(read.getCigar(), currentOperatorIndex, CigarOperator.I, isFirstOp, isFirstBaseOfOp);
    isNextToSoftClip = isBeforeOp(read.getCigar(), currentOperatorIndex, CigarOperator.S, isLastOp, isLastBaseOfOp) || isAfterOp(read.getCigar(), currentOperatorIndex, CigarOperator.S, isFirstOp, isFirstBaseOfOp);
    return true;
}

Example 32 with CigarElement

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

the class LocusIteratorByStateBaseTest method makePermutationTest.

private LIBSTest makePermutationTest(final List<CigarElement> elements) {
    CigarElement last = null;
    boolean hasMatch = false;
    // starts with D => bad
    if (startsWithDeletion(elements))
        return null;
    // ends with D => bad
    if (elements.get(elements.size() - 1).getOperator() == CigarOperator.D)
        return null;
    // make sure it's valid
    String cigar = "";
    for (final CigarElement ce : elements) {
        if (ce.getOperator() == CigarOperator.N)
            // TODO -- don't support N
            return null;
        // abort on a bad cigar
        if (last != null) {
            if (ce.getOperator() == last.getOperator())
                return null;
            if (isIndel(ce) && isIndel(last))
                return null;
        }
        cigar += ce.getLength() + ce.getOperator().toString();
        last = ce;
        hasMatch = hasMatch || ce.getOperator() == CigarOperator.M;
    }
    if (!hasMatch && elements.size() == 1 && !(last.getOperator() == CigarOperator.I || last.getOperator() == CigarOperator.S))
        return null;
    return new LIBSTest(cigar);
}

Example 33 with CigarElement

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

the class ReferenceConfidenceModel method createReferenceHaplotype.

/**
     * Create a reference haplotype for an active region
     *
     * @param activeRegion the active region
     * @param refBases the ref bases
     * @param paddedReferenceLoc the location spanning of the refBases -- can be longer than activeRegion.getLocation()
     * @return a reference haplotype
     */
public static Haplotype createReferenceHaplotype(final AssemblyRegion activeRegion, final byte[] refBases, final SimpleInterval paddedReferenceLoc) {
    Utils.nonNull(activeRegion, "null region");
    Utils.nonNull(refBases, "null refBases");
    Utils.nonNull(paddedReferenceLoc, "null paddedReferenceLoc");
    final int alignmentStart = activeRegion.getExtendedSpan().getStart() - paddedReferenceLoc.getStart();
    if (alignmentStart < 0) {
        throw new IllegalStateException("Bad alignment start in createReferenceHaplotype " + alignmentStart);
    }
    final Haplotype refHaplotype = new Haplotype(refBases, true);
    refHaplotype.setAlignmentStartHapwrtRef(alignmentStart);
    final Cigar c = new Cigar();
    c.add(new CigarElement(refHaplotype.getBases().length, CigarOperator.M));
    refHaplotype.setCigar(c);
    return refHaplotype;
}

Example 34 with CigarElement

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

the class ReadThreadingGraph method mergeDanglingHead.

/**
     * Actually merge the dangling head if possible
     *
     * @param danglingHeadMergeResult   the result from generating a Cigar for the dangling head against the reference
     * @return 1 if merge was successful, 0 otherwise
     */
@VisibleForTesting
int mergeDanglingHead(final DanglingChainMergeHelper danglingHeadMergeResult) {
    final List<CigarElement> elements = danglingHeadMergeResult.cigar.getCigarElements();
    final CigarElement firstElement = elements.get(0);
    Utils.validateArg(firstElement.getOperator() == CigarOperator.M, "The first Cigar element must be an M");
    final int indexesToMerge = bestPrefixMatch(danglingHeadMergeResult.referencePathString, danglingHeadMergeResult.danglingPathString, firstElement.getLength());
    if (indexesToMerge <= 0) {
        return 0;
    }
    // we can't push back the reference path
    if (indexesToMerge >= danglingHeadMergeResult.referencePath.size() - 1) {
        return 0;
    }
    // but we can manipulate the dangling path if we need to
    if (indexesToMerge >= danglingHeadMergeResult.danglingPath.size() && !extendDanglingPathAgainstReference(danglingHeadMergeResult, indexesToMerge - danglingHeadMergeResult.danglingPath.size() + 2)) {
        return 0;
    }
    addEdge(danglingHeadMergeResult.referencePath.get(indexesToMerge + 1), danglingHeadMergeResult.danglingPath.get(indexesToMerge), ((MyEdgeFactory) getEdgeFactory()).createEdge(false, 1));
    return 1;
}

Example 35 with CigarElement

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

the class ReadThreadingGraph method mergeDanglingTail.

/**
     * Actually merge the dangling tail if possible
     *
     * @param danglingTailMergeResult   the result from generating a Cigar for the dangling tail against the reference
     * @return 1 if merge was successful, 0 otherwise
     */
@VisibleForTesting
final int mergeDanglingTail(final DanglingChainMergeHelper danglingTailMergeResult) {
    final List<CigarElement> elements = danglingTailMergeResult.cigar.getCigarElements();
    final CigarElement lastElement = elements.get(elements.size() - 1);
    Utils.validateArg(lastElement.getOperator() == CigarOperator.M, "The last Cigar element must be an M");
    final int lastRefIndex = danglingTailMergeResult.cigar.getReferenceLength() - 1;
    final int matchingSuffix = Math.min(longestSuffixMatch(danglingTailMergeResult.referencePathString, danglingTailMergeResult.danglingPathString, lastRefIndex), lastElement.getLength());
    if (matchingSuffix == 0) {
        return 0;
    }
    final int altIndexToMerge = Math.max(danglingTailMergeResult.cigar.getReadLength() - matchingSuffix - 1, 0);
    // there is an important edge condition that we need to handle here: Smith-Waterman correctly calculates that there is a
    // deletion, that deletion is left-aligned such that the LCA node is part of that deletion, and the rest of the dangling
    // tail is a perfect match to the suffix of the reference path.  In this case we need to push the reference index to merge
    // down one position so that we don't incorrectly cut a base off of the deletion.
    final boolean firstElementIsDeletion = elements.get(0).getOperator() == CigarOperator.D;
    final boolean mustHandleLeadingDeletionCase = firstElementIsDeletion && (elements.get(0).getLength() + matchingSuffix == lastRefIndex + 1);
    final int refIndexToMerge = lastRefIndex - matchingSuffix + 1 + (mustHandleLeadingDeletionCase ? 1 : 0);
    // merge back to the LCA, which results in a cycle in the graph.  So we do not want to merge in such a case.
    if (refIndexToMerge == 0) {
        return 0;
    }
    // it's safe to merge now
    addEdge(danglingTailMergeResult.danglingPath.get(altIndexToMerge), danglingTailMergeResult.referencePath.get(refIndexToMerge), ((MyEdgeFactory) getEdgeFactory()).createEdge(false, 1));
    return 1;
}