Package htsjdk.samtools

Examples of htsjdk.samtools.Cigar

207
208
209
210
211
212
213
214
215
216
217
218
    }


    private List<Haplotype> assemble(final ReadThreadingAssembler assembler, final byte[] refBases, final GenomeLoc loc, final List<GATKSAMRecord> reads) {
        final Haplotype refHaplotype = new Haplotype(refBases, true);
        final Cigar c = new Cigar();
        c.add(new CigarElement(refHaplotype.getBases().length, CigarOperator.M));
        refHaplotype.setCigar(c);

        final ActiveRegion activeRegion = new ActiveRegion(loc, null, true, genomeLocParser, 0);
        activeRegion.addAll(reads);
//        logger.warn("Assembling " + activeRegion + " with " + engine);
View Full Code Here
146
147
148
149
150
151
152
153
154
155
156
157
158
159
        final String notPaddedsRef = SW_PAD + new String(notPaddedRef) + SW_PAD;
        final String notpaddedsHap = SW_PAD + new String(notPaddedHap) + SW_PAD;
        final SmithWaterman paddedAlignment = new SWPairwiseAlignment( paddedsRef.getBytes(), paddedsHap.getBytes(), CigarUtils.NEW_SW_PARAMETERS );
        final SmithWaterman notPaddedAlignment = new SWPairwiseAlignment( notPaddedsRef.getBytes(), notpaddedsHap.getBytes(), CigarUtils.NEW_SW_PARAMETERS );
        //Now verify that the two sequences have the same alignment and not match positions.
        Cigar rawPadded = paddedAlignment.getCigar();
        Cigar notPadded= notPaddedAlignment.getCigar();
        List<CigarElement> paddedC=rawPadded.getCigarElements();
        List<CigarElement> notPaddedC=notPadded.getCigarElements();
        Assert.assertEquals(paddedC.size(),notPaddedC.size());
        for(int i=0;i<notPaddedC.size();i++)
        {
            CigarElement pc=paddedC.get(i);
            CigarElement npc=notPaddedC.get(i);
View Full Code Here
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
        for ( final VariantContext vc : stateForTesting )
            addVC(vc);
    }

    protected void processCigarForInitialEvents() {
        final Cigar cigar = haplotype.getCigar();
        final byte[] alignment = haplotype.getBases();

        int refPos = haplotype.getAlignmentStartHapwrtRef();
        if( refPos < 0 ) {
            return;
        } // Protection against SW failures

        final List<VariantContext> proposedEvents = new ArrayList<>();

        int alignmentPos = 0;

        for( int cigarIndex = 0; cigarIndex < cigar.numCigarElements(); cigarIndex++ ) {
            final CigarElement ce = cigar.getCigarElement(cigarIndex);
            final int elementLength = ce.getLength();
            switch( ce.getOperator() ) {
                case I:
                {
                    if( refPos > 0 ) { // protect against trying to create insertions/deletions at the beginning of a contig
                        final List<Allele> insertionAlleles = new ArrayList<Allele>();
                        final int insertionStart = refLoc.getStart() + refPos - 1;
                        final byte refByte = ref[refPos-1];
                        if( BaseUtils.isRegularBase(refByte) ) {
                            insertionAlleles.add( Allele.create(refByte, true) );
                        }
                        if( cigarIndex == 0 || cigarIndex == cigar.getCigarElements().size() - 1 ) {
                            // if the insertion isn't completely resolved in the haplotype, skip it
                            // note this used to emit SYMBOLIC_UNASSEMBLED_EVENT_ALLELE but that seems dangerous
                        } else {
                            byte[] insertionBases = new byte[]{};
                            insertionBases = ArrayUtils.add(insertionBases, ref[refPos - 1]); // add the padding base
View Full Code Here
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
        // compute the smith-waterman alignment of read -> haplotype
        final SWPairwiseAlignment swPairwiseAlignment = new SWPairwiseAlignment(haplotype.getBases(), originalRead.getReadBases(), CigarUtils.NEW_SW_PARAMETERS);
        if ( swPairwiseAlignment.getAlignmentStart2wrt1() == -1 )
            // sw can fail (reasons not clear) so if it happens just don't realign the read
            return originalRead;
        final Cigar swCigar = consolidateCigar(swPairwiseAlignment.getCigar());

        // since we're modifying the read we need to clone it
        final GATKSAMRecord read = (GATKSAMRecord)originalRead.clone();

        // only informative reads are given the haplotype tag to enhance visualization
        if ( isInformative )
            read.setAttribute(HAPLOTYPE_TAG, haplotype.hashCode());

        // compute here the read starts w.r.t. the reference from the SW result and the hap -> ref cigar
        final Cigar extendedHaplotypeCigar = haplotype.getConsolidatedPaddedCigar(1000);
        final int readStartOnHaplotype = calcFirstBaseMatchingReferenceInCigar(extendedHaplotypeCigar, swPairwiseAlignment.getAlignmentStart2wrt1());
        final int readStartOnReference = referenceStart + haplotype.getAlignmentStartHapwrtRef() + readStartOnHaplotype;
        read.setAlignmentStart(readStartOnReference);
        read.resetSoftStartAndEnd();

        // compute the read -> ref alignment by mapping read -> hap -> ref from the
        // SW of read -> hap mapped through the given by hap -> ref
        final Cigar haplotypeToRef = trimCigarByBases(extendedHaplotypeCigar, swPairwiseAlignment.getAlignmentStart2wrt1(), extendedHaplotypeCigar.getReadLength() - 1);
        final Cigar readToRefCigarRaw = applyCigarToCigar(swCigar, haplotypeToRef);
        final Cigar readToRefCigarClean = cleanUpCigar(readToRefCigarRaw);
        final Cigar readToRefCigar = leftAlignIndel(readToRefCigarClean, haplotype.getBases(),
                originalRead.getReadBases(), swPairwiseAlignment.getAlignmentStart2wrt1(), 0, true);

        read.setCigar(readToRefCigar);

        if ( readToRefCigar.getReadLength() != read.getReadLength() )
            throw new IllegalStateException("Cigar " + readToRefCigar + " with read length " + readToRefCigar.getReadLength()
                    + " != read length " + read.getReadLength() + " for read " + read.format() + "\nhapToRef " + haplotypeToRef + " length " + haplotypeToRef.getReadLength() + "/" + haplotypeToRef.getReferenceLength()
                    + "\nreadToHap " + swCigar + " length " + swCigar.getReadLength() + "/" + swCigar.getReferenceLength());

        return read;
    }
View Full Code Here
285
286
287
288
289
290
291
292
293
294
295
296
297
        MismatchCount mc = new MismatchCount();

        int readIdx = 0;
        final int endOnRead = startOnRead + nReadBases - 1; // index of the last base on read we want to count (note we are including soft-clipped bases with this math)
        final byte[] readSeq = r.getReadBases();
        final Cigar c = r.getCigar();
        final byte[] readQuals = r.getBaseQualities();
        for (final CigarElement ce : c.getCigarElements()) {

            if (readIdx > endOnRead)
                break;

            final int elementLength = ce.getLength();
View Full Code Here
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
     * @return number of continuous alignment blocks (i.e. 'M' elements of the cigar; all indel and clipping elements are ignored).
     */
    @Ensures("result >= 0")
    public static int getNumAlignmentBlocks(final SAMRecord r) {
        if ( r == null ) throw new IllegalArgumentException("read cannot be null");
        final Cigar cigar = r.getCigar();
        if (cigar == null) return 0;

        int n = 0;
        for (final CigarElement e : cigar.getCigarElements()) {
            if (ALIGNED_TO_GENOME_OPERATORS.contains(e.getOperator()))
                n++;
        }

        return n;
View Full Code Here
375
376
377
378
379
380
381
382
383
384
385
386
387
388
     * @param r a non-null GATKSAMRecord
     * @return the number of bases aligned to the genome in R, including soft clipped bases
     */
    public static int getNumAlignedBasesCountingSoftClips(final GATKSAMRecord r) {
        int n = 0;
        final Cigar cigar = r.getCigar();
        if (cigar == null) return 0;

        for (final CigarElement e : cigar.getCigarElements())
            if (ALIGNED_TO_GENOME_PLUS_SOFTCLIPS.contains(e.getOperator()))
                n += e.getLength();

        return n;
    }
View Full Code Here
398
399
400
401
402
403
404
405
406
407
408
409
410
411
    @Ensures("result >= 0")
    public static int getNumHardClippedBases(final SAMRecord r) {
        if ( r == null ) throw new IllegalArgumentException("Read cannot be null");

        int n = 0;
        final Cigar cigar = r.getCigar();
        if (cigar == null) return 0;

        for (final CigarElement e : cigar.getCigarElements())
            if (e.getOperator() == CigarOperator.H)
                n += e.getLength();

        return n;
    }
View Full Code Here
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
        // fast check to determine if there's anything worth doing before we create new Cigar and actually do some work
        if ( ! needsConsolidation(c) )
            return c;

        final Cigar returnCigar = new Cigar();
        int sumLength = 0;
        CigarElement lastElement = null;

        for( final CigarElement cur : c.getCigarElements() ) {
            if ( cur.getLength() == 0 )
                continue; // don't add elements of 0 length

            if ( lastElement != null && lastElement.getOperator() != cur.getOperator() ) {
                returnCigar.add(new CigarElement(sumLength, lastElement.getOperator()));
                sumLength = 0;
            }

            sumLength += cur.getLength();
            lastElement = cur;
        }

        if ( sumLength > 0 ) {
            returnCigar.add(new CigarElement(sumLength, lastElement.getOperator()));
        }

        return returnCigar;
    }
View Full Code Here
119
120
121
122
123
124
125
126
127
128
129
        if ( getCigar() == null ) throw new IllegalArgumentException("Cannot trim haplotype without a cigar " + this);

        final int newStart = loc.getStart() - this.genomeLocation.getStart();
        final int newStop = newStart + loc.size() - 1;
        final byte[] newBases = AlignmentUtils.getBasesCoveringRefInterval(newStart, newStop, getBases(), 0, getCigar());
        final Cigar newCigar = AlignmentUtils.trimCigarByReference(getCigar(), newStart, newStop);

        if ( newBases == null || AlignmentUtils.startsOrEndsWithInsertionOrDeletion(newCigar) )
            // we cannot meaningfully chop down the haplotype, so return null
            return null;
View Full Code Here
TOP

Related Classes of htsjdk.samtools.Cigar

  • bcbio.gatk.tools.walkers.LeftAlignVariants
  • bcbio.variation.annotate.ReadMeanLen
  • bcbio.variation.annotate.ReadMeanPos
  • bcbio.variation.annotate.ReadPosEndDist
  • org.broad.igv.ga4gh.Ga4ghAlignment
  • org.broad.igv.sam.reader.SamUtils
  • org.broadinstitute.gatk.engine.alignment.bwa.java.AlignmentMatchSequence
  • org.broadinstitute.gatk.engine.datasources.reads.TheoreticalMinimaBenchmark
  • org.broadinstitute.gatk.engine.filters.BadCigarFilter
  • org.broadinstitute.gatk.engine.filters.MalformedReadFilterUnitTest
    • Copyright © 2018 www.massapicom. All rights reserved.
    All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.