Package htsjdk.samtools

Examples of htsjdk.samtools.SAMSequenceRecord

49
50
51
52
53
54
55
56
57
58
59
    private SAMFileHeader samFileHeader;

    @BeforeClass
    public void init() {
        SAMSequenceDictionary sequenceDictionary = new SAMSequenceDictionary();
        sequenceDictionary.addSequence(new SAMSequenceRecord("1",500));
        sequenceDictionary.addSequence(new SAMSequenceRecord("2",500));

        samFileHeader = new SAMFileHeader();
        samFileHeader.setSequenceDictionary(sequenceDictionary);
    }
View Full Code Here
172
173
174
175
176
177
178
179
180
181
182
            // Write the ann and amb for this reference sequence.
            annFile = File.createTempFile("bwt",".ann");
            ambFile = File.createTempFile("bwt",".amb");

            SAMSequenceDictionary dictionary = new SAMSequenceDictionary();
            dictionary.addSequence(new SAMSequenceRecord("autogenerated",normalizedReferenceSequence.length));

            ANNWriter annWriter = new ANNWriter(annFile);
            annWriter.write(dictionary);
            annWriter.close();
View Full Code Here
102
103
104
105
106
107
108
109
110
111
112
113
114
115
        return getReferenceBases(genomeLocParser.createGenomeLoc(read));

    }

    protected byte[] getReferenceBases( GenomeLoc genomeLoc ) {
        SAMSequenceRecord sequenceInfo = reference.getSequenceDictionary().getSequence(genomeLoc.getContig());

        long start = genomeLoc.getStart();
        long stop = Math.min( genomeLoc.getStop(), sequenceInfo.getSequenceLength() );

        // Read with no aligned bases?  Return an empty array.
        if(stop - start + 1 == 0)
            return new byte[0];
View Full Code Here
47
48
49
50
51
52
53
54
55
56
57
58
59
    private static Logger logger = Logger.getLogger(SequenceDictionaryUtilsUnitTest.class);


    @DataProvider( name = "SequenceDictionaryDataProvider" )
    public Object[][] generateSequenceDictionaryTestData() {
        final SAMSequenceRecord CHRM_HG19 = new SAMSequenceRecord("chrM", 16571);
        final SAMSequenceRecord CHR_NONSTANDARD1 = new SAMSequenceRecord("NonStandard1", 8675309);
        final SAMSequenceRecord CHR_NONSTANDARD2 = new SAMSequenceRecord("NonStandard2", 8675308);

        final Class NO_COMMON_CONTIGS_EXCEPTION = UserException.IncompatibleSequenceDictionaries.class;
        final Class UNEQUAL_COMMON_CONTIGS_EXCEPTION = UserException.IncompatibleSequenceDictionaries.class;
        final Class NON_CANONICAL_HUMAN_ORDER_EXCEPTION = UserException.LexicographicallySortedSequenceDictionary.class;
        final Class OUT_OF_ORDER_EXCEPTION = UserException.IncompatibleSequenceDictionaries.class;
View Full Code Here
82
83
84
85
86
87
88
89
90
91
    @Test(dataProvider = "fastas", enabled = true && ! DEBUG)
    public void testCachingIndexedFastaReaderSequential1(File fasta, int cacheSize, int querySize) throws FileNotFoundException {
        final CachingIndexedFastaSequenceFile caching = new CachingIndexedFastaSequenceFile(fasta, getCacheSize(cacheSize), true, false);

        SAMSequenceRecord contig = caching.getSequenceDictionary().getSequence(0);
        logger.warn(String.format("Checking contig %s length %d with cache size %d and query size %d",
                contig.getSequenceName(), contig.getSequenceLength(), cacheSize, querySize));
        testSequential(caching, fasta, querySize);
    }
View Full Code Here
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
    private void testSequential(final CachingIndexedFastaSequenceFile caching, final File fasta, final int querySize) throws FileNotFoundException {
        Assert.assertTrue(caching.isPreservingCase(), "testSequential only works for case preserving CachingIndexedFastaSequenceFile readers");

        final IndexedFastaSequenceFile uncached = new IndexedFastaSequenceFile(fasta);

        SAMSequenceRecord contig = uncached.getSequenceDictionary().getSequence(0);
        for ( int i = 0; i < contig.getSequenceLength(); i += STEP_SIZE ) {
            int start = i;
            int stop = start + querySize;
            if ( stop <= contig.getSequenceLength() ) {
                ReferenceSequence cachedVal = caching.getSubsequenceAt(contig.getSequenceName(), start, stop);
                ReferenceSequence uncachedVal = uncached.getSubsequenceAt(contig.getSequenceName(), start, stop);

                Assert.assertEquals(cachedVal.getName(), uncachedVal.getName());
                Assert.assertEquals(cachedVal.getContigIndex(), uncachedVal.getContigIndex());
                Assert.assertEquals(cachedVal.getBases(), uncachedVal.getBases());
            }
        }

        // asserts for efficiency.  We are going to make contig.length / STEP_SIZE queries
        // at each of range: start -> start + querySize against a cache with size of X.
        // we expect to hit the cache each time range falls within X.  We expect a hit
        // on the cache if range is within X.  Which should happen at least (X - query_size * 2) / STEP_SIZE
        // times.
        final int minExpectedHits = (int)Math.floor((Math.min(caching.getCacheSize(), contig.getSequenceLength()) - querySize * 2.0) / STEP_SIZE);
        caching.printEfficiency(Priority.WARN);
        Assert.assertTrue(caching.getCacheHits() >= minExpectedHits, "Expected at least " + minExpectedHits + " cache hits but only got " + caching.getCacheHits());

    }
View Full Code Here
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
    @Test(dataProvider = "fastas", enabled = true && ! DEBUG)
    public void testCachingIndexedFastaReaderTwoStage(File fasta, int cacheSize, int querySize) throws FileNotFoundException {
        final IndexedFastaSequenceFile uncached = new IndexedFastaSequenceFile(fasta);
        final CachingIndexedFastaSequenceFile caching = new CachingIndexedFastaSequenceFile(fasta, getCacheSize(cacheSize), true, false);

        SAMSequenceRecord contig = uncached.getSequenceDictionary().getSequence(0);

        int middleStart = (contig.getSequenceLength() - querySize) / 2;
        int middleStop = middleStart + querySize;

        logger.warn(String.format("Checking contig %s length %d with cache size %d and query size %d with intermediate query",
                contig.getSequenceName(), contig.getSequenceLength(), cacheSize, querySize));

        for ( int i = 0; i < contig.getSequenceLength(); i += 10 ) {
            int start = i;
            int stop = start + querySize;
            if ( stop <= contig.getSequenceLength() ) {
                ReferenceSequence grabMiddle = caching.getSubsequenceAt(contig.getSequenceName(), middleStart, middleStop);
                ReferenceSequence cachedVal = caching.getSubsequenceAt(contig.getSequenceName(), start, stop);
                ReferenceSequence uncachedVal = uncached.getSubsequenceAt(contig.getSequenceName(), start, stop);

                Assert.assertEquals(cachedVal.getName(), uncachedVal.getName());
                Assert.assertEquals(cachedVal.getContigIndex(), uncachedVal.getContigIndex());
                Assert.assertEquals(cachedVal.getBases(), uncachedVal.getBases());
            }
View Full Code Here
119
120
121
122
123
124
125
126
127
128
129
130
    @Test
    public void testContigHasColon() {
        SAMFileHeader header = new SAMFileHeader();
        header.setSortOrder(htsjdk.samtools.SAMFileHeader.SortOrder.coordinate);
        SAMSequenceDictionary dict = new SAMSequenceDictionary();
        SAMSequenceRecord rec = new SAMSequenceRecord("c:h:r1", 10);
        rec.setSequenceLength(10);
        dict.addSequence(rec);
        header.setSequenceDictionary(dict);

        final GenomeLocParser myGenomeLocParser = new GenomeLocParser(header.getSequenceDictionary());
        GenomeLoc loc = myGenomeLocParser.parseGenomeLoc("c:h:r1:4-5");
View Full Code Here
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
        GenomeLoc cur = context.getLocation();

        if ( verbose && showSkipped ) {
            for(long i = context.getSkippedBases(); i >= 0; i--) {
                SAMSequenceDictionary dictionary = getToolkit().getReferenceDataSource().getReference().getSequenceDictionary();
                SAMSequenceRecord contig = dictionary.getSequence(cur.getContig());
                if(cur.getStop() < contig.getSequenceLength())
                    cur = getToolkit().getGenomeLocParser().incPos(cur,1);
                else
                    cur = getToolkit().getGenomeLocParser().createGenomeLoc(dictionary.getSequence(contig.getSequenceIndex()+1).getSequenceName(),1,1);
                out.printf("%s: skipped%n", cur);

            }
        }
View Full Code Here
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
    @BeforeClass
    public void init() {
        header = ArtificialSAMUtils.createArtificialSamHeader(3, 1, ArtificialSAMUtils.DEFAULT_READ_LENGTH * 2);
        GenomeLocParser genomeLocParser = new GenomeLocParser(header.getSequenceDictionary());
        SAMSequenceRecord record;

        record = header.getSequence(0);
        firstContig = genomeLocParser.createGenomeLoc(record.getSequenceName(), 1, record.getSequenceLength());
        record = header.getSequence(1);
        secondContig = genomeLocParser.createGenomeLoc(record.getSequenceName(), 1, record.getSequenceLength());

        readMapped = createMappedRead("mapped", 1);

        readUnmappedFlag = createMappedRead("unmappedFlagged", 2);
        readUnmappedFlag.setReadUnmappedFlag(true);
View Full Code Here
TOP

Related Classes of htsjdk.samtools.SAMSequenceRecord

  • org.broad.igv.plugin.mongovariant.VariantReviewSource
  • org.broadinstitute.gatk.engine.alignment.bwa.BWTFiles
  • org.broadinstitute.gatk.engine.datasources.providers.ReferenceView
  • org.broadinstitute.gatk.engine.datasources.providers.ReferenceViewTemplate
  • org.broadinstitute.gatk.engine.datasources.reads.BAMScheduler
  • org.broadinstitute.gatk.engine.datasources.reads.IntervalOverlapFilteringIteratorUnitTest
  • org.broadinstitute.gatk.engine.iterators.VerifyingSamIteratorUnitTest
  • org.broadinstitute.gatk.engine.refdata.tracks.IndexDictionaryUtils
  • org.broadinstitute.gatk.tools.walkers.qc.CountRODs
  • org.broadinstitute.gatk.utils.fasta.CachingIndexedFastaSequenceFile
    • 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.