Search in sources :

Example 46 with VariantContextBuilder

use of htsjdk.variant.variantcontext.VariantContextBuilder in project gatk by broadinstitute.

the class HomRefBlockUnitTest method testIsContiguous.

@Test(dataProvider = "ContiguousData")
public void testIsContiguous(final String contig, final int pos, final boolean expected) {
    final VariantContext vc = getVariantContext();
    final HomRefBlock band = getHomRefBlock(vc);
    final VariantContext testVC = new VariantContextBuilder(vc).chr(contig).start(pos).stop(pos).make();
    Assert.assertEquals(band.isContiguous(testVC), expected);
}
Also used : VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) VariantContext(htsjdk.variant.variantcontext.VariantContext) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 47 with VariantContextBuilder

use of htsjdk.variant.variantcontext.VariantContextBuilder in project gatk-protected by broadinstitute.

the class ReadThreadingAssemblerUnitTest method testAssembleRefAndInsertion.

@Test(dataProvider = "AssembleIntervalsWithVariantData")
public void testAssembleRefAndInsertion(final ReadThreadingAssembler assembler, final SimpleInterval loc, final int nReadsToUse, final int variantSite) {
    final byte[] refBases = seq.getSubsequenceAt(loc.getContig(), loc.getStart(), loc.getEnd()).getBases();
    for (int insertionLength = 1; insertionLength < 10; insertionLength++) {
        final Allele refBase = Allele.create(refBases[variantSite], false);
        final Allele altBase = Allele.create(new String(refBases).substring(variantSite, variantSite + insertionLength + 1), true);
        final VariantContextBuilder vcb = new VariantContextBuilder("x", loc.getContig(), variantSite, variantSite + insertionLength, Arrays.asList(refBase, altBase));
        testAssemblyWithVariant(assembler, refBases, loc, nReadsToUse, vcb.make());
    }
}
Also used : Allele(htsjdk.variant.variantcontext.Allele) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 48 with VariantContextBuilder

use of htsjdk.variant.variantcontext.VariantContextBuilder in project gatk-protected by broadinstitute.

the class ReadThreadingAssemblerUnitTest method testAssembleRefAndDeletion.

@Test(dataProvider = "AssembleIntervalsWithVariantData")
public void testAssembleRefAndDeletion(final ReadThreadingAssembler assembler, final SimpleInterval loc, final int nReadsToUse, final int variantSite) {
    final byte[] refBases = seq.getSubsequenceAt(loc.getContig(), loc.getStart(), loc.getEnd()).getBases();
    for (int deletionLength = 1; deletionLength < 10; deletionLength++) {
        final Allele refBase = Allele.create(new String(refBases).substring(variantSite, variantSite + deletionLength + 1), true);
        final Allele altBase = Allele.create(refBase.getBases()[0], false);
        final VariantContextBuilder vcb = new VariantContextBuilder("x", loc.getContig(), variantSite, variantSite + deletionLength, Arrays.asList(refBase, altBase));
        testAssemblyWithVariant(assembler, refBases, loc, nReadsToUse, vcb.make());
    }
}
Also used : Allele(htsjdk.variant.variantcontext.Allele) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 49 with VariantContextBuilder

use of htsjdk.variant.variantcontext.VariantContextBuilder in project gatk-protected by broadinstitute.

the class ReadThreadingAssemblerUnitTest method testAssembleRefAndSNP.

@Test(dataProvider = "AssembleIntervalsWithVariantData")
public void testAssembleRefAndSNP(final ReadThreadingAssembler assembler, final SimpleInterval loc, final int nReadsToUse, final int variantSite) {
    final byte[] refBases = seq.getSubsequenceAt(loc.getContig(), loc.getStart(), loc.getEnd()).getBases();
    final Allele refBase = Allele.create(refBases[variantSite], true);
    final Allele altBase = Allele.create((byte) (refBase.getBases()[0] == 'A' ? 'C' : 'A'), false);
    final VariantContextBuilder vcb = new VariantContextBuilder("x", loc.getContig(), variantSite, variantSite, Arrays.asList(refBase, altBase));
    testAssemblyWithVariant(assembler, refBases, loc, nReadsToUse, vcb.make());
}
Also used : Allele(htsjdk.variant.variantcontext.Allele) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 50 with VariantContextBuilder

use of htsjdk.variant.variantcontext.VariantContextBuilder in project gatk by broadinstitute.

the class LiftOverVcf method doWork.

@Override
protected Object doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsReadable(REFERENCE_SEQUENCE);
    IOUtil.assertFileIsReadable(CHAIN);
    IOUtil.assertFileIsWritable(OUTPUT);
    IOUtil.assertFileIsWritable(REJECT);
    ////////////////////////////////////////////////////////////////////////
    // Setup the inputs
    ////////////////////////////////////////////////////////////////////////
    final LiftOver liftOver = new LiftOver(CHAIN);
    final VCFFileReader in = new VCFFileReader(INPUT, false);
    logger.info("Loading up the target reference genome.");
    final ReferenceSequenceFileWalker walker = new ReferenceSequenceFileWalker(REFERENCE_SEQUENCE);
    final Map<String, byte[]> refSeqs = new HashMap<>();
    for (final SAMSequenceRecord rec : walker.getSequenceDictionary().getSequences()) {
        refSeqs.put(rec.getSequenceName(), walker.get(rec.getSequenceIndex()).getBases());
    }
    CloserUtil.close(walker);
    ////////////////////////////////////////////////////////////////////////
    // Setup the outputs
    ////////////////////////////////////////////////////////////////////////
    final VCFHeader inHeader = in.getFileHeader();
    final VCFHeader outHeader = new VCFHeader(inHeader);
    outHeader.setSequenceDictionary(walker.getSequenceDictionary());
    final VariantContextWriter out = new VariantContextWriterBuilder().setOption(Options.INDEX_ON_THE_FLY).setOutputFile(OUTPUT).setReferenceDictionary(walker.getSequenceDictionary()).build();
    out.writeHeader(outHeader);
    final VariantContextWriter rejects = new VariantContextWriterBuilder().setOutputFile(REJECT).unsetOption(Options.INDEX_ON_THE_FLY).build();
    final VCFHeader rejectHeader = new VCFHeader(in.getFileHeader());
    for (final VCFFilterHeaderLine line : FILTERS) rejectHeader.addMetaDataLine(line);
    rejects.writeHeader(rejectHeader);
    ////////////////////////////////////////////////////////////////////////
    // Read the input VCF, lift the records over and write to the sorting
    // collection.
    ////////////////////////////////////////////////////////////////////////
    long failedLiftover = 0, failedAlleleCheck = 0, total = 0;
    logger.info("Lifting variants over and sorting.");
    final SortingCollection<VariantContext> sorter = SortingCollection.newInstance(VariantContext.class, new VCFRecordCodec(outHeader), outHeader.getVCFRecordComparator(), MAX_RECORDS_IN_RAM, TMP_DIR);
    ProgressLogger progress = new ProgressLogger(logger, 1000000, "read");
    for (final VariantContext ctx : in) {
        ++total;
        final Interval source = new Interval(ctx.getContig(), ctx.getStart(), ctx.getEnd(), false, ctx.getContig() + ":" + ctx.getStart() + "-" + ctx.getEnd());
        final Interval target = liftOver.liftOver(source, 1.0);
        if (target == null) {
            rejects.add(new VariantContextBuilder(ctx).filter(FILTER_CANNOT_LIFTOVER).make());
            failedLiftover++;
        } else {
            // Fix the alleles if we went from positive to negative strand
            final List<Allele> alleles = new ArrayList<>();
            for (final Allele oldAllele : ctx.getAlleles()) {
                if (target.isPositiveStrand() || oldAllele.isSymbolic()) {
                    alleles.add(oldAllele);
                } else {
                    alleles.add(Allele.create(SequenceUtil.reverseComplement(oldAllele.getBaseString()), oldAllele.isReference()));
                }
            }
            // Build the new variant context
            final VariantContextBuilder builder = new VariantContextBuilder(ctx.getSource(), target.getContig(), target.getStart(), target.getEnd(), alleles);
            builder.id(ctx.getID());
            builder.attributes(ctx.getAttributes());
            builder.genotypes(ctx.getGenotypes());
            builder.filters(ctx.getFilters());
            builder.log10PError(ctx.getLog10PError());
            // Check that the reference allele still agrees with the reference sequence
            boolean mismatchesReference = false;
            for (final Allele allele : builder.getAlleles()) {
                if (allele.isReference()) {
                    final byte[] ref = refSeqs.get(target.getContig());
                    final String refString = StringUtil.bytesToString(ref, target.getStart() - 1, target.length());
                    if (!refString.equalsIgnoreCase(allele.getBaseString())) {
                        mismatchesReference = true;
                    }
                    break;
                }
            }
            if (mismatchesReference) {
                rejects.add(new VariantContextBuilder(ctx).filter(FILTER_MISMATCHING_REF_ALLELE).make());
                failedAlleleCheck++;
            } else {
                sorter.add(builder.make());
            }
        }
        progress.record(ctx.getContig(), ctx.getStart());
    }
    final NumberFormat pfmt = new DecimalFormat("0.0000%");
    final String pct = pfmt.format((failedLiftover + failedAlleleCheck) / (double) total);
    logger.info("Processed ", total, " variants.");
    logger.info(Long.toString(failedLiftover), " variants failed to liftover.");
    logger.info(Long.toString(failedAlleleCheck), " variants lifted over but had mismatching reference alleles after lift over.");
    logger.info(pct, " of variants were not successfully lifted over and written to the output.");
    rejects.close();
    in.close();
    ////////////////////////////////////////////////////////////////////////
    // Write the sorted outputs to the final output file
    ////////////////////////////////////////////////////////////////////////
    sorter.doneAdding();
    progress = new ProgressLogger(logger, 1000000, "written");
    logger.info("Writing out sorted records to final VCF.");
    for (final VariantContext ctx : sorter) {
        out.add(ctx);
        progress.record(ctx.getContig(), ctx.getStart());
    }
    out.close();
    sorter.cleanup();
    return null;
}
Also used : LiftOver(htsjdk.samtools.liftover.LiftOver) HashMap(java.util.HashMap) DecimalFormat(java.text.DecimalFormat) VCFFileReader(htsjdk.variant.vcf.VCFFileReader) ArrayList(java.util.ArrayList) VariantContext(htsjdk.variant.variantcontext.VariantContext) SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord) ProgressLogger(org.broadinstitute.hellbender.utils.runtime.ProgressLogger) VCFRecordCodec(htsjdk.variant.vcf.VCFRecordCodec) VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter) VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine) ReferenceSequenceFileWalker(htsjdk.samtools.reference.ReferenceSequenceFileWalker) VCFHeader(htsjdk.variant.vcf.VCFHeader) Allele(htsjdk.variant.variantcontext.Allele) VariantContextWriterBuilder(htsjdk.variant.variantcontext.writer.VariantContextWriterBuilder) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) NumberFormat(java.text.NumberFormat)

Aggregations

VariantContextBuilder (htsjdk.variant.variantcontext.VariantContextBuilder)57 VariantContext (htsjdk.variant.variantcontext.VariantContext)40 Test (org.testng.annotations.Test)26 Allele (htsjdk.variant.variantcontext.Allele)25 BaseTest (org.broadinstitute.hellbender.utils.test.BaseTest)21 Genotype (htsjdk.variant.variantcontext.Genotype)12 ReferenceContext (org.broadinstitute.hellbender.engine.ReferenceContext)8 SimpleInterval (org.broadinstitute.hellbender.utils.SimpleInterval)7 GenotypeBuilder (htsjdk.variant.variantcontext.GenotypeBuilder)6 GenotypesContext (htsjdk.variant.variantcontext.GenotypesContext)6 File (java.io.File)6 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)5 VCFHeader (htsjdk.variant.vcf.VCFHeader)5 CommandLineProgramTest (org.broadinstitute.hellbender.CommandLineProgramTest)5 BeforeTest (org.testng.annotations.BeforeTest)5 SAMSequenceRecord (htsjdk.samtools.SAMSequenceRecord)4 VariantContextWriter (htsjdk.variant.variantcontext.writer.VariantContextWriter)4 VCFHeaderLine (htsjdk.variant.vcf.VCFHeaderLine)3 ArrayList (java.util.ArrayList)3 Collectors (java.util.stream.Collectors)3