Search in sources :

Example 11 with CachingIndexedFastaSequenceFile

use of org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile in project gatk-protected by broadinstitute.

the class HaplotypeCallerEngineUnitTest method testIsActive.

@Test
public void testIsActive() throws IOException {
    final File testBam = new File(NA12878_20_21_WGS_bam);
    final File reference = new File(b37_reference_20_21);
    final SimpleInterval shardInterval = new SimpleInterval("20", 10000000, 10001000);
    final SimpleInterval paddedShardInterval = new SimpleInterval(shardInterval.getContig(), shardInterval.getStart() - 100, shardInterval.getEnd() + 100);
    final HaplotypeCallerArgumentCollection hcArgs = new HaplotypeCallerArgumentCollection();
    // We expect isActive() to return 1.0 for the sites below, and 0.0 for all other sites
    final List<SimpleInterval> expectedActiveSites = Arrays.asList(new SimpleInterval("20", 9999996, 9999996), new SimpleInterval("20", 9999997, 9999997), new SimpleInterval("20", 10000117, 10000117), new SimpleInterval("20", 10000211, 10000211), new SimpleInterval("20", 10000439, 10000439), new SimpleInterval("20", 10000598, 10000598), new SimpleInterval("20", 10000694, 10000694), new SimpleInterval("20", 10000758, 10000758), new SimpleInterval("20", 10001019, 10001019));
    try (final ReadsDataSource reads = new ReadsDataSource(testBam.toPath());
        final ReferenceDataSource ref = new ReferenceFileSource(reference);
        final CachingIndexedFastaSequenceFile referenceReader = new CachingIndexedFastaSequenceFile(reference)) {
        final HaplotypeCallerEngine hcEngine = new HaplotypeCallerEngine(hcArgs, reads.getHeader(), referenceReader);
        List<ReadFilter> hcFilters = HaplotypeCallerEngine.makeStandardHCReadFilters();
        hcFilters.forEach(filter -> filter.setHeader(reads.getHeader()));
        ReadFilter hcCombinedFilter = hcFilters.get(0);
        for (int i = 1; i < hcFilters.size(); ++i) {
            hcCombinedFilter = hcCombinedFilter.and(hcFilters.get(i));
        }
        final Iterator<GATKRead> readIter = new ReadFilteringIterator(reads.query(paddedShardInterval), hcCombinedFilter);
        final LocusIteratorByState libs = new LocusIteratorByState(readIter, DownsamplingMethod.NONE, false, ReadUtils.getSamplesFromHeader(reads.getHeader()), reads.getHeader(), false);
        libs.forEachRemaining(pileup -> {
            final SimpleInterval pileupInterval = new SimpleInterval(pileup.getLocation());
            final ReferenceContext pileupRefContext = new ReferenceContext(ref, pileupInterval);
            final ActivityProfileState isActiveResult = hcEngine.isActive(pileup, pileupRefContext, new FeatureContext(null, pileupInterval));
            final double expectedIsActiveValue = expectedActiveSites.contains(pileupInterval) ? 1.0 : 0.0;
            Assert.assertEquals(isActiveResult.isActiveProb(), expectedIsActiveValue, "Wrong isActive probability for site " + pileupInterval);
        });
    }
}
Also used : ReadFilteringIterator(org.broadinstitute.hellbender.utils.iterators.ReadFilteringIterator) GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) ActivityProfileState(org.broadinstitute.hellbender.utils.activityprofile.ActivityProfileState) CachingIndexedFastaSequenceFile(org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile) LocusIteratorByState(org.broadinstitute.hellbender.utils.locusiterator.LocusIteratorByState) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) ReadFilter(org.broadinstitute.hellbender.engine.filters.ReadFilter) CachingIndexedFastaSequenceFile(org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile) File(java.io.File) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 12 with CachingIndexedFastaSequenceFile

use of org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile in project gatk by broadinstitute.

the class HaplotypeCaller method getReferenceReader.

private static CachingIndexedFastaSequenceFile getReferenceReader(ReferenceInputArgumentCollection referenceArguments) {
    final CachingIndexedFastaSequenceFile referenceReader;
    final File reference = new File(referenceArguments.getReferenceFileName());
    try {
        referenceReader = new CachingIndexedFastaSequenceFile(reference);
    } catch (FileNotFoundException e) {
        throw new UserException.CouldNotReadInputFile(reference, e);
    }
    return referenceReader;
}
Also used : FileNotFoundException(java.io.FileNotFoundException) CachingIndexedFastaSequenceFile(org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile) UserException(org.broadinstitute.hellbender.exceptions.UserException) CachingIndexedFastaSequenceFile(org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile) File(java.io.File) ReferenceSequenceFile(htsjdk.samtools.reference.ReferenceSequenceFile)

Example 13 with CachingIndexedFastaSequenceFile

use of org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile in project gatk by broadinstitute.

the class HaplotypeCallerEngineUnitTest method testIsActive.

@Test
public void testIsActive() throws IOException {
    final File testBam = new File(NA12878_20_21_WGS_bam);
    final File reference = new File(b37_reference_20_21);
    final SimpleInterval shardInterval = new SimpleInterval("20", 10000000, 10001000);
    final SimpleInterval paddedShardInterval = new SimpleInterval(shardInterval.getContig(), shardInterval.getStart() - 100, shardInterval.getEnd() + 100);
    final HaplotypeCallerArgumentCollection hcArgs = new HaplotypeCallerArgumentCollection();
    // We expect isActive() to return 1.0 for the sites below, and 0.0 for all other sites
    final List<SimpleInterval> expectedActiveSites = Arrays.asList(new SimpleInterval("20", 9999996, 9999996), new SimpleInterval("20", 9999997, 9999997), new SimpleInterval("20", 10000117, 10000117), new SimpleInterval("20", 10000211, 10000211), new SimpleInterval("20", 10000439, 10000439), new SimpleInterval("20", 10000598, 10000598), new SimpleInterval("20", 10000694, 10000694), new SimpleInterval("20", 10000758, 10000758), new SimpleInterval("20", 10001019, 10001019));
    try (final ReadsDataSource reads = new ReadsDataSource(testBam.toPath());
        final ReferenceDataSource ref = new ReferenceFileSource(reference);
        final CachingIndexedFastaSequenceFile referenceReader = new CachingIndexedFastaSequenceFile(reference)) {
        final HaplotypeCallerEngine hcEngine = new HaplotypeCallerEngine(hcArgs, reads.getHeader(), referenceReader);
        List<ReadFilter> hcFilters = HaplotypeCallerEngine.makeStandardHCReadFilters();
        hcFilters.forEach(filter -> filter.setHeader(reads.getHeader()));
        ReadFilter hcCombinedFilter = hcFilters.get(0);
        for (int i = 1; i < hcFilters.size(); ++i) {
            hcCombinedFilter = hcCombinedFilter.and(hcFilters.get(i));
        }
        final Iterator<GATKRead> readIter = new ReadFilteringIterator(reads.query(paddedShardInterval), hcCombinedFilter);
        final LocusIteratorByState libs = new LocusIteratorByState(readIter, DownsamplingMethod.NONE, false, ReadUtils.getSamplesFromHeader(reads.getHeader()), reads.getHeader(), false);
        libs.forEachRemaining(pileup -> {
            final SimpleInterval pileupInterval = new SimpleInterval(pileup.getLocation());
            final ReferenceContext pileupRefContext = new ReferenceContext(ref, pileupInterval);
            final ActivityProfileState isActiveResult = hcEngine.isActive(pileup, pileupRefContext, new FeatureContext(null, pileupInterval));
            final double expectedIsActiveValue = expectedActiveSites.contains(pileupInterval) ? 1.0 : 0.0;
            Assert.assertEquals(isActiveResult.isActiveProb(), expectedIsActiveValue, "Wrong isActive probability for site " + pileupInterval);
        });
    }
}
Also used : ReadFilteringIterator(org.broadinstitute.hellbender.utils.iterators.ReadFilteringIterator) GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) ActivityProfileState(org.broadinstitute.hellbender.utils.activityprofile.ActivityProfileState) CachingIndexedFastaSequenceFile(org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile) LocusIteratorByState(org.broadinstitute.hellbender.utils.locusiterator.LocusIteratorByState) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) ReadFilter(org.broadinstitute.hellbender.engine.filters.ReadFilter) CachingIndexedFastaSequenceFile(org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile) File(java.io.File) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 14 with CachingIndexedFastaSequenceFile

use of org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile in project gatk by broadinstitute.

the class AssemblyRegionUnitTest method init.

@BeforeClass
public void init() throws FileNotFoundException {
    // sequence
    seq = new CachingIndexedFastaSequenceFile(new File(hg19MiniReference));
    contig = "1";
    contigLength = seq.getSequence(contig).length();
    header = ArtificialReadUtils.createArtificialSamHeader(seq.getSequenceDictionary());
}
Also used : CachingIndexedFastaSequenceFile(org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile) CachingIndexedFastaSequenceFile(org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile) File(java.io.File) IndexedFastaSequenceFile(htsjdk.samtools.reference.IndexedFastaSequenceFile) BeforeClass(org.testng.annotations.BeforeClass)

Example 15 with CachingIndexedFastaSequenceFile

use of org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile in project gatk by broadinstitute.

the class BaseTest method initGenomeLocParser.

@BeforeClass
public void initGenomeLocParser() throws FileNotFoundException {
    hg19ReferenceReader = new CachingIndexedFastaSequenceFile(new File(hg19MiniReference));
    hg19Header = new SAMFileHeader();
    hg19Header.setSequenceDictionary(hg19ReferenceReader.getSequenceDictionary());
    hg19GenomeLocParser = new GenomeLocParser(hg19ReferenceReader);
}
Also used : CachingIndexedFastaSequenceFile(org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile) SAMFileHeader(htsjdk.samtools.SAMFileHeader) GenomeLocParser(org.broadinstitute.hellbender.utils.GenomeLocParser) CachingIndexedFastaSequenceFile(org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile) File(java.io.File) BeforeClass(org.testng.annotations.BeforeClass)

Aggregations

CachingIndexedFastaSequenceFile (org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile)15 File (java.io.File)14 BeforeClass (org.testng.annotations.BeforeClass)6 SAMFileHeader (htsjdk.samtools.SAMFileHeader)5 IndexedFastaSequenceFile (htsjdk.samtools.reference.IndexedFastaSequenceFile)5 BaseTest (org.broadinstitute.hellbender.utils.test.BaseTest)5 Test (org.testng.annotations.Test)5 ReferenceSequenceFile (htsjdk.samtools.reference.ReferenceSequenceFile)4 FileNotFoundException (java.io.FileNotFoundException)3 UserException (org.broadinstitute.hellbender.exceptions.UserException)3 GenomeLocParser (org.broadinstitute.hellbender.utils.GenomeLocParser)3 ReadFilter (org.broadinstitute.hellbender.engine.filters.ReadFilter)2 SimpleInterval (org.broadinstitute.hellbender.utils.SimpleInterval)2 ActivityProfileState (org.broadinstitute.hellbender.utils.activityprofile.ActivityProfileState)2 ReadFilteringIterator (org.broadinstitute.hellbender.utils.iterators.ReadFilteringIterator)2 LocusIteratorByState (org.broadinstitute.hellbender.utils.locusiterator.LocusIteratorByState)2 GATKRead (org.broadinstitute.hellbender.utils.read.GATKRead)2 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)1 SAMSequenceRecord (htsjdk.samtools.SAMSequenceRecord)1 NDNCigarReadTransformer (org.broadinstitute.hellbender.transformers.NDNCigarReadTransformer)1