Java Code Examples for htsjdk.variant.vcf.VCFHeader#setSequenceDictionary()

private static VariantContextWriter getVariantContextWriter(final File outputFile,
                                                            final File referenceSequenceFileName,
                                                            final String sample,
                                                            final String source,
                                                            final ReferenceSequenceFile ref) {
    final VariantContextWriter variantContextWriter = new VariantContextWriterBuilder()
            .setReferenceDictionary(ref.getSequenceDictionary())
            .setOutputFile(outputFile).build();

    final Set<VCFHeaderLine> lines = new LinkedHashSet<>();
    lines.add(new VCFHeaderLine("reference", referenceSequenceFileName.getAbsolutePath()));
    lines.add(new VCFHeaderLine("source", source));
    lines.add(new VCFHeaderLine("fileDate", new Date().toString()));

    lines.add(VCFStandardHeaderLines.getFormatLine(VCFConstants.GENOTYPE_PL_KEY));
    lines.add(VCFStandardHeaderLines.getFormatLine(VCFConstants.GENOTYPE_ALLELE_DEPTHS));
    lines.add(VCFStandardHeaderLines.getFormatLine(VCFConstants.DEPTH_KEY));

    final VCFHeader header = new VCFHeader(lines, Collections.singletonList(sample));
    header.setSequenceDictionary(ref.getSequenceDictionary());
    variantContextWriter.writeHeader(header);
    return variantContextWriter;
}
 

public SageVCF(@NotNull final IndexedFastaSequenceFile reference, @NotNull final SageConfig config) {
    writer = new VariantContextWriterBuilder().setOutputFile(config.outputFile())
            .modifyOption(Options.INDEX_ON_THE_FLY, true)
            .modifyOption(Options.USE_ASYNC_IO, false)
            .setReferenceDictionary(reference.getSequenceDictionary())
            .build();
    refContextEnrichment = new SomaticRefContextEnrichment(reference, this::writeToFile);

    final VCFHeader header = refContextEnrichment.enrichHeader(header(config));
    header.setSequenceDictionary(reference.getSequenceDictionary());
    writer.writeHeader(header);
}
 

/**
 * Writes an appropriate VCF header, given our arguments, to the provided writer
 *
 * @param vcfWriter writer to which the header should be written
 */
public void writeHeader( final VariantContextWriter vcfWriter, final SAMSequenceDictionary sequenceDictionary,
                         final Set<VCFHeaderLine>  defaultToolHeaderLines) {
    Utils.nonNull(vcfWriter);

    final Set<VCFHeaderLine> headerInfo = new HashSet<>();
    headerInfo.addAll(defaultToolHeaderLines);

    headerInfo.addAll(genotypingEngine.getAppropriateVCFInfoHeaders());
    // all annotation fields from VariantAnnotatorEngine
    headerInfo.addAll(annotationEngine.getVCFAnnotationDescriptions());
    // all callers need to add these standard annotation header lines
    headerInfo.add(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.DOWNSAMPLED_KEY));
    headerInfo.add(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.MLE_ALLELE_COUNT_KEY));
    headerInfo.add(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.MLE_ALLELE_FREQUENCY_KEY));
    // all callers need to add these standard FORMAT field header lines
    VCFStandardHeaderLines.addStandardFormatLines(headerInfo, true,
            VCFConstants.GENOTYPE_KEY,
            VCFConstants.GENOTYPE_QUALITY_KEY,
            VCFConstants.DEPTH_KEY,
            VCFConstants.GENOTYPE_PL_KEY);

    if ( ! hcArgs.doNotRunPhysicalPhasing ) {
        headerInfo.add(GATKVCFHeaderLines.getFormatLine(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_ID_KEY));
        headerInfo.add(GATKVCFHeaderLines.getFormatLine(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_GT_KEY));
    }

    // FILTER fields are added unconditionally as it's not always 100% certain the circumstances
    // where the filters are used.  For example, in emitting all sites the lowQual field is used
    headerInfo.add(GATKVCFHeaderLines.getFilterLine(GATKVCFConstants.LOW_QUAL_FILTER_NAME));

    if ( emitReferenceConfidence() ) {
        headerInfo.addAll(referenceConfidenceModel.getVCFHeaderLines());
    }

    final VCFHeader vcfHeader = new VCFHeader(headerInfo, sampleSet);
    vcfHeader.setSequenceDictionary(sequenceDictionary);
    vcfWriter.writeHeader(vcfHeader);
}
 

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);
    }
}
 

@Override
protected int doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsReadable(SEQUENCE_DICTIONARY);
    IOUtil.assertFileIsWritable(OUTPUT);

    final SAMSequenceDictionary samSequenceDictionary = SAMSequenceDictionaryExtractor.extractDictionary(SEQUENCE_DICTIONARY.toPath());

    final VCFFileReader fileReader = new VCFFileReader(INPUT, false);
    final VCFHeader fileHeader = fileReader.getFileHeader();

    final VariantContextWriterBuilder builder = new VariantContextWriterBuilder()
            .setReferenceDictionary(samSequenceDictionary)
            .clearOptions();
    if (CREATE_INDEX)
        builder.setOption(Options.INDEX_ON_THE_FLY);

    final VariantContextWriter vcfWriter = builder.setOutputFile(OUTPUT).build();
    fileHeader.setSequenceDictionary(samSequenceDictionary);
    vcfWriter.writeHeader(fileHeader);

    final ProgressLogger progress = new ProgressLogger(log, 10000);
    final CloseableIterator<VariantContext> iterator = fileReader.iterator();
    while (iterator.hasNext()) {
        final VariantContext context = iterator.next();
        vcfWriter.add(context);
        progress.record(context.getContig(), context.getStart());
    }

    CloserUtil.close(iterator);
    CloserUtil.close(fileReader);
    vcfWriter.close();

    return 0;
}
 

/**
 * We want the headers that are used to create the merged header to have the same sequence dictionary
 * that was  returned from the data source and used during validation (which may or may not be the
 * one that was embedded in the input file itself), so get the embedded one from the data source and
 * update it to include the actual sequence dictionary.
 */
private VCFHeader getHeaderWithUpdatedSequenceDictionary(final FeatureDataSource<VariantContext> dataSource) {
    final VCFHeader header = (VCFHeader) dataSource.getHeader();
    if (header.getSequenceDictionary() == null && dataSource.getSequenceDictionary() != null) {
        header.setSequenceDictionary(dataSource.getSequenceDictionary());
    }
    return header;
}
 

@Override
public void onTraversalStart() {
    VCFHeader inputHeader = getHeaderForVariants();
    VCFHeader outputHeader = inputHeader == null ?
            new VCFHeader() :
            new VCFHeader(inputHeader.getMetaDataInInputOrder(), inputHeader.getGenotypeSamples()) ;
    getDefaultToolVCFHeaderLines().forEach(line -> outputHeader.addMetaDataLine(line));
    sourceDictionary = getBestAvailableSequenceDictionary();

    // If -replace is set, do not need to check the sequence dictionary for validity here -- it will still be
    // checked in our normal sequence dictionary validation. Warn and require opt-in via -replace if we're about to
    // clobber a valid sequence dictionary. Check the input file directly via the header rather than using the
    // engine, since it might dig one up from an index.
    if (!replace) {
        SAMSequenceDictionary oldDictionary =
                inputHeader == null ? null : inputHeader.getSequenceDictionary();
        if (oldDictionary != null && !oldDictionary.getSequences().isEmpty())  {
            throw new CommandLineException.BadArgumentValue(
                    String.format(
                            "The input variant file %s already contains a sequence dictionary. " +
                                    "Use %s to force the dictionary to be replaced.",
                            getDrivingVariantsFeatureInput().getName(),
                            REPLACE_ARGUMENT_NAME
                    )
            );
        }

    }

    outputHeader.setSequenceDictionary(sourceDictionary);
    vcfWriter = createVCFWriter(new File(outFile));
    vcfWriter.writeHeader(outputHeader);
}
 

public void writeHeader(final VariantContextWriter vcfWriter, final Set<VCFHeaderLine> defaultToolHeaderLines) {
    final Set<VCFHeaderLine> headerInfo = new HashSet<>();
    headerInfo.add(new VCFHeaderLine("MutectVersion", MUTECT_VERSION));
    headerInfo.add(new VCFHeaderLine(FilterMutectCalls.FILTERING_STATUS_VCF_KEY, "Warning: unfiltered Mutect 2 calls.  Please run " + FilterMutectCalls.class.getSimpleName() + " to remove false positives."));
    headerInfo.addAll(annotationEngine.getVCFAnnotationDescriptions(false));
    headerInfo.addAll(defaultToolHeaderLines);
    STANDARD_MUTECT_INFO_FIELDS.stream().map(GATKVCFHeaderLines::getInfoLine).forEach(headerInfo::add);

    VCFStandardHeaderLines.addStandardFormatLines(headerInfo, true,
            VCFConstants.GENOTYPE_KEY,
            VCFConstants.GENOTYPE_ALLELE_DEPTHS,
            VCFConstants.GENOTYPE_QUALITY_KEY,
            VCFConstants.DEPTH_KEY,
            VCFConstants.GENOTYPE_PL_KEY);
    headerInfo.add(GATKVCFHeaderLines.getFormatLine(GATKVCFConstants.ALLELE_FRACTION_KEY));
    headerInfo.add(GATKVCFHeaderLines.getFormatLine(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_ID_KEY));
    headerInfo.add(GATKVCFHeaderLines.getFormatLine(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_GT_KEY));
    headerInfo.add(VCFStandardHeaderLines.getFormatLine(VCFConstants.PHASE_SET_KEY));

    ReadUtils.getSamplesFromHeader(header).stream().filter(this::isTumorSample)
            .forEach(s -> headerInfo.add(new VCFHeaderLine(TUMOR_SAMPLE_KEY_IN_VCF_HEADER, s)));

    normalSamples.forEach(sample -> headerInfo.add(new VCFHeaderLine(NORMAL_SAMPLE_KEY_IN_VCF_HEADER, sample)));
    if (emitReferenceConfidence()) {
        headerInfo.addAll(referenceConfidenceModel.getVCFHeaderLines());
        headerInfo.add(GATKVCFHeaderLines.getFormatLine(GATKVCFConstants.TUMOR_LOG_10_ODDS_KEY));
    }

    final VCFHeader vcfHeader = new VCFHeader(headerInfo, samplesList.asListOfSamples());
    vcfHeader.setSequenceDictionary(header.getSequenceDictionary());
    vcfWriter.writeHeader(vcfHeader);
}
 

@VisibleForTesting
static VCFHeader getVcfHeader(final SAMSequenceDictionary referenceSequenceDictionary) {
    final Set<VCFHeaderLine> headerLines = new HashSet<>(GATKSVVCFHeaderLines.getSymbAltAlleleLines());
    headerLines.addAll(GATKSVVCFHeaderLines.getInfoLines());
    headerLines.add(VCFStandardHeaderLines.getInfoLine(VCFConstants.END_KEY));
    headerLines.addAll(GATKSVVCFHeaderLines.getFormatLines());
    headerLines.addAll(GATKSVVCFHeaderLines.getFilterLines());
    final VCFHeader header = new VCFHeader(new VCFHeader( headerLines ));
    header.setSequenceDictionary(referenceSequenceDictionary);
    return header;
}
 

private static VCFHeader getHeader() {
    final Set<VCFHeaderLine> headerlines = new LinkedHashSet<>();
    VCFStandardHeaderLines.addStandardFormatLines(headerlines, true,
                                                  VCFConstants.GENOTYPE_KEY,
                                                  VCFConstants.GENOTYPE_QUALITY_KEY,
                                                  VCFConstants.GENOTYPE_PL_KEY, VCFConstants.DEPTH_KEY);
    final SAMSequenceDictionary dict = new SAMSequenceDictionary(
            Collections.singletonList(new SAMSequenceRecord("1", 100)));
    final VCFHeader header = new VCFHeader(headerlines, Collections.singleton(SAMPLE));
    header.setSequenceDictionary(dict);
    return header;
}
 

@Test(dataProvider="UpdateGoodSequenceDictionaryData")
public void testGoodUpdateSequenceDictionary(
        final File inputVariantsFile,
        final File inputSourceFile,
        final File inputReferenceFile,
        final String masterSequenceDictionary,
        final boolean replace,
        final boolean disableSequenceDictionaryValidation) throws FileNotFoundException, URISyntaxException {
    final SAMSequenceDictionary resultingDictionary =
            updateSequenceDictionary(inputVariantsFile, inputSourceFile, inputReferenceFile, masterSequenceDictionary, replace, disableSequenceDictionaryValidation);

    // get the original sequence dictionary from the source for comparison
    SAMSequenceDictionary sourceDictionary =
            SAMSequenceDictionaryExtractor.extractDictionary(
                    inputSourceFile == null ?
                            inputReferenceFile.toPath() : inputSourceFile.toPath()
            );

    // Some sequence dictionary sources will contain optional attributes (i.e., if the source is a .dict file,
    // or if only a reference is presented to the tool using -R, which will in turn cause the framework to retrieve
    // the dictionary from the .dict file accompanying the reference, the dictionary will likely include md5 and
    // UR attributes). However, htsjdk doesn't propagate these attributes to the VCF header properly
    // (https://github.com/samtools/htsjdk/issues/730), and many are stripped out. In order to ensure the
    // roundtrip comparison succeeds, roundtrip it through a VCF header to match what htsjdk will have written out.
    VCFHeader sourceVCFHeader = new VCFHeader();
    sourceVCFHeader.setSequenceDictionary(sourceDictionary);
    sourceDictionary = sourceVCFHeader.getSequenceDictionary();
    Assert.assertEquals(sourceDictionary, resultingDictionary);
}
 

private static File createInputVCF(final String sampleName) {
    final String contig = "chr20";
    final SAMSequenceDictionary dict = new SAMSequenceDictionary(
            Collections.singletonList(new SAMSequenceRecord(contig, 64444167)));

    final VCFFormatHeaderLine formatField = new VCFFormatHeaderLine(SAMPLE_NAME_KEY, 1, VCFHeaderLineType.String,
                                                                    "the name of the sample this genotype came from");
    final Set<VCFHeaderLine> headerLines = new HashSet<>();
    headerLines.add(formatField);
    headerLines.add(new VCFFormatHeaderLine(ANOTHER_ATTRIBUTE_KEY, 1, VCFHeaderLineType.Integer, "Another value"));
    headerLines.add(VCFStandardHeaderLines.getFormatLine("GT"));

    final File out = createTempFile(sampleName +"_", ".vcf");
    try (final VariantContextWriter writer = GATKVariantContextUtils.createVCFWriter(out.toPath(), dict, false,
                                                                                     Options.INDEX_ON_THE_FLY)) {
        final VCFHeader vcfHeader = new VCFHeader(headerLines, Collections.singleton(sampleName));
        vcfHeader.setSequenceDictionary(dict);
        writer.writeHeader(vcfHeader);
        final Allele Aref = Allele.create("A", true);
        final Allele C = Allele.create("C");
        final List<Allele> alleles = Arrays.asList(Aref, C);
        final VariantContext variant = new VariantContextBuilder("invented", contig, INTERVAL.get(0).getStart(), INTERVAL.get(0).getStart(), alleles)
                .genotypes(new GenotypeBuilder(sampleName, alleles).attribute(SAMPLE_NAME_KEY, sampleName)
                                   .attribute(ANOTHER_ATTRIBUTE_KEY, 10).make())
                .make();
        writer.add(variant);
        return out;
    }
}
 

@Override
protected int doWork() {
    final ProgressLogger progress = new ProgressLogger(log, 10000);
    final List<String> sampleList = new ArrayList<String>();
    INPUT = IOUtil.unrollFiles(INPUT, IOUtil.VCF_EXTENSIONS);
    final Collection<CloseableIterator<VariantContext>> iteratorCollection = new ArrayList<CloseableIterator<VariantContext>>(INPUT.size());
    final Collection<VCFHeader> headers = new HashSet<VCFHeader>(INPUT.size());
    VariantContextComparator variantContextComparator = null;
    SAMSequenceDictionary sequenceDictionary = null;

    if (SEQUENCE_DICTIONARY != null) {
        sequenceDictionary = SamReaderFactory.makeDefault().referenceSequence(REFERENCE_SEQUENCE).open(SEQUENCE_DICTIONARY).getFileHeader().getSequenceDictionary();
    }

    for (final File file : INPUT) {
        IOUtil.assertFileIsReadable(file);
        final VCFFileReader fileReader = new VCFFileReader(file, false);
        final VCFHeader fileHeader = fileReader.getFileHeader();
        if (fileHeader.getContigLines().isEmpty()) {
            if (sequenceDictionary == null) {
                throw new IllegalArgumentException(SEQ_DICT_REQUIRED);
            } else {
                fileHeader.setSequenceDictionary(sequenceDictionary);
            }
        }

        if (variantContextComparator == null) {
            variantContextComparator = fileHeader.getVCFRecordComparator();
        } else {
            if (!variantContextComparator.isCompatible(fileHeader.getContigLines())) {
                throw new IllegalArgumentException(
                        "The contig entries in input file " + file.getAbsolutePath() + " are not compatible with the others.");
            }
        }

        if (sequenceDictionary == null) sequenceDictionary = fileHeader.getSequenceDictionary();

        if (sampleList.isEmpty()) {
            sampleList.addAll(fileHeader.getSampleNamesInOrder());
        } else {
            if (!sampleList.equals(fileHeader.getSampleNamesInOrder())) {
                throw new IllegalArgumentException("Input file " + file.getAbsolutePath() + " has sample entries that don't match the other files.");
            }
        }
        
        // add comments in the first header
        if (headers.isEmpty()) {
            COMMENT.stream().forEach(C -> fileHeader.addMetaDataLine(new VCFHeaderLine("MergeVcfs.comment", C)));
        }

        headers.add(fileHeader);
        iteratorCollection.add(fileReader.iterator());
    }

    if (CREATE_INDEX && sequenceDictionary == null) {
        throw new PicardException(String.format("Index creation failed. %s", SEQ_DICT_REQUIRED));
    }

    final VariantContextWriterBuilder builder = new VariantContextWriterBuilder()
            .setOutputFile(OUTPUT)
            .setReferenceDictionary(sequenceDictionary);

    if (CREATE_INDEX) {
        builder.setOption(Options.INDEX_ON_THE_FLY);
    } else {
        builder.unsetOption(Options.INDEX_ON_THE_FLY);
    }
    final VariantContextWriter writer = builder.build();

    writer.writeHeader(new VCFHeader(VCFUtils.smartMergeHeaders(headers, false), sampleList));

    final MergingIterator<VariantContext> mergingIterator = new MergingIterator<VariantContext>(variantContextComparator, iteratorCollection);
    while (mergingIterator.hasNext()) {
        final VariantContext context = mergingIterator.next();
        writer.add(context);
        progress.record(context.getContig(), context.getStart());
    }

    CloserUtil.close(mergingIterator);
    writer.close();
    return 0;
}
 

/**
 * Create a VCF header.
 *
 * @param sequenceDictionary sequence dictionary for the reads
 * @return a VCF header
 */
public VCFHeader makeVCFHeader( final SAMSequenceDictionary sequenceDictionary, final Set<VCFHeaderLine>  defaultToolHeaderLines ) {
    final Set<VCFHeaderLine> headerInfo = new HashSet<>();
    headerInfo.addAll(defaultToolHeaderLines);

    headerInfo.addAll(genotypingEngine.getAppropriateVCFInfoHeaders());
    // all annotation fields from VariantAnnotatorEngine
    headerInfo.addAll(annotationEngine.getVCFAnnotationDescriptions(emitReferenceConfidence()));
    // all callers need to add these standard annotation header lines
    headerInfo.add(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.MLE_ALLELE_COUNT_KEY));
    headerInfo.add(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.MLE_ALLELE_FREQUENCY_KEY));
    // all callers need to add these standard FORMAT field header lines
    VCFStandardHeaderLines.addStandardFormatLines(headerInfo, true,
            VCFConstants.GENOTYPE_KEY,
            VCFConstants.GENOTYPE_QUALITY_KEY,
            VCFConstants.DEPTH_KEY,
            VCFConstants.GENOTYPE_PL_KEY);

    if (hcArgs.standardArgs.genotypeArgs.supportVariants != null) {
        headerInfo.add(VCFStandardHeaderLines.getInfoLine(VCFConstants.ALLELE_COUNT_KEY));
        headerInfo.add(VCFStandardHeaderLines.getInfoLine(VCFConstants.ALLELE_FREQUENCY_KEY));
        headerInfo.add(VCFStandardHeaderLines.getInfoLine(VCFConstants.ALLELE_NUMBER_KEY));
        headerInfo.add(GATKVCFHeaderLines.getFormatLine(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY));
        headerInfo.add(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.GENOTYPE_PRIOR_KEY));
    }

    if ( ! hcArgs.doNotRunPhysicalPhasing ) {
        headerInfo.add(GATKVCFHeaderLines.getFormatLine(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_ID_KEY));
        headerInfo.add(GATKVCFHeaderLines.getFormatLine(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_GT_KEY));
        headerInfo.add(VCFStandardHeaderLines.getFormatLine(VCFConstants.PHASE_SET_KEY));
    }

    // FILTER fields are added unconditionally as it's not always 100% certain the circumstances
    // where the filters are used.  For example, in emitting all sites the lowQual field is used
    headerInfo.add(GATKVCFHeaderLines.getFilterLine(GATKVCFConstants.LOW_QUAL_FILTER_NAME));

    if ( emitReferenceConfidence() ) {
        headerInfo.addAll(referenceConfidenceModel.getVCFHeaderLines());
    }

    final VCFHeader vcfHeader = new VCFHeader(headerInfo, sampleSet);
    vcfHeader.setSequenceDictionary(sequenceDictionary);
    return vcfHeader;
}