Search in sources :

Example 56 with VariantContext

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

the class GermlineProbabilityCalculator method calculateAnnotations.

public static Map<String, Object> calculateAnnotations(List<VariantContext> germlineResourceVariants, final List<Allele> altAlleles, final double[] tumorLog10Odds, final Optional<double[]> normalLog10Odds, final double afOfAllelesNotInGermlineResource, final double log10PriorProbOfSomaticEvent) {
    final double[] normalLog10OddsOrFlat = normalLog10Odds.orElseGet(() -> MathUtils.applyToArray(tumorLog10Odds, x -> 0));
    final Optional<VariantContext> germlineVC = germlineResourceVariants.isEmpty() ? Optional.empty() : // assume only one VC per site
    Optional.of(germlineResourceVariants.get(0));
    final double[] populationAlleleFrequencies = getGermlineAltAlleleFrequencies(altAlleles, germlineVC, afOfAllelesNotInGermlineResource);
    // note the minus sign required because Mutect has the convention that this is log odds of allele *NOT* being in the normal
    final double[] germlineLog10Posteriors = new IndexRange(0, altAlleles.size()).mapToDouble(n -> log10PosteriorProbabilityOfGermlineVariant(-normalLog10OddsOrFlat[n], tumorLog10Odds[n], populationAlleleFrequencies[n], log10PriorProbOfSomaticEvent));
    return ImmutableMap.of(POPULATION_AF_VCF_ATTRIBUTE, populationAlleleFrequencies, GERMLINE_POSTERIORS_VCF_ATTRIBUTE, germlineLog10Posteriors);
}
Also used : IntStream(java.util.stream.IntStream) Allele(htsjdk.variant.variantcontext.Allele) java.util(java.util) MathUtils(org.broadinstitute.hellbender.utils.MathUtils) Doubles(com.google.cloud.dataflow.sdk.repackaged.com.google.common.primitives.Doubles) ImmutableMap(com.google.common.collect.ImmutableMap) GATKProtectedVariantContextUtils(org.broadinstitute.hellbender.utils.GATKProtectedVariantContextUtils) VariantContext(htsjdk.variant.variantcontext.VariantContext) VisibleForTesting(com.google.common.annotations.VisibleForTesting) IndexRange(org.broadinstitute.hellbender.utils.IndexRange) VCFConstants(htsjdk.variant.vcf.VCFConstants) IndexRange(org.broadinstitute.hellbender.utils.IndexRange) VariantContext(htsjdk.variant.variantcontext.VariantContext)

Example 57 with VariantContext

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

the class CalculateMixingFractions method apply.

@Override
public void apply(final VariantContext vc, final ReadsContext readsContext, final ReferenceContext refContext, final FeatureContext fc) {
    if (!isBiallelicSingletonHetSnp(vc)) {
        return;
    }
    final Optional<String> variantSample = StreamSupport.stream(vc.getGenotypes().spliterator(), false).filter(genotype -> genotype.isHet()).map(genotype -> genotype.getSampleName()).findFirst();
    if (!variantSample.isPresent()) {
        return;
    }
    final List<GATKRead> reads = new ArrayList<>();
    final List<Integer> offsets = new ArrayList<>();
    for (final GATKRead read : readsContext) {
        if (read.failsVendorQualityCheck()) {
            continue;
        }
        final AlignmentStateMachine asm = new AlignmentStateMachine(read);
        while (asm.stepForwardOnGenome() != null && asm.getGenomePosition() < vc.getStart()) {
        }
        if (asm.getGenomePosition() == vc.getStart()) {
            reads.add(read);
            offsets.add(asm.getReadOffset());
        }
    }
    final ReadPileup pileup = new ReadPileup(vc, reads, offsets);
    final byte altBase = vc.getAlternateAllele(0).getBases()[0];
    final long altCount = StreamSupport.stream(pileup.spliterator(), false).filter(pe -> pe.getBase() == altBase).count();
    final long totalCount = pileup.size();
    sampleCounts.get(variantSample.get()).addCounts(altCount, totalCount);
}
Also used : DocumentedFeature(org.broadinstitute.barclay.help.DocumentedFeature) htsjdk.variant.vcf(htsjdk.variant.vcf) CommandLineProgramProperties(org.broadinstitute.barclay.argparser.CommandLineProgramProperties) java.util(java.util) AlignmentStateMachine(org.broadinstitute.hellbender.utils.locusiterator.AlignmentStateMachine) ReadPileup(org.broadinstitute.hellbender.utils.pileup.ReadPileup) Argument(org.broadinstitute.barclay.argparser.Argument) StandardArgumentDefinitions(org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions) GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) Collectors(java.util.stream.Collectors) VariantProgramGroup(org.broadinstitute.hellbender.cmdline.programgroups.VariantProgramGroup) File(java.io.File) org.broadinstitute.hellbender.engine(org.broadinstitute.hellbender.engine) GATKProtectedVariantContextUtils(org.broadinstitute.hellbender.utils.GATKProtectedVariantContextUtils) VariantContext(htsjdk.variant.variantcontext.VariantContext) StreamSupport(java.util.stream.StreamSupport) MutableLong(org.apache.commons.lang.mutable.MutableLong) GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) ReadPileup(org.broadinstitute.hellbender.utils.pileup.ReadPileup) AlignmentStateMachine(org.broadinstitute.hellbender.utils.locusiterator.AlignmentStateMachine)

Example 58 with VariantContext

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

the class ReadThreadingAssembler method composeGivenHaplotypes.

/**
     * Create the list of artificial GGA-mode haplotypes by injecting each of the provided alternate alleles into the reference haplotype
     *
     * @param refHaplotype the reference haplotype
     * @param givenHaplotypes the list of alternate alleles in VariantContexts
     * @param activeRegionWindow the window containing the reference haplotype
     *
     * @return a non-null list of haplotypes
     */
private static List<Haplotype> composeGivenHaplotypes(final Haplotype refHaplotype, final List<VariantContext> givenHaplotypes, final SimpleInterval activeRegionWindow) {
    Utils.nonNull(refHaplotype, "the reference haplotype cannot be null");
    Utils.nonNull(givenHaplotypes, "given haplotypes cannot be null");
    Utils.nonNull(activeRegionWindow, "active region window cannot be null");
    Utils.validateArg(activeRegionWindow.size() == refHaplotype.length(), "inconsistent reference haplotype and active region window");
    final Set<Haplotype> returnHaplotypes = new LinkedHashSet<>();
    final int activeRegionStart = refHaplotype.getAlignmentStartHapwrtRef();
    for (final VariantContext compVC : givenHaplotypes) {
        Utils.validateArg(GATKVariantContextUtils.overlapsRegion(compVC, activeRegionWindow), " some variant provided does not overlap with active region window");
        for (final Allele compAltAllele : compVC.getAlternateAlleles()) {
            final Haplotype insertedRefHaplotype = refHaplotype.insertAllele(compVC.getReference(), compAltAllele, activeRegionStart + compVC.getStart() - activeRegionWindow.getStart(), compVC.getStart());
            if (insertedRefHaplotype != null) {
                // can be null if the requested allele can't be inserted into the haplotype
                returnHaplotypes.add(insertedRefHaplotype);
            }
        }
    }
    return new ArrayList<>(returnHaplotypes);
}
Also used : Allele(htsjdk.variant.variantcontext.Allele) VariantContext(htsjdk.variant.variantcontext.VariantContext) Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype)

Example 59 with VariantContext

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

the class ASEReadCounter method apply.

@Override
public void apply(AlignmentContext alignmentContext, ReferenceContext referenceContext, FeatureContext featureContext) {
    final String contig = alignmentContext.getContig();
    final long position = alignmentContext.getPosition();
    final char refAllele = (char) referenceContext.getBase();
    final List<VariantContext> VCs = featureContext.getValues(variants);
    if (VCs != null && VCs.size() > 1) {
        throw new UserException("More then one variant context at position: " + contig + ":" + position);
    }
    if (VCs == null || VCs.isEmpty()) {
        return;
    }
    final VariantContext vc = VCs.get(0);
    if (!vc.isBiallelic()) {
        logger.warn("Ignoring site: cannot run ASE on non-biallelic sites: " + vc.toString());
        return;
    }
    if (vc.getHetCount() < 1) {
        logger.warn("Ignoring site: variant is not het at postion: " + contig + ":" + position);
        return;
    }
    if (vc.getNAlleles() == 1 || vc.getAlternateAllele(0).getBases().length == 0) {
        throw new UserException("The file of variant sites must contain heterozygous sites and cannot be a GVCF file containing <NON_REF> alleles.");
    }
    final char altAllele = (char) vc.getAlternateAllele(0).getBases()[0];
    final String siteID = vc.getID();
    final ReadPileup pileup = filterPileup(alignmentContext.getBasePileup(), countType);
    // count up the depths of all and QC+ bases
    final String line = calculateLineForSite(pileup, siteID, refAllele, altAllele);
    if (line != null) {
        outputStream.println(line);
    }
}
Also used : ReadPileup(org.broadinstitute.hellbender.utils.pileup.ReadPileup) VariantContext(htsjdk.variant.variantcontext.VariantContext) UserException(org.broadinstitute.hellbender.exceptions.UserException)

Example 60 with VariantContext

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

the class FeatureDataSourceUnitTest method testGetSequenceDictionary.

@Test
public void testGetSequenceDictionary() {
    try (FeatureDataSource<VariantContext> featureSource = new FeatureDataSource<>(QUERY_TEST_VCF, "CustomName")) {
        final SAMSequenceDictionary dict = featureSource.getSequenceDictionary();
        Assert.assertEquals(dict.size(), 4);
        Assert.assertEquals(dict.getSequences().stream().map(s -> s.getSequenceName()).collect(Collectors.toList()), Arrays.asList("1", "2", "3", "4"));
    }
}
Also used : VariantContext(htsjdk.variant.variantcontext.VariantContext) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Aggregations

VariantContext (htsjdk.variant.variantcontext.VariantContext)237 Test (org.testng.annotations.Test)119 File (java.io.File)98 BaseTest (org.broadinstitute.hellbender.utils.test.BaseTest)68 Allele (htsjdk.variant.variantcontext.Allele)55 Genotype (htsjdk.variant.variantcontext.Genotype)49 VariantContextBuilder (htsjdk.variant.variantcontext.VariantContextBuilder)43 Collectors (java.util.stream.Collectors)43 SimpleInterval (org.broadinstitute.hellbender.utils.SimpleInterval)43 CommandLineProgramTest (org.broadinstitute.hellbender.CommandLineProgramTest)41 FeatureDataSource (org.broadinstitute.hellbender.engine.FeatureDataSource)36 StandardArgumentDefinitions (org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions)35 java.util (java.util)33 IOException (java.io.IOException)25 Assert (org.testng.Assert)25 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)24 VCFFileReader (htsjdk.variant.vcf.VCFFileReader)22 UserException (org.broadinstitute.hellbender.exceptions.UserException)22 Target (org.broadinstitute.hellbender.tools.exome.Target)22 DataProvider (org.testng.annotations.DataProvider)22