Search in sources :

Example 96 with CigarElement

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

the class SWPairwiseAlignment method printAlignment.

public void printAlignment(final byte[] ref, final byte[] read, final int width) {
    final StringBuilder bread = new StringBuilder();
    final StringBuilder bref = new StringBuilder();
    final StringBuilder match = new StringBuilder();
    int i = 0;
    int j = 0;
    final int offset = getAlignmentStart2wrt1();
    Cigar cigar = getCigar();
    if (overhangStrategy != OverhangStrategy.SOFTCLIP) {
        // otherwise offset is never negative
        if (offset < 0) {
            for (; j < (-offset); j++) {
                bread.append((char) read[j]);
                bref.append(' ');
                match.append(' ');
            }
            // at negative offsets, our cigar's first element carries overhanging bases
            // that we have just printed above. Tweak the first element to
            // exclude those bases. Here we create a new list of cigar elements, so the original
            // list/original cigar are unchanged (they are unmodifiable anyway!)
            final List<CigarElement> tweaked = new ArrayList<>();
            tweaked.addAll(cigar.getCigarElements());
            tweaked.set(0, new CigarElement(cigar.getCigarElement(0).getLength() + offset, cigar.getCigarElement(0).getOperator()));
            cigar = new Cigar(tweaked);
        }
    }
    if (offset > 0) {
        // note: the way this implementation works, cigar will ever start from S *only* if read starts before the ref, i.e. offset = 0
        for (; i < getAlignmentStart2wrt1(); i++) {
            bref.append((char) ref[i]);
            bread.append(' ');
            match.append(' ');
        }
    }
    for (final CigarElement e : cigar.getCigarElements()) {
        switch(e.getOperator()) {
            case M:
                for (int z = 0; z < e.getLength(); z++, i++, j++) {
                    bref.append((i < ref.length) ? (char) ref[i] : ' ');
                    bread.append((j < read.length) ? (char) read[j] : ' ');
                    match.append((i < ref.length && j < read.length) ? (ref[i] == read[j] ? '.' : '*') : ' ');
                }
                break;
            case I:
                for (int z = 0; z < e.getLength(); z++, j++) {
                    bref.append('-');
                    bread.append((char) read[j]);
                    match.append('I');
                }
                break;
            case S:
                for (int z = 0; z < e.getLength(); z++, j++) {
                    bref.append(' ');
                    bread.append((char) read[j]);
                    match.append('S');
                }
                break;
            case D:
                for (int z = 0; z < e.getLength(); z++, i++) {
                    bref.append((char) ref[i]);
                    bread.append('-');
                    match.append('D');
                }
                break;
            default:
                throw new GATKException("Unexpected Cigar element:" + e.getOperator());
        }
    }
    for (; i < ref.length; i++) bref.append((char) ref[i]);
    for (; j < read.length; j++) bread.append((char) read[j]);
    int pos = 0;
    final int maxlength = Math.max(match.length(), Math.max(bread.length(), bref.length()));
    while (pos < maxlength) {
        print_cautiously(match, pos, width);
        print_cautiously(bread, pos, width);
        print_cautiously(bref, pos, width);
        System.out.println();
        pos += width;
    }
}
Also used : Cigar(htsjdk.samtools.Cigar) ArrayList(java.util.ArrayList) GATKException(org.broadinstitute.hellbender.exceptions.GATKException) CigarElement(htsjdk.samtools.CigarElement)

Example 97 with CigarElement

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

the class HaplotypeUnitTest method testComplexDeletionAllele.

@Test
public void testComplexDeletionAllele() {
    final String bases = "ATCG" + "CCGGCCGGCC" + "ATCGATCG" + "AGGGGGA" + "AGGC";
    final ArrayList<CigarElement> h1CigarList = new ArrayList<>();
    h1CigarList.add(new CigarElement(4, CigarOperator.M));
    h1CigarList.add(new CigarElement(10, CigarOperator.I));
    h1CigarList.add(new CigarElement(8, CigarOperator.M));
    h1CigarList.add(new CigarElement(3, CigarOperator.D));
    h1CigarList.add(new CigarElement(7 + 4, CigarOperator.M));
    final Cigar h1Cigar = new Cigar(h1CigarList);
    String h1bases = "A" + "CCGGCCGGCC" + "ATCGATCG" + "AGGGGGA" + "AGGC";
    basicInsertTest("ATCG", "A", 0, h1Cigar, bases, h1bases);
    h1bases = "ATCG" + "CCGGCCGGCC" + "ATAAAG" + "AGGGGGA" + "AGGC";
    basicInsertTest("CGATC", "AAA", 6, h1Cigar, bases, h1bases);
    h1bases = "ATCG" + "CCGGCCGGCC" + "ATCGATCG" + "AGA" + "AGGC";
    basicInsertTest("GGGGG", "G", 16, h1Cigar, bases, h1bases);
}
Also used : Cigar(htsjdk.samtools.Cigar) ArrayList(java.util.ArrayList) CigarElement(htsjdk.samtools.CigarElement) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 98 with CigarElement

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

the class HaplotypeUnitTest method testComplexSNPAllele.

@Test
public void testComplexSNPAllele() {
    final String bases = "ATCG" + "CCGGCCGGCC" + "ATCGATCG" + "AGGGGGA" + "AGGC";
    final ArrayList<CigarElement> h1CigarList = new ArrayList<>();
    h1CigarList.add(new CigarElement(4, CigarOperator.M));
    h1CigarList.add(new CigarElement(10, CigarOperator.I));
    h1CigarList.add(new CigarElement(8, CigarOperator.M));
    h1CigarList.add(new CigarElement(3, CigarOperator.D));
    h1CigarList.add(new CigarElement(7 + 4, CigarOperator.M));
    final Cigar h1Cigar = new Cigar(h1CigarList);
    String h1bases = "AGCG" + "CCGGCCGGCC" + "ATCGATCG" + "AGGGGGA" + "AGGC";
    basicInsertTest("T", "G", 1, h1Cigar, bases, h1bases);
    h1bases = "ATCG" + "CCGGCCGGCC" + "ATCTATCG" + "AGGGGGA" + "AGGC";
    basicInsertTest("G", "T", 7, h1Cigar, bases, h1bases);
    h1bases = "ATCG" + "CCGGCCGGCC" + "ATCGATCG" + "AGCGGGA" + "AGGC";
    basicInsertTest("G", "C", 17, h1Cigar, bases, h1bases);
}
Also used : Cigar(htsjdk.samtools.Cigar) ArrayList(java.util.ArrayList) CigarElement(htsjdk.samtools.CigarElement) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 99 with CigarElement

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

the class HaplotypeUnitTest method testSimpleInsertionAllele.

@Test
public void testSimpleInsertionAllele() {
    final String bases = "ACTGGTCAACTGGTCAACTGGTCAACTGGTCA";
    final ArrayList<CigarElement> h1CigarList = new ArrayList<>();
    h1CigarList.add(new CigarElement(bases.length(), CigarOperator.M));
    final Cigar h1Cigar = new Cigar(h1CigarList);
    String h1bases = "AACTTCTGGTCAACTGGTCAACTGGTCAACTGGTCA";
    basicInsertTest("A", "AACTT", 0, h1Cigar, bases, h1bases);
    h1bases = "ACTGGTCAACTTACTGGTCAACTGGTCAACTGGTCA";
    basicInsertTest("A", "AACTT", 7, h1Cigar, bases, h1bases);
    h1bases = "ACTGGTCAACTGGTCAAACTTCTGGTCAACTGGTCA";
    basicInsertTest("A", "AACTT", 16, h1Cigar, bases, h1bases);
}
Also used : Cigar(htsjdk.samtools.Cigar) ArrayList(java.util.ArrayList) CigarElement(htsjdk.samtools.CigarElement) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 100 with CigarElement

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

the class PileupElement method getNextIndelCigarElement.

/**
     * Helpful function to get the immediately following cigar element, for an insertion or deletion
     *
     * if this state precedes a deletion (i.e., next position on genome) or insertion (immediately between
     * this and the next position) returns the CigarElement corresponding to this event.  Otherwise returns
     * null.
     *
     * @return a CigarElement, or null if the next alignment state isn't an insertion or deletion.
     */
private CigarElement getNextIndelCigarElement() {
    if (isBeforeDeletionStart()) {
        final CigarElement element = getNextOnGenomeCigarElement();
        Utils.validate(element != null && element.getOperator() == CigarOperator.D, () -> "Immediately before deletion but the next cigar element isn't a deletion " + element);
        return element;
    } else if (isBeforeInsertion()) {
        final CigarElement element = getBetweenNextPosition().get(0);
        Utils.validate(element.getOperator() == CigarOperator.I, () -> "Immediately before insertion but the next cigar element isn't an insertion " + element);
        return element;
    } else {
        return null;
    }
}
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