Java Code Examples for htsjdk.samtools.SAMSequenceDictionary#assertSameDictionary()

/**
 * Interprets intervals from the input INTERVALS, if there's a file with that name, it opens the file, otherwise
 * it tries to parse it, checks that their dictionaries all agree with the input SequenceDictionary and returns the
 * intersection of all the lists.
 *
 * @param sequenceDictionary a {@link SAMSequenceDictionary} to parse the intervals against.
 * @return the intersection of the interval lists pointed to by the input parameter INTERVALS
 */
private IntervalList getIntervals(final SAMSequenceDictionary sequenceDictionary) {

    IntervalList regionOfInterest = null;

    for (final File intervalListFile : INTERVALS) {

        if (!intervalListFile.exists()) {
            throw new IllegalArgumentException("Input file " + intervalListFile + " doesn't seem to exist. ");
        }

        log.info("Reading IntervalList ", intervalListFile, ".");
        final IntervalList temp = IntervalList.fromFile(intervalListFile);
        sequenceDictionary.assertSameDictionary(temp.getHeader().getSequenceDictionary());

        if (regionOfInterest == null) {
            regionOfInterest = temp;
        } else {
            log.info("Intersecting interval_list: ", intervalListFile, ".");
            regionOfInterest = IntervalList.intersection(regionOfInterest, temp);
        }
    }
    return regionOfInterest;
}
 

@Test
public void testBestSequenceDictionary_FromVariantReference() throws Exception {
    final GATKTool tool = new TestGATKToolWithFeatures();
    final CommandLineParser clp = new CommandLineArgumentParser(tool);
    final File vcfFile = new File(publicTestDir + "org/broadinstitute/hellbender/engine/example_variants_noSequenceDict.vcf");
    final String[] args = {"-V", vcfFile.getCanonicalPath(), "-R", hg19MiniReference};
    clp.parseArguments(System.out, args);
    tool.onStartup();
    // make sure we DON'T get the seq dictionary from the VCF index, and instead get the one from
    // the reference when its better
    final SAMSequenceDictionary toolDict = tool.getBestAvailableSequenceDictionary();
    Assert.assertFalse(toolDict.getSequences().stream().allMatch(seqRec -> seqRec.getSequenceLength() == 0));

    SAMSequenceDictionary refDict = new ReferenceFileSource(IOUtils.getPath(hg19MiniReference)).getSequenceDictionary();
    toolDict.assertSameDictionary(refDict);
    refDict.assertSameDictionary(toolDict);
    Assert.assertEquals(toolDict, refDict);
}
 

private void collectFileReadersAndHeaders(final List<String> sampleList, SAMSequenceDictionary samSequenceDictionary) {
    for (final File input : INPUT) {
        final VCFFileReader in = new VCFFileReader(input, false);
        final VCFHeader header = in.getFileHeader();
        final SAMSequenceDictionary dict = in.getFileHeader().getSequenceDictionary();
        if (dict == null || dict.isEmpty()) {
            if (null == samSequenceDictionary) {
                throw new IllegalArgumentException("Sequence dictionary was missing or empty for the VCF: " + input.getAbsolutePath() + " Please add a sequence dictionary to this VCF or specify SEQUENCE_DICTIONARY.");
            }
            header.setSequenceDictionary(samSequenceDictionary);
        } else {
            if (null == samSequenceDictionary) {
                samSequenceDictionary = dict;
            } else {
                try {
                    samSequenceDictionary.assertSameDictionary(dict);
                } catch (final AssertionError e) {
                    throw new IllegalArgumentException(e);
                }
            }
        }
        if (sampleList.isEmpty()) {
            sampleList.addAll(header.getSampleNamesInOrder());
        } else {
            if (!sampleList.equals(header.getSampleNamesInOrder())) {
                throw new IllegalArgumentException("Input file " + input.getAbsolutePath() + " has sample names that don't match the other files.");
            }
        }
        inputReaders.add(in);
        inputHeaders.add(header);
    }
}
 

/** Validates that all headers contain the same set of genotyped samples and that files are in order by position of first record. */
private static void assertSameSamplesAndValidOrdering(final List<Path> inputFiles, final boolean disableContigOrderingCheck) {
    final VCFHeader firstHeader = getHeader(inputFiles.get(0));
    final SAMSequenceDictionary dict = firstHeader.getSequenceDictionary();
    if ( dict == null) {
        throw new UserException.BadInput("The first VCF specified is missing the required sequence dictionary. " +
                                                 "This is required to perform validation.  You can skip this validation " +
                                                 "using --"+IGNORE_SAFETY_CHECKS_LONG_NAME +" but ignoring safety checks " +
                                                 "can result in invalid output.");
    }
    final VariantContextComparator comparator = new VariantContextComparator(dict);
    final List<String> samples = firstHeader.getGenotypeSamples();

    Path lastFile = null;
    VariantContext lastContext = null;

    for (final Path f : inputFiles) {
        final FeatureReader<VariantContext> in = getReaderFromVCFUri(f, 0);
        final VCFHeader header = (VCFHeader)in.getHeader();
        dict.assertSameDictionary(header.getSequenceDictionary());
        final List<String> theseSamples = header.getGenotypeSamples();

        if (!samples.equals(theseSamples)) {
            final SortedSet<String> s1 = new TreeSet<>(samples);
            final SortedSet<String> s2 = new TreeSet<>(theseSamples);
            s1.removeAll(theseSamples);
            s2.removeAll(samples);

            throw new IllegalArgumentException("VCFs do not have identical sample lists." +
                    " Samples unique to first file: " + s1 + ". Samples unique to " + f.toUri().toString() + ": " + s2 + ".");
        }

        try(final CloseableIterator<VariantContext> variantIterator = in.iterator()) {
            if (variantIterator.hasNext()) {
                final VariantContext currentContext = variantIterator.next();
                if (lastContext != null) {
                    if ( disableContigOrderingCheck ) {
                        if ( lastContext.getContig().equals(currentContext.getContig()) && lastContext.getStart() >= currentContext.getStart() ) {
                            throw new IllegalArgumentException(
                                    "First record in file " + f.toUri().toString() + " is not after first record in " +
                                            "previous file " + lastFile.toUri().toString());
                        }
                    }
                    else {
                        if ( comparator.compare(lastContext, currentContext) >= 0 ) {
                            throw new IllegalArgumentException(
                                    "First record in file " + f.toUri().toString() + " is not after first record in " +
                                            "previous file " + lastFile.toUri().toString());
                        }
                    }
                }

                lastContext = currentContext;
                lastFile = f;
            }
        } catch (final IOException e) {
            throw new UserException.CouldNotReadInputFile(f, e.getMessage(), e);
        }

        CloserUtil.close(in);
    }
}