Examples of org.broadinstitute.gatk.utils.sam.GATKSAMRecord

Package org.broadinstitute.gatk.utils.sam

Examples of org.broadinstitute.gatk.utils.sam.GATKSAMRecord

org.broadinstitute.gatk.utils.sam.GATKSAMRecord
@author ebanks, depristoGATKSAMRecord this class extends the samtools BAMRecord class (and SAMRecord) and caches important (and oft-accessed) data that's not already cached by the SAMRecord class IMPORTANT NOTE: Because ReadGroups are not set through the SAMRecord, if they are ever modified externally then one must also invoke the setReadGroup() method here to ensure that the cache is kept up-to-date. WARNING -- GATKSAMRecords cache several values (that are expensive to compute) that depending on the inferred insert size and alignment starts and stops of this read and its mate. Changing these values in any way will invalidate the cached value. However, we do not monitor those setter functions, so modifying a GATKSAMRecord in any way may result in stale cached values.


        return tests.toArray(new Object[][]{});
    }


    private GATKSAMRecord modifyBaseQualities(final GATKSAMRecord read, final int startOffset, final int length) throws Exception {
        final GATKSAMRecord readWithLowQuals = (GATKSAMRecord)read.clone();
        final byte[] withLowQuals = Arrays.copyOf(read.getBaseQualities(), read.getBaseQualities().length);
        for ( int i = startOffset; i < startOffset + length; i++ )
            withLowQuals[i] = (byte)(read.getBaseQualities()[i] + (i % 2 == 0 ? -1 : 0));
        readWithLowQuals.setBaseQualities(withLowQuals);
        return readWithLowQuals;
    }

View Full Code Here

    private GATKSAMRecord makeOverlappingRead(final String leftFlank, final int leftQual, final String overlapBases,
                                              final byte[] overlapQuals, final String rightFlank, final int rightQual,
                                              final int alignmentStart) {
        final String bases = leftFlank + overlapBases + rightFlank;
        final int readLength = bases.length();
        final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "myRead", 0, alignmentStart, readLength);
        final byte[] leftQuals = Utils.dupBytes((byte) leftQual, leftFlank.length());
        final byte[] rightQuals = Utils.dupBytes((byte) rightQual, rightFlank.length());
        final byte[] quals = Utils.concat(leftQuals, overlapQuals, rightQuals);
        read.setCigarString(readLength + "M");
        read.setReadBases(bases.getBytes());
        for ( final EventType type : EventType.values() )
            read.setBaseQualities(quals, type);
        read.setReadGroup(rgForMerged);
        read.setMappingQuality(60);
        return read;
    }

View Full Code Here

        return read;
    }


    @Test(enabled = !DEBUG, dataProvider = "MergeFragmentsTest")
    public void testMergingTwoReads(final String name, final GATKSAMRecord read1, final GATKSAMRecord read2, final GATKSAMRecord expectedMerged) {
        final GATKSAMRecord actual = FragmentUtils.mergeOverlappingPairedFragments(read1, read2);


        if ( expectedMerged == null ) {
            Assert.assertNull(actual, "Expected reads not to merge, but got non-null result from merging");
        } else {
            Assert.assertTrue(actual.isStrandless(), "Merged reads should be strandless");
            Assert.assertNotNull(actual, "Expected reads to merge, but got null result from merging");
            // I really care about the bases, the quals, the CIGAR, and the read group tag
            Assert.assertEquals(actual.getCigarString(), expectedMerged.getCigarString());
            Assert.assertEquals(actual.getReadBases(), expectedMerged.getReadBases());
            Assert.assertEquals(actual.getReadGroup(), expectedMerged.getReadGroup());
            Assert.assertEquals(actual.getMappingQuality(), expectedMerged.getMappingQuality());
            for ( final EventType type : EventType.values() )
                Assert.assertEquals(actual.getBaseQualities(type), expectedMerged.getBaseQualities(type), "Failed base qualities for event type " + type);
        }
    }

View Full Code Here

    public void testHardClippingBeforeMerge() {
        final String common = Utils.dupString("A", 10);
        final byte[] commonQuals = Utils.dupBytes((byte)30, common.length());
        final String adapter    = "NNNN";


        final GATKSAMRecord read1 = makeOverlappingRead(adapter, 30, common, commonQuals, "", 30, 10);
        final GATKSAMRecord read2 = makeOverlappingRead("", 30, common, commonQuals, adapter, 30, 10);
        final GATKSAMRecord expectedMerged = makeOverlappingRead("", 30, common, commonQuals, "", 30, 10);
        read1.setCigarString("4S" + common.length() + "M");
        read1.setProperPairFlag(true);
        read1.setReadPairedFlag(true);
        read1.setFirstOfPairFlag(true);
        read1.setReadNegativeStrandFlag(true);
        read1.setMateNegativeStrandFlag(false);
        read1.setMateAlignmentStart(read2.getAlignmentStart());
        read2.setCigarString(common.length() + "M4S");
        read2.setProperPairFlag(true);
        read2.setReadPairedFlag(true);
        read2.setFirstOfPairFlag(false);
        read2.setReadNegativeStrandFlag(false);
        read2.setMateNegativeStrandFlag(true);
        read2.setMateAlignmentStart(read1.getAlignmentStart());


        final int insertSize = common.length() - 1;
        read1.setInferredInsertSize(-insertSize);
        read2.setInferredInsertSize(insertSize);


        final GATKSAMRecord actual = FragmentUtils.mergeOverlappingPairedFragments(read1, read2);
        Assert.assertEquals(actual.getCigarString(), expectedMerged.getCigarString());
        Assert.assertEquals(actual.getReadBases(), expectedMerged.getReadBases());
        Assert.assertEquals(actual.getReadGroup(), expectedMerged.getReadGroup());
        Assert.assertEquals(actual.getMappingQuality(), expectedMerged.getMappingQuality());
        for ( final EventType type : EventType.values() )
            Assert.assertEquals(actual.getBaseQualities(type), expectedMerged.getBaseQualities(type), "Failed base qualities for event type " + type);
    }

View Full Code Here


    @Test(enabled = true)
    public void testHardClippingBeforeMergeResultingInCompletelyContainedSecondRead() {
        final String adapter    = "NNNN";


        final GATKSAMRecord read1 = makeOverlappingRead(adapter, 30, Utils.dupString("A", 10), Utils.dupBytes((byte)30, 10), "", 30, 10);
        final GATKSAMRecord read2 = makeOverlappingRead("", 30, Utils.dupString("A", 7), Utils.dupBytes((byte)30, 7), adapter, 30, 10);
        read1.setCigarString("4S10M");
        read1.setProperPairFlag(true);
        read1.setFirstOfPairFlag(true);
        read1.setReadNegativeStrandFlag(true);
        read1.setMateAlignmentStart(10);
        read2.setCigarString("7M4S");
        read2.setProperPairFlag(true);
        read2.setFirstOfPairFlag(false);
        read2.setReadNegativeStrandFlag(false);


        final int insertSize = 7 - 1;
        read1.setInferredInsertSize(insertSize);
        read2.setInferredInsertSize(-insertSize);


        final GATKSAMRecord actual = FragmentUtils.mergeOverlappingPairedFragments(read1, read2);
        Assert.assertNull(actual);
    }

View Full Code Here

    @Test(dataProvider = "MergeFragmentsOffContig")
    public void testMergeFragmentsOffContig(final int pre1, final int post1, final int pre2, final int post2) {
        final int contigSize = 10;
        final SAMFileHeader header = ArtificialSAMUtils.createArtificialSamHeader(1, 0, contigSize);


        final GATKSAMRecord read1 = createReadOffContig(header, false, pre1, post1);
        final GATKSAMRecord read2 = createReadOffContig(header, true, pre2, post2);


        final GATKSAMRecord merged = FragmentUtils.mergeOverlappingPairedFragments(read1, read2);
    }

View Full Code Here

    }


    private GATKSAMRecord createReadOffContig(final SAMFileHeader header, final boolean negStrand, final int pre, final int post) {
        final int contigLen = header.getSequence(0).getSequenceLength();
        final int readLen = pre + contigLen + post;
        final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "read1", 0, 1, readLen);
        read.setAlignmentStart(1);
        read.setCigar(TextCigarCodec.getSingleton().decode(pre + "S" + contigLen + "M" + post + "S"));
        read.setBaseQualities(Utils.dupBytes((byte) 30, readLen));
        read.setReadBases(Utils.dupBytes((byte)'A', readLen));
        read.setMappingQuality(60);
        read.setMateAlignmentStart(1);
        read.setProperPairFlag(true);
        read.setReadPairedFlag(true);
        read.setInferredInsertSize(30);
        read.setReadNegativeStrandFlag(negStrand);
        read.setMateNegativeStrandFlag(! negStrand);
        read.setReadGroup(new GATKSAMReadGroupRecord("foo"));
        return read;
    }

View Full Code Here

        final String allOverlappingBases = "ACGTACGTGGAACCTTAG";
        for ( int overlapSize = 1; overlapSize < allOverlappingBases.length(); overlapSize++ ) {
            final String overlappingBases = allOverlappingBases.substring(0, overlapSize);
            final byte[] overlappingBaseQuals = new byte[overlapSize];
            for ( int i = 0; i < overlapSize; i++ ) overlappingBaseQuals[i] = highQuality;
            final GATKSAMRecord read1  = makeOverlappingRead(leftFlank, highQuality, overlappingBases, overlappingBaseQuals, "", highQuality, 1);
            final GATKSAMRecord read2  = makeOverlappingRead("", highQuality, overlappingBases, overlappingBaseQuals, rightFlank, highQuality, leftFlank.length() + 1);
            tests.add(new Object[]{read1, read2, overlapSize});
        }


        return tests.toArray(new Object[][]{});
    }

View Full Code Here

        private FragmentUtilsTest(String name, int readLen, int leftStart, int rightStart,
                                  boolean leftIsFirst, boolean leftIsNegative) {
            super(FragmentUtilsTest.class, String.format("%s-leftIsFirst:%b-leftIsNegative:%b", name, leftIsFirst, leftIsNegative));


            List<GATKSAMRecord> pair = ArtificialSAMUtils.createPair(header, "readpair", readLen, leftStart, rightStart, leftIsFirst, leftIsNegative);
            GATKSAMRecord left = pair.get(0);
            GATKSAMRecord right = pair.get(1);


            for ( int pos = leftStart; pos < rightStart + readLen; pos++) {
                boolean posCoveredByLeft = pos >= left.getAlignmentStart() && pos <= left.getAlignmentEnd();
                boolean posCoveredByRight = pos >= right.getAlignmentStart() && pos <= right.getAlignmentEnd();


                if ( posCoveredByLeft || posCoveredByRight ) {
                    List<GATKSAMRecord> reads = new ArrayList<GATKSAMRecord>();
                    List<Integer> offsets = new ArrayList<Integer>();


                    if ( posCoveredByLeft ) {
                        reads.add(left);
                        offsets.add(pos - left.getAlignmentStart());
                    }


                    if ( posCoveredByRight ) {
                        reads.add(right);
                        offsets.add(pos - right.getAlignmentStart());
                    }


                    boolean shouldBeFragment = posCoveredByLeft && posCoveredByRight;
                    ReadBackedPileup pileup = new ReadBackedPileupImpl(null, reads, offsets);
                    TestState testState = new TestState(shouldBeFragment ? 0 : 1, shouldBeFragment ? 1 : 0, pileup, null);
                    statesForPileup.add(testState);
                }


                TestState testState = left.getAlignmentEnd() >= right.getAlignmentStart() ? new TestState(0, 1, null, pair) : new TestState(2, 0, null, pair);
                statesForReads.add(testState);
            }
        }

View Full Code Here

        if ( read.getReadPairedFlag() ) {
            // paired
            final GenomeLoc readMateLoc = engine.getGenomeLocParser().createGenomeLoc(read.getMateReferenceName(), read.getMateAlignmentStart(), read.getMateAlignmentStart());


            for (List<GATKSAMRecord> reads : readSets) {
                GATKSAMRecord key = reads.get(0);


                // read and key start at the same place, and either the this read and the key
                // share a mate location or the read is flagged as a duplicate
                if ( read.getAlignmentStart() == key.getAlignmentStart() && key.getReadPairedFlag() && ( key.getDuplicateReadFlag() || read.getDuplicateReadFlag() ) ) {
                    // at least one has to be marked as a duplicate
                    final GenomeLoc keyMateLoc = engine.getGenomeLocParser().createGenomeLoc(key.getMateReferenceName(), key.getMateAlignmentStart(), key.getMateAlignmentStart());
                    if ( readMateLoc.compareTo(keyMateLoc) == 0 ) {
                        // we are at the same position as the dup and have the same mat pos, it's a dup
                        if (DEBUG) logger.debug(String.format("  => Adding read to dups list: %s %d %s vs. %s", read, reads.size(), readMateLoc, keyMateLoc));
                        return reads;
                    }
                }
            }
        } else {
            for (List<GATKSAMRecord> reads : readSets) {
                GATKSAMRecord key = reads.get(0);
                boolean v = (! key.getReadPairedFlag()) && read.getAlignmentStart() == key.getAlignmentStart() && ( key.getDuplicateReadFlag() || read.getDuplicateReadFlag() ) && read.getReadLength() == key.getReadLength();
                //System.out.printf("%s %s %b %b %d %d %d %d => %b%n",
                //        read.getReadPairedFlag(), key.getReadPairedFlag(), read.getDuplicateReadFlag(), key.getDuplicateReadFlag(),
                //        read.getAlignmentStart(), key.getAlignmentStart(), read.getReadLength(), key.getReadLength(), v);
                if ( v ) {
                    //System.out.printf("Returning reads...%n");

View Full Code Here

0 1 2 3 4 5 6 7 8 9

TOP

Related Classes of org.broadinstitute.gatk.utils.sam.GATKSAMRecord

bcbio.gatk.tools.walkers.annotator.DepthOfCoverage

bcbio.gatk.tools.walkers.annotator.FisherStrand

org.broadinstitute.gatk.engine.datasources.providers.LocusViewTemplate

org.broadinstitute.gatk.engine.downsampling.AlleleBiasedDownsamplingUtilsUnitTest

org.broadinstitute.gatk.engine.filters.BadCigarFilterUnitTest

org.broadinstitute.gatk.engine.filters.MalformedReadFilterUnitTest

org.broadinstitute.gatk.engine.io.stubs.SAMFileWriterStub

org.broadinstitute.gatk.engine.traversals.TraverseDuplicates

org.broadinstitute.gatk.tools.walkers.annotator.DepthPerAlleleBySample

org.broadinstitute.gatk.tools.walkers.annotator.HaplotypeScore

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.