Search in sources :

Example 51 with CigarElement

use of htsjdk.samtools.CigarElement in project ASCIIGenome by dariober.

the class TextReadTest method getNumAlignedBasesFromCigar.

@Test
public void getNumAlignedBasesFromCigar() {
    // A reminder of how cigar works.
    SAMRecord rec = new SAMRecord(null);
    rec.setCigarString("2M2D4M");
    int noRefBases = rec.getCigar().getReferenceLength();
    assertEquals(8, noRefBases);
    // N same as D
    rec.setCigarString("2M2N4M");
    noRefBases = rec.getCigar().getReferenceLength();
    assertEquals(8, noRefBases);
    rec.setCigarString("2M2I4M");
    noRefBases = rec.getCigar().getReferenceLength();
    assertEquals(6, noRefBases);
    rec.setCigarString("2S10M4H");
    noRefBases = rec.getCigar().getReferenceLength();
    assertEquals(10, noRefBases);
    // Cigar oprators consuming:
    // Hard clipping doesn't consume bases on read or ref
    rec.setCigarString("2H");
    CigarElement el = rec.getCigar().getCigarElement(0);
    assertFalse(el.getOperator().consumesReadBases());
    assertFalse(el.getOperator().consumesReferenceBases());
    // Clips the read not the ref.
    rec.setCigarString("2S");
    el = rec.getCigar().getCigarElement(0);
    assertTrue(el.getOperator().consumesReadBases());
    assertFalse(el.getOperator().consumesReferenceBases());
    // This makes a gap in the read when aligned to ref
    rec.setCigarString("2D");
    el = rec.getCigar().getCigarElement(0);
    assertFalse(el.getOperator().consumesReadBases());
    assertTrue(el.getOperator().consumesReferenceBases());
    // Opposite of D
    rec.setCigarString("2I");
    el = rec.getCigar().getCigarElement(0);
    assertTrue(el.getOperator().consumesReadBases());
    assertFalse(el.getOperator().consumesReferenceBases());
}
Also used : SAMRecord(htsjdk.samtools.SAMRecord) CigarElement(htsjdk.samtools.CigarElement) Test(org.junit.Test)

Example 52 with CigarElement

use of htsjdk.samtools.CigarElement in project ASCIIGenome by dariober.

the class TextRead method getSoftUnclippedAlignmentStart.

public int getSoftUnclippedAlignmentStart(SAMRecord rec) {
    List<CigarElement> cigar = rec.getCigar().getCigarElements();
    if (cigar.size() == 0) {
        return rec.getUnclippedStart();
    }
    if (cigar.size() == 1 && cigar.get(0).getOperator().equals(CigarOperator.SOFT_CLIP)) {
        return rec.getUnclippedStart();
    }
    int offset = 0;
    for (int i = 0; i < cigar.size(); i++) {
        CigarElement op = cigar.get(i);
        if (op.getOperator().equals(CigarOperator.HARD_CLIP)) {
            continue;
        } else if (op.getOperator().equals(CigarOperator.SOFT_CLIP)) {
            offset += op.getLength();
        } else {
            break;
        }
    }
    int start = rec.getAlignmentStart() - offset;
    // }
    return start;
}
Also used : CigarElement(htsjdk.samtools.CigarElement)

Example 53 with CigarElement

use of htsjdk.samtools.CigarElement in project ASCIIGenome by dariober.

the class TextRead method setTextPositionsOfSkippedBases.

/**
 *List of positions on screen where the skipped bases (cigar op: N) start
 * @throws IOException
 * @throws InvalidGenomicCoordsException
 */
private void setTextPositionsOfSkippedBases() throws InvalidGenomicCoordsException, IOException {
    int genomicPosition = this.getSamRecord().getAlignmentStart();
    List<CigarElement> cigar = this.getSamRecord().getCigar().getCigarElements();
    for (CigarElement el : cigar) {
        if (el.getOperator().equals(CigarOperator.SKIPPED_REGION)) {
            int[] textPositions = new int[2];
            // +1 because textPosition is 1-based
            textPositions[0] = Utils.getIndexOfclosestValue(genomicPosition, this.gc.getMapping()) + 1;
            textPositions[1] = Utils.getIndexOfclosestValue(genomicPosition + el.getLength(), this.gc.getMapping()) + 1;
            this.textPositionsOfSkippedBases.add(textPositions);
        }
        ;
        if (el.getOperator().consumesReferenceBases()) {
            genomicPosition += el.getLength();
        }
    }
}
Also used : CigarElement(htsjdk.samtools.CigarElement)

Example 54 with CigarElement

use of htsjdk.samtools.CigarElement in project ASCIIGenome by dariober.

the class Track method isSNVRead.

/**
 *Return true if samrecord contains a mismatch or insertion/deletion in the target region.
 */
private boolean isSNVRead(SAMRecord rec, boolean variantOnly) {
    boolean passed = false;
    int varFrom = this.getFeatureFilter().getVariantFrom();
    int varTo = this.getFeatureFilter().getVariantTo();
    if (this.getFeatureFilter().getVariantChrom().equals(rec.getReferenceName()) && varFrom <= rec.getAlignmentEnd() && rec.getAlignmentStart() <= varTo) {
        // Variant read filter is set and this read overlaps it.
        if (!variantOnly) {
            // No need to check whether read is variant.
            return true;
        }
        int readPos = 0;
        int refPos = rec.getAlignmentStart();
        for (CigarElement cigar : rec.getCigar().getCigarElements()) {
            if (cigar.getOperator().equals(CigarOperator.SOFT_CLIP)) {
                readPos += cigar.getLength();
            } else if (cigar.getOperator().equals(CigarOperator.M) || cigar.getOperator().equals(CigarOperator.EQ) || cigar.getOperator().equals(CigarOperator.X)) {
                for (int i = 0; i < cigar.getLength(); i++) {
                    if (refPos >= varFrom && refPos <= varTo && rec.getReadLength() > 0) {
                        byte readBase = rec.getReadBases()[readPos];
                        byte refBase = this.getFeatureFilter().getFaSeq()[refPos - varFrom];
                        if (readBase != refBase) {
                            passed = true;
                            break;
                        }
                    }
                    readPos++;
                    refPos++;
                }
            } else if (cigar.getOperator().equals(CigarOperator.DELETION)) {
                // ^^^^
                for (int i = 0; i < cigar.getLength(); i++) {
                    if (refPos >= varFrom && refPos <= varTo) {
                        passed = true;
                        break;
                    }
                    refPos++;
                }
            } else if (cigar.getOperator().equals(CigarOperator.INSERTION)) {
                // ^
                for (int i = 0; i < cigar.getLength(); i++) {
                    if (refPos >= varFrom && refPos <= varTo) {
                        passed = true;
                        break;
                    }
                    readPos++;
                }
            } else if (cigar.getOperator().equals(CigarOperator.HARD_CLIP)) {
            // 
            } else if (cigar.getOperator().equals(CigarOperator.SKIPPED_REGION)) {
                // Same deletion but it's not a mismatch
                refPos += cigar.getLength();
            } else if (cigar.getOperator().equals(CigarOperator.PADDING)) {
            // Not sure what to do with this...
            }
            if (passed) {
                break;
            }
        }
    } else {
        // Variant read filter is set and this read does not overlap it.
        passed = false;
    }
    return passed;
}
Also used : CigarElement(htsjdk.samtools.CigarElement)

Example 55 with CigarElement

use of htsjdk.samtools.CigarElement in project gridss by PapenfussLab.

the class CigarUtil method convertOpenEndedOperations.

/**
 * Converts indel operations that do not contain anchoring alignments on both sides to
 * a soft cliped version of the operator
 * @param list
 */
private static void convertOpenEndedOperations(List<CigarElement> list) {
    for (int i = list.size() - 1; i >= 0; i--) {
        CigarElement ce = list.get(i);
        if (ce.getOperator() == CigarOperator.SOFT_CLIP || ce.getOperator() == CigarOperator.HARD_CLIP) {
        } else if (ce.getOperator() == CigarOperator.INSERTION) {
            list.set(i, new CigarElement(ce.getLength(), CigarOperator.SOFT_CLIP));
        } else if (ce.getOperator() == CigarOperator.DELETION) {
            list.remove(i);
        } else {
            break;
        }
    }
    for (int i = 0; i < list.size(); i++) {
        CigarElement ce = list.get(i);
        if (ce.getOperator() == CigarOperator.SOFT_CLIP || ce.getOperator() == CigarOperator.HARD_CLIP) {
        } else if (ce.getOperator() == CigarOperator.INSERTION) {
            list.set(i, new CigarElement(ce.getLength(), CigarOperator.SOFT_CLIP));
        } else if (ce.getOperator() == CigarOperator.DELETION) {
            list.remove(i);
            i--;
        } else {
            break;
        }
    }
}
Also used : 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