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

Example 26 with CigarOperator

use of htsjdk.samtools.CigarOperator in project jvarkit by lindenb.

the class ReadClipper method clip.

public SAMRecord clip(final SAMRecord rec, final Interval fragment) {
    if (rec.getReadUnmappedFlag()) {
        return rec;
    }
    if (!fragment.getContig().equals(rec.getContig())) {
        return rec;
    }
    if (rec.getAlignmentEnd() < fragment.getStart()) {
        return rec;
    }
    if (rec.getAlignmentStart() > fragment.getEnd()) {
        return rec;
    }
    Cigar cigar = rec.getCigar();
    if (cigar == null) {
        LOG.warning("cigar missing in " + rec);
        return rec;
    }
    final List<BaseOp> bases = new ArrayList<>(cigar.getCigarElements().stream().mapToInt(C -> C.getLength()).sum());
    // expand cigar
    int refPos = rec.getUnclippedStart();
    for (final CigarElement ce : cigar.getCigarElements()) {
        final CigarOperator op = ce.getOperator();
        if (op.equals(CigarOperator.P))
            continue;
        for (int x = 0; x < ce.getLength(); ++x) {
            bases.add(new BaseOp(op, op.consumesReferenceBases() || op.isClipping() ? refPos : -1));
            if (op.consumesReferenceBases() || op.isClipping()) {
                refPos++;
            }
        }
    }
    /* 5' side */
    int newStart = rec.getAlignmentStart();
    int x = 0;
    while (x < bases.size()) {
        final BaseOp b = bases.get(x);
        if (b.refPos != -1 && b.isMatch() && b.refPos >= fragment.getStart()) {
            newStart = b.refPos;
            break;
        } else if (b.isDeletion()) {
            bases.remove(x);
            continue;
        } else if (!b.op.isClipping()) {
            b.op = CigarOperator.S;
            ++x;
        } else {
            ++x;
        }
    }
    /* 3' side */
    x = bases.size() - 1;
    while (x >= 0) {
        final BaseOp b = bases.get(x);
        if (b.refPos != -1 && b.isMatch() && b.refPos <= fragment.getEnd()) {
            break;
        } else if (b.isDeletion()) {
            bases.remove(x);
            --x;
            continue;
        } else if (!b.op.isClipping()) {
            b.op = CigarOperator.S;
            --x;
        } else {
            --x;
        }
    }
    // build new cigar
    boolean found_M = false;
    final List<CigarElement> newcigarlist = new ArrayList<>();
    x = 0;
    while (x < bases.size()) {
        final CigarOperator op = bases.get(x).op;
        if (op.equals(CigarOperator.M) || op.equals(CigarOperator.EQ)) {
            found_M = true;
        }
        int x2 = x;
        while (x2 < bases.size() && op.equals(bases.get(x2).op)) {
            ++x2;
        }
        newcigarlist.add(new CigarElement(x2 - x, op));
        x = x2;
    }
    if (!found_M) {
        SAMUtils.makeReadUnmappedWithOriginalTags(rec);
        if (this.programGroup != null) {
            rec.setAttribute("PG", programGroup);
        }
        return rec;
    }
    cigar = new Cigar(newcigarlist);
    rec.setCigar(cigar);
    rec.setAlignmentStart(newStart);
    if (this.programGroup != null) {
        rec.setAttribute("PG", programGroup);
    }
    return rec;
}

Example 27 with CigarOperator

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

the class PileupElementUnitTest method testIsImmediatelyAfter_afterFirstD.

@Test
public void testIsImmediatelyAfter_afterFirstD() {
    final GATKRead read = ArtificialReadUtils.createArtificialRead("10M10D10M");
    final PileupElement pe = PileupElement.createPileupForReadAndOffset(read, 10);
    Assert.assertFalse(pe.atEndOfCurrentCigar(), "atEndOfCurrentCigar");
    Assert.assertTrue(pe.atStartOfCurrentCigar(), "atStartOfCurrentCigar");
    Assert.assertTrue(pe.isImmediatelyAfter(D));
    for (final CigarOperator op : CigarOperator.values()) {
        if (op != D) {
            Assert.assertFalse(pe.isImmediatelyAfter(op), "isImmediatelyAfter " + op);
        }
    }
    for (final CigarOperator op : CigarOperator.values()) {
        Assert.assertFalse(pe.isImmediatelyBefore(op));
    }
}

Example 28 with CigarOperator

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

the class PileupElementUnitTest method testImmediateBeforeAndAfterTest_simple.

@Test(dataProvider = "PrevAndNextTest_simple")
public void testImmediateBeforeAndAfterTest_simple(final GATKRead read, final CigarOperator middleOp) {
    final AlignmentStateMachine state = new AlignmentStateMachine(read);
    state.stepForwardOnGenome();
    //before the 'middleOp'
    final PileupElement pe1 = state.makePileupElement();
    Assert.assertEquals(pe1.getAdjacentOperator(PileupElement.Direction.PREV), null, "PREV");
    Assert.assertEquals(pe1.getAdjacentOperator(PileupElement.Direction.NEXT), middleOp, "NEXT");
    for (final CigarOperator op : CigarOperator.values()) {
        Assert.assertEquals(pe1.isImmediatelyBefore(op), middleOp == op, op.toString());
        Assert.assertFalse(pe1.isImmediatelyAfter(op), op.toString());
    }
    state.stepForwardOnGenome();
    //after the 'middleOp'
    final PileupElement pe2 = state.makePileupElement();
    Assert.assertEquals(pe2.getAdjacentOperator(PileupElement.Direction.PREV), middleOp, "PREV");
    Assert.assertEquals(pe2.getAdjacentOperator(PileupElement.Direction.NEXT), null, "NEXT");
    for (final CigarOperator op : CigarOperator.values()) {
        Assert.assertFalse(pe2.isImmediatelyBefore(op), op.toString());
        Assert.assertEquals(pe2.isImmediatelyAfter(op), op == middleOp, op.toString());
    }
}

Example 29 with CigarOperator

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

the class PileupElementUnitTest method makePrevAndNextTest.

@DataProvider(name = "PrevAndNextTest")
public Object[][] makePrevAndNextTest() {
    final List<Object[]> tests = new LinkedList<>();
    final List<CigarOperator> operators = Arrays.asList(CigarOperator.I, CigarOperator.P, CigarOperator.S);
    for (final CigarOperator firstOp : Arrays.asList(M)) {
        for (final CigarOperator lastOp : Arrays.asList(M, D)) {
            for (final int nIntermediate : Arrays.asList(1, 2, 3)) {
                for (final List<CigarOperator> combination : Utils.makePermutations(operators, nIntermediate, false)) {
                    final int readLength = 2 + combination.size();
                    final GATKRead read = ArtificialReadUtils.createArtificialRead(header, "read", 0, 1, readLength);
                    read.setBases(Utils.dupBytes((byte) 'A', readLength));
                    read.setBaseQualities(Utils.dupBytes((byte) 30, readLength));
                    String cigar = "1" + firstOp;
                    for (final CigarOperator op : combination) {
                        cigar += "1" + op;
                    }
                    cigar += "1" + lastOp;
                    read.setCigar(cigar);
                    tests.add(new Object[] { read, firstOp, lastOp, combination });
                }
            }
        }
    }
    return tests.toArray(new Object[][] {});
}

Example 30 with CigarOperator

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

the class PileupElementUnitTest method testIsImmediatelyAfter_insideMs.

@Test
public void testIsImmediatelyAfter_insideMs() {
    final GATKRead read = ArtificialReadUtils.createArtificialRead("10M10D10M");
    final PileupElement pe = PileupElement.createPileupForReadAndOffset(read, 5);
    Assert.assertFalse(pe.atStartOfCurrentCigar(), "atStartOfCurrentCigar");
    Assert.assertFalse(pe.atEndOfCurrentCigar(), "atEndOfCurrentCigar");
    for (final CigarOperator op : CigarOperator.values()) {
        Assert.assertFalse(pe.isImmediatelyAfter(op), "isImmediatelyAfter " + op);
    }
    for (final CigarOperator op : CigarOperator.values()) {
        Assert.assertFalse(pe.isImmediatelyBefore(op), "isImmediatelyBefore " + op);
    }
}