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