htsjdk.tribble.CloseableTribbleIterator Java Examples

/**
 * Refill our cache from disk after a cache miss. Will prefetch Features overlapping an additional
 * queryLookaheadBases bases after the end of the provided interval, in addition to those overlapping
 * the interval itself.
 * <p>
 * Calling this has the side effect of invalidating (closing) any currently-open iteration over
 * this data source.
 *
 * @param interval the query interval that produced a cache miss
 */
private void refillQueryCache(final Locatable interval) {
    // Tribble documentation states that having multiple iterators open simultaneously over the same FeatureReader
    // results in undefined behavior
    closeOpenIterationIfNecessary();

    // Expand the end of our query by the configured number of bases, in anticipation of probable future
    // queries with slightly larger start/stop positions.
    //
    // Note that it doesn't matter if we go off the end of the contig in the process, since
    // our reader's query operation is not aware of (and does not care about) contig boundaries.
    // Note: we use addExact to blow up on overflow rather than propagate negative results downstream
    final SimpleInterval queryInterval = new SimpleInterval(interval.getContig(), interval.getStart(), Math.addExact(interval.getEnd(), queryLookaheadBases));

    // Query iterator over our reader will be immediately closed after re-populating our cache
    try (final CloseableTribbleIterator<T> queryIter = featureReader.query(queryInterval.getContig(), queryInterval.getStart(), queryInterval.getEnd())) {
        queryCache.fill(queryIter, queryInterval);
    } catch (final IOException e) {
        throw new GATKException("Error querying file " + featureInput + " over interval " + interval, e);
    }
}
 

/** Create a collection based on the contents of an input file and a given config file.  The config file must be the same as
 * is ingested by {@link XsvLocatableTableCodec}.
 *
 * @param input readable path to use for the xsv file.  Must be readable.  Never {@code null}.
 * @param inputConfigFile config file for specifying the format of the xsv file.  Must be readable.  Never {@code null}.
 * @param headersOfInterest Only preserve these headers.  Only a warning will be logged if all of these are not
 *                          present in the input file.  This parameter should not include the locatable columns
 *                          defined by the config file, which are always preserved.
 *                          Use {@code null} to indicate "all headers are of interest".
 * @return never {@code null}
 */
public static AnnotatedIntervalCollection create(final Path input, final Path inputConfigFile, final Set<String> headersOfInterest) {

    IOUtils.assertFileIsReadable(input);
    IOUtils.assertFileIsReadable(inputConfigFile);

    final AnnotatedIntervalCodec codec = new AnnotatedIntervalCodec(inputConfigFile);
    final List<AnnotatedInterval> regions = new ArrayList<>();
    
    if (codec.canDecode(input.toUri().toString())) {
        try (final FeatureReader<AnnotatedInterval> reader = AbstractFeatureReader.getFeatureReader(input.toUri().toString(), codec, false)){

            // This cast is an artifact of the tribble framework.
            @SuppressWarnings("unchecked")
            final AnnotatedIntervalHeader header = (AnnotatedIntervalHeader) reader.getHeader();

            final List<String> finalAnnotations = determineCollectionAnnotations(headersOfInterest, header.getAnnotations());

            final CloseableTribbleIterator<AnnotatedInterval> it = reader.iterator();
            StreamSupport.stream(it.spliterator(), false)
                    .filter(r -> r != null)
                    .map(r -> copyAnnotatedInterval(r,  new HashSet<>(finalAnnotations)))
                    .forEach(r -> regions.add(r));

            return new AnnotatedIntervalCollection(header.getSamFileHeader(), finalAnnotations, regions);

        } catch ( final IOException ex ) {
            throw new GATKException("Error - IO problem with file " + input, ex);
        }
    } else {
        throw new UserException.BadInput("Could not parse xsv file: " + input.toUri().toString());
    }
}
 

@Override
public CloseableTribbleIterator<VariantContext> query(final String chr, final int start, final int end) {
    final SimpleInterval queryInterval = new SimpleInterval(chr, start, end);
    if( !interval.equals(queryInterval)){
        throw new GATKException("Cannot call query with different interval, expected:" + this.interval + " queried with: " + queryInterval);
    }
    if( query != null ){
        final CloseableTribbleIterator<VariantContext> tmp = query;
        query = null;
        return tmp;
    } else {
        throw new GATKException("Cannot call query twice on this wrapper.");
    }
}
 

private static List<GencodeFuncotation> createGencodeFuncotations(final String contig, final int start, final int end, final String ref, final String alt, final String referenceFileName, final ReferenceDataSource referenceDataSource, final FeatureReader<GencodeGtfFeature> featureReader, final String transcriptFastaFile, final String transcriptGtfFile, final TranscriptSelectionMode transcriptSelectionMode) {
    final SimpleInterval variantInterval = new SimpleInterval( contig, start, end );
    final VariantContext variantContext = createVariantContext(contig, start, end, ref, alt, referenceFileName);

    final ReferenceContext referenceContext = new ReferenceContext(referenceDataSource, variantInterval );

    // Get our gene feature iterator:
    final CloseableTribbleIterator<GencodeGtfFeature> gtfFeatureIterator;
    try {
        gtfFeatureIterator = featureReader.query(variantContext.getContig(), variantContext.getStart(), variantContext.getEnd());
    }
    catch (final IOException ex) {
        throw new GATKException("Could not finish the test!", ex);
    }
    final List<Feature> featureList = Collections.singletonList(gtfFeatureIterator.next());

    final String gencode_test = "GENCODE_TEST";
    final GencodeFuncotationFactory gencodeFactory = new GencodeFuncotationFactory(Paths.get(transcriptFastaFile),
    "TEST", gencode_test, transcriptSelectionMode, new HashSet<>(), new LinkedHashMap<>(),
            new FeatureInput<>(transcriptGtfFile, gencode_test, Collections.emptyMap()), "TEST");

    final FeatureContext featureContext = FuncotatorTestUtils.createFeatureContext(Collections.singletonList(gencodeFactory), "FuncotationMapUnitTest",
            variantInterval, 0, 0, 0, null);

    return gencodeFactory.createFuncotations(variantContext, referenceContext, featureContext).stream()
        .map(f -> (GencodeFuncotation) f).collect(Collectors.toList());
}
 

private static GencodeGtfExonFeature getExonForVariant( final CloseableTribbleIterator<GencodeGtfFeature> gtfFeatureIterator,
                                                        final SimpleInterval variantLocation ) {

    final Optional<GencodeGtfExonFeature> exonOption = gtfFeatureIterator.stream()
                                                .map( f -> ((GencodeGtfGeneFeature) f))
                                                .flatMap(g -> g.getTranscripts().stream())
                                                .flatMap(t -> t.getExons().stream())
                                                .filter( x -> x.overlaps(variantLocation) )
                                                .findFirst();

    return exonOption.orElseThrow( () -> new GATKException("Could not get an exon associated with this variant's interval: " + variantLocation.toString()) );
}
 

@Test(dataProvider = "getRenameCombinations")
public void testRenamingSamples(final Map<String, String> renamingMap, final int threads, final int batchSize) throws IOException {
    final LinkedHashMap<String, String> sampleMap = new LinkedHashMap<>(renamingMap);
    sampleMap.replaceAll( (newSampleName, originalSampleName)-> createInputVCF(originalSampleName).getAbsolutePath());

    final File sampleMapFile = getSampleMapFile(sampleMap);

    final String workspace = createTempDir("workspace").getAbsolutePath();
    Files.delete(Paths.get(workspace));
    final ArgumentsBuilder args = new ArgumentsBuilder()
            .add(GenomicsDBImport.SAMPLE_NAME_MAP_LONG_NAME, sampleMapFile.getAbsolutePath())
            .add(GenomicsDBImport.WORKSPACE_ARG_LONG_NAME, new File(workspace).getAbsolutePath())
            .add(GenomicsDBImport.VCF_INITIALIZER_THREADS_LONG_NAME, String.valueOf(threads))
            .add(GenomicsDBImport.BATCHSIZE_ARG_LONG_NAME, String.valueOf(batchSize))
            .addInterval(INTERVAL.get(0));

    runCommandLine(args);
    final Set<String> expectedSampleNames = sampleMap.keySet();
    try(final FeatureReader<VariantContext> reader = getGenomicsDBFeatureReader(workspace, b37_reference_20_21)) {
        final CloseableTribbleIterator<VariantContext> iterator = reader.iterator();
        Assert.assertTrue(iterator.hasNext(), "expected to see a variant");
        Assert.assertTrue(expectedSampleNames.size() > 0);
        Assert.assertEquals(expectedSampleNames.size(), renamingMap.size());
        iterator.forEachRemaining(vc -> {
            Assert.assertEquals(vc.getSampleNames().size(), expectedSampleNames.size());
            Assert.assertEqualsNoOrder(vc.getSampleNames().toArray(), expectedSampleNames.toArray());
            expectedSampleNames.forEach( sample -> {
                Assert.assertEquals(vc.getGenotype(sample).getAnyAttribute(SAMPLE_NAME_KEY), renamingMap.get(sample));
                //check another attribute just to make sure we're not mangling things
                Assert.assertEquals(VariantContextGetters.getAttributeAsInt(vc.getGenotype(sample), ANOTHER_ATTRIBUTE_KEY, -1), 10);
            });
        });
    }

}
 

public LongISLNDReadAlignmentMap(final Collection<String> readAlignmentMapFiles) throws IOException {
    for (final String readAlignmentMapFileName : readAlignmentMapFiles) {
        log.info("Reading in read map from " + readAlignmentMapFileName);
        final AbstractFeatureReader<BEDFeature, LineIterator> featureReader = AbstractFeatureReader.getFeatureReader(readAlignmentMapFileName, new BEDCodec(), false);
        try {
            final CloseableTribbleIterator<BEDFeature> featureIterator = featureReader.iterator();
            while (featureIterator.hasNext()) {
                final BEDFeature feature = featureIterator.next();
                readAlignmentMap.put(feature.getName(), new LongISLNDReadMapRecord(feature).toReadMapRecord());
            }
        } finally {
            featureReader.close();
        }
    }
}
 

@Override
protected int doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsReadable(SEQUENCE_DICTIONARY);
    IOUtil.assertFileIsWritable(OUTPUT);
    try {
        // create a new header that we will assign the dictionary provided by the SAMSequenceDictionaryExtractor to.
        final SAMFileHeader header = new SAMFileHeader();
        final SAMSequenceDictionary samSequenceDictionary = SAMSequenceDictionaryExtractor.extractDictionary(SEQUENCE_DICTIONARY.toPath());
        header.setSequenceDictionary(samSequenceDictionary);
        // set the sort order to be sorted by coordinate, which is actually done below
        // by getting the .uniqued() intervals list before we write out the file
        header.setSortOrder(SAMFileHeader.SortOrder.coordinate);
        final IntervalList intervalList = new IntervalList(header);

        final FeatureReader<BEDFeature> bedReader = AbstractFeatureReader.getFeatureReader(INPUT.getAbsolutePath(), new BEDCodec(), false);
        final CloseableTribbleIterator<BEDFeature> iterator = bedReader.iterator();
        final ProgressLogger progressLogger = new ProgressLogger(LOG, (int) 1e6);

        while (iterator.hasNext()) {
            final BEDFeature bedFeature = iterator.next();
            final String sequenceName = bedFeature.getContig();
            final int start = bedFeature.getStart();
            final int end = bedFeature.getEnd();
            // NB: do not use an empty name within an interval
            final String name;
            if (bedFeature.getName().isEmpty()) {
                name = null;
            } else {
                name = bedFeature.getName();
            }

            final SAMSequenceRecord sequenceRecord = header.getSequenceDictionary().getSequence(sequenceName);

            // Do some validation
            if (null == sequenceRecord) {
                if (DROP_MISSING_CONTIGS) {
                    LOG.info(String.format("Dropping interval with missing contig: %s:%d-%d", sequenceName, bedFeature.getStart(), bedFeature.getEnd()));
                    missingIntervals++;
                    missingRegion += bedFeature.getEnd() - bedFeature.getStart();
                    continue;
                }
                throw new PicardException(String.format("Sequence '%s' was not found in the sequence dictionary", sequenceName));
            } else if (start < 1) {
                throw new PicardException(String.format("Start on sequence '%s' was less than one: %d", sequenceName, start));
            } else if (sequenceRecord.getSequenceLength() < start) {
                throw new PicardException(String.format("Start on sequence '%s' was past the end: %d < %d", sequenceName, sequenceRecord.getSequenceLength(), start));
            } else if ((end == 0 && start != 1 ) //special case for 0-length interval at the start of a contig
                    || end < 0 ) {
                throw new PicardException(String.format("End on sequence '%s' was less than one: %d", sequenceName, end));
            } else if (sequenceRecord.getSequenceLength() < end) {
                throw new PicardException(String.format("End on sequence '%s' was past the end: %d < %d", sequenceName, sequenceRecord.getSequenceLength(), end));
            } else if (end < start - 1) {
                throw new PicardException(String.format("On sequence '%s', end < start-1: %d <= %d", sequenceName, end, start));
            }

            final boolean isNegativeStrand = bedFeature.getStrand() == Strand.NEGATIVE;
            final Interval interval = new Interval(sequenceName, start, end, isNegativeStrand, name);
            intervalList.add(interval);

            progressLogger.record(sequenceName, start);
        }
        CloserUtil.close(bedReader);

        if (DROP_MISSING_CONTIGS) {
            if (missingRegion == 0) {
                LOG.info("There were no missing regions.");
            } else {
                LOG.warn(String.format("There were %d missing regions with a total of %d bases", missingIntervals, missingRegion));
            }
        }
        // Sort and write the output
        IntervalList out = intervalList;
        if (SORT) {
            out = out.sorted();
        }
        if (UNIQUE) {
            out = out.uniqued();
        }
        out.write(OUTPUT);
        LOG.info(String.format("Wrote %d intervals spanning a total of %d bases", out.getIntervals().size(),out.getBaseCount()));

    } catch (final IOException e) {
        throw new RuntimeException(e);
    }

    return 0;
}
 

@Override
public CloseableTribbleIterator<VariantContext> iterator() {
    throw new UnsupportedOperationException("iterator() not supported, this should not have been called and indicates an issue with GenomicsDB integration");
}
 

@Test (dataProvider = "provideMuc16SnpDataForGetVariantClassification")
void testGetVariantClassificationForCodingRegions(final int chromosomeNumber,
                                  final int start,
                                  final int end,
                                  final GencodeFuncotation.VariantType variantType,
                                  final String ref,
                                  final String alt,
                                  final GencodeFuncotation.VariantClassification expectedVariantClassification) {

    // This test can only deal with variants in coding regions.
    // So we ignore any tests that are expected outside of coding regions.
    // i.e. expectedVariantClassification is one of:
    //     { INTRON, FIVE_PRIME_UTR, THREE_PRIME_UTR, IGR, FIVE_PRIME_FLANK, DE_NOVO_START_IN_FRAME, DE_NOVO_START_OUT_FRAME, RNA, LINCRNA }
    // We test these cases in another unit test.
    if ((expectedVariantClassification == GencodeFuncotation.VariantClassification.INTRON) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.FIVE_PRIME_UTR) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.THREE_PRIME_UTR) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.IGR) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.FIVE_PRIME_FLANK) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.DE_NOVO_START_IN_FRAME) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.DE_NOVO_START_OUT_FRAME) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.RNA) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.LINCRNA) )
    {
        return;
    }

    final String contig = "chr" + Integer.toString(chromosomeNumber);
    final SimpleInterval variantInterval = new SimpleInterval( contig, start, end );

    final Allele refAllele = Allele.create(ref, true);
    final Allele altAllele = Allele.create(alt);

    final VariantContextBuilder variantContextBuilder = new VariantContextBuilder(
            FuncotatorReferenceTestUtils.retrieveHg19Chr19Ref(),
            contig,
            start,
            end,
            Arrays.asList(refAllele, altAllele)
    );
    final VariantContext variantContext = variantContextBuilder.make();

    // Get our gene feature iterator:
    final CloseableTribbleIterator<GencodeGtfFeature> gtfFeatureIterator;
    try {
        gtfFeatureIterator = gencodeHg19FeatureReader.query(contig, start, end);
    }
    catch (final IOException ex) {
        throw new GATKException("Could not finish the test!", ex);
    }

    // Get the gene.
    // We know the first gene is the right one - the gene in question is the MUC16 gene:
    final GencodeGtfGeneFeature             gene = (GencodeGtfGeneFeature) gtfFeatureIterator.next();
    final GencodeGtfTranscriptFeature transcript = getMuc16Transcript(gene);
    final GencodeGtfExonFeature             exon = getExonForVariant( gene, variantInterval );

    final ReferenceContext referenceContext = new ReferenceContext( refDataSourceHg19Ch19, variantInterval );

    final List<? extends Locatable> exonPositionList = GencodeFuncotationFactory.getSortedCdsAndStartStopPositions(transcript);

    final ReferenceDataSource muc16TranscriptDataSource = ReferenceDataSource.of(new File(FuncotatorTestConstants.GENCODE_DATA_SOURCE_FASTA_PATH_HG19).toPath());
    final Map<String, GencodeFuncotationFactory.MappedTranscriptIdInfo> muc16TranscriptIdMap = GencodeFuncotationFactory. createTranscriptIdMap(muc16TranscriptDataSource);

    final SequenceComparison seqComp =
            GencodeFuncotationFactory.createSequenceComparison(
                    variantContext,
                    altAllele,
                    referenceContext,
                    transcript,
                    exonPositionList,
                    muc16TranscriptIdMap,
                    muc16TranscriptDataSource,
                    true);

    final GencodeFuncotation.VariantClassification varClass = GencodeFuncotationFactory.createVariantClassification(
            variantContext,
            altAllele,
            variantType,
            exon,
            transcript.getExons().size(),
            seqComp
    );

    Assert.assertEquals( varClass, expectedVariantClassification );
}
 

@Test (dataProvider = "provideMuc16SnpDataForGetVariantClassificationWithOutOfCdsData")
void testMuc16SnpCreateFuncotations(final int chromosomeNumber,
                                    final int start,
                                    final int end,
                                    final GencodeFuncotation.VariantType expectedVariantType,
                                    final String ref,
                                    final String alt,
                                    final GencodeFuncotation.VariantClassification expectedVariantClassification) {

    final String contig = "chr" + Integer.toString(chromosomeNumber);
    final SimpleInterval variantInterval = new SimpleInterval( contig, start, end );

    final Allele refAllele = Allele.create(ref, true);
    final Allele altAllele = Allele.create(alt);

    final VariantContextBuilder variantContextBuilder = new VariantContextBuilder(
            FuncotatorReferenceTestUtils.retrieveHg19Chr19Ref(),
            contig,
            start,
            end,
            Arrays.asList(refAllele, altAllele)
    );
    final VariantContext variantContext = variantContextBuilder.make();

    // Get our gene feature iterator:
    final CloseableTribbleIterator<GencodeGtfFeature> gtfFeatureIterator;
    try {
        gtfFeatureIterator = gencodeHg19FeatureReader.query(contig, start, end);
    }
    catch (final IOException ex) {
        throw new GATKException("Could not finish the test!", ex);
    }

    // Get the gene.
    // We know the first gene is the right one - the gene in question is the MUC16 gene:
    final GencodeGtfGeneFeature gene = (GencodeGtfGeneFeature) gtfFeatureIterator.next();
    final ReferenceContext referenceContext = new ReferenceContext(refDataSourceHg19Ch19, variantInterval );
    final Set<String> requestedTranscriptIds = getValidTranscriptsForGene("MUC16");

    // Create a factory for our funcotations:
    try (final GencodeFuncotationFactory funcotationFactory = new GencodeFuncotationFactory(
            IOUtils.getPath(FuncotatorTestConstants.GENCODE_DATA_SOURCE_FASTA_PATH_HG19),
            "VERSION",
            GencodeFuncotationFactory.DEFAULT_NAME,
            FuncotatorArgumentDefinitions.TRANSCRIPT_SELECTION_MODE_DEFAULT_VALUE,
            requestedTranscriptIds,
            new LinkedHashMap<>(), createFeatureInputForMuc16Ds(GencodeFuncotationFactory.DEFAULT_NAME), "HG19")) {

        // Generate our funcotations:
        final List<Feature> featureList = new ArrayList<>();
        featureList.add( gene );
        final List<Funcotation> funcotations = funcotationFactory.createFuncotationsOnVariant(variantContext, referenceContext, featureList);

        // Make sure we get what we expected:
        Assert.assertEquals(funcotations.size(), 1);

        final GencodeFuncotation funcotation = (GencodeFuncotation)funcotations.get(0);

        Assert.assertEquals(funcotation.getVariantClassification(), expectedVariantClassification);
        Assert.assertEquals(funcotation.getVariantType(), expectedVariantType);
    }
}
 

@Test(dataProvider = "provideForCreateNonBasicFuncotations")
void createNonBasicFuncotations(final int start, final int end) {

    final String contig = "chr19";

    final SimpleInterval variantInterval = new SimpleInterval( contig, start, end );

    final Allele refAllele = Allele.create("C", true);
    final Allele altAllele = Allele.create("G", false);

    final VariantContextBuilder variantContextBuilder = new VariantContextBuilder(
            FuncotatorReferenceTestUtils.retrieveHg19Chr19Ref(),
            contig,
            start,
            end,
            Arrays.asList(refAllele, altAllele)
    );
    final VariantContext variantContext = variantContextBuilder.make();

    // Get our gene feature iterator:
    final CloseableTribbleIterator<GencodeGtfFeature> gtfFeatureIterator;
    try {
        gtfFeatureIterator = muc16NonBasicFeatureReader.query(contig, start, end);
    }
    catch (final IOException ex) {
        throw new GATKException("Could not finish the test!", ex);
    }

    // Get the gene.
    // We know the first gene is the right one - the gene in question is the MUC16 gene:
    final GencodeGtfGeneFeature gene = (GencodeGtfGeneFeature) gtfFeatureIterator.next();
    final ReferenceContext referenceContext = new ReferenceContext(refDataSourceHg19Ch19, variantInterval );

    // Create a factory for our funcotations:
    try (final GencodeFuncotationFactory funcotationFactory = new GencodeFuncotationFactory(
                                                                IOUtils.getPath(FuncotatorTestConstants.GENCODE_DATA_SOURCE_FASTA_PATH_HG19),
                                                                "VERSION",
                                                                GencodeFuncotationFactory.DEFAULT_NAME,
                                                                FuncotatorArgumentDefinitions.TRANSCRIPT_SELECTION_MODE_DEFAULT_VALUE,
                                                                new HashSet<>(),
                                                                new LinkedHashMap<>(),
                                                                createFeatureInputForMuc16Ds(GencodeFuncotationFactory.DEFAULT_NAME),
                                                                "HG19")) {

        // Generate our funcotations:
        final List<Feature> featureList = new ArrayList<>();
        featureList.add( gene );
        final List<Funcotation> funcotations = funcotationFactory.createFuncotationsOnVariant(variantContext, referenceContext, featureList);

        // Make sure we get what we expected (a single IGR funcotation):
        Assert.assertEquals(funcotations.size(), 1);
        Assert.assertTrue(funcotations.get(0) instanceof GencodeFuncotation);

        final GencodeFuncotation gencodeFuncotation = (GencodeFuncotation)funcotations.get(0);
        Assert.assertEquals(gencodeFuncotation.getVariantClassification(), GencodeFuncotation.VariantClassification.IGR);
        Assert.assertNull(gencodeFuncotation.getHugoSymbol());
    }
}