Search in sources :

Example 11 with IndexRange

use of org.broadinstitute.hellbender.utils.IndexRange 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 12 with IndexRange

use of org.broadinstitute.hellbender.utils.IndexRange in project gatk by broadinstitute.

the class SegmentUtils method trimInterval.

/**
     * Given an interval and collections of targets and SNPs, returns a trimmed interval produced by removing the empty
     * portions at the start and the end of the original interval that do not overlap the targets and SNPs that overlap
     * with the original interval.  If this procedure would remove the entire interval, the original interval is
     * returned instead.  Note that this method will not expand an interval to the start of the first overlapping target
     * and the end of the last overlapping target; it will only shrink the interval or leave it alone.  This is to
     * avoid overlapping segments (which would occur if a SNP breakpoint fell within a target and the interval
     * were expanded, for example).
     */
public static SimpleInterval trimInterval(final Locatable interval, final TargetCollection<? extends Locatable> targets, final TargetCollection<? extends Locatable> snps) {
    Utils.nonNull(interval, "The interval cannot be null.");
    Utils.nonNull(targets, "The collection of targets cannot be null.");
    Utils.nonNull(snps, "The collection of SNPs cannot be null.");
    final IndexRange targetRange = targets.indexRange(interval);
    final IndexRange snpRange = snps.indexRange(interval);
    final int numTargetsInInterval = targetRange.size();
    final int numSNPsInInterval = snpRange.size();
    int start = interval.getStart();
    int end = interval.getEnd();
    if (numTargetsInInterval == 0 && numSNPsInInterval > 0) {
        //if there are no targets overlapping interval, use SNPs to determine trimmed interval
        start = snps.target(snpRange.from).getStart();
        end = snps.target(snpRange.to - 1).getEnd();
    } else if (numTargetsInInterval > 0) {
        //if interval start does not fall within first target, use start of first target as start of trimmed interval
        start = Math.max(start, targets.target(targetRange.from).getStart());
        //if interval end does not fall within last target, use end of last target as end of trimmed interval
        end = Math.min(end, targets.target(targetRange.to - 1).getEnd());
        if (numSNPsInInterval > 0) {
            //if there are also SNPs within interval, check to see if they give a larger trimmed interval
            start = Math.min(start, snps.target(snpRange.from).getStart());
            end = Math.max(end, snps.target(snpRange.to - 1).getEnd());
        }
    }
    if (start < interval.getStart() || end > interval.getEnd() || end < start) {
        throw new GATKException.ShouldNeverReachHereException("Something went wrong in trimming interval.");
    }
    return new SimpleInterval(interval.getContig(), start, end);
}
Also used : IndexRange(org.broadinstitute.hellbender.utils.IndexRange) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval)

Example 13 with IndexRange

use of org.broadinstitute.hellbender.utils.IndexRange in project gatk by broadinstitute.

the class XHMMSegmentCallerIntegrationTest method assertOutputHasConsistentNumberOfTargets.

private void assertOutputHasConsistentNumberOfTargets(final List<HiddenStateSegmentRecord<CopyNumberTriState, Target>> outputRecords, final TargetCollection<Target> targets) {
    for (final HiddenStateSegmentRecord<CopyNumberTriState, Target> nextRecord : outputRecords) {
        final IndexRange indexRange = targets.indexRange(nextRecord.getSegment());
        Assert.assertEquals(indexRange.to - indexRange.from, nextRecord.getSegment().getTargetCount());
    }
}
Also used : IndexRange(org.broadinstitute.hellbender.utils.IndexRange) Target(org.broadinstitute.hellbender.tools.exome.Target) CopyNumberTriState(org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriState)

Example 14 with IndexRange

use of org.broadinstitute.hellbender.utils.IndexRange in project gatk by broadinstitute.

the class XHMMSegmentCallerIntegrationTest method assertSampleSegmentsCoverAllTargets.

private void assertSampleSegmentsCoverAllTargets(final List<HiddenStateSegmentRecord<CopyNumberTriState, Target>> sampleRecords, final TargetCollection<Target> targets) {
    int next = 0;
    for (final HiddenStateSegmentRecord<CopyNumberTriState, Target> record : sampleRecords) {
        final IndexRange range = targets.indexRange(record.getSegment());
        Assert.assertEquals(range.from, next);
        next = range.to;
    }
}
Also used : IndexRange(org.broadinstitute.hellbender.utils.IndexRange) Target(org.broadinstitute.hellbender.tools.exome.Target) CopyNumberTriState(org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriState)

Example 15 with IndexRange

use of org.broadinstitute.hellbender.utils.IndexRange in project gatk by broadinstitute.

the class XHMMSegmentCallerIntegrationTest method assertOutputIsInOrder.

private void assertOutputIsInOrder(final List<HiddenStateSegmentRecord<CopyNumberTriState, Target>> outputRecords, final TargetCollection<Target> targets) {
    for (int i = 1; i < outputRecords.size(); i++) {
        final HiddenStateSegmentRecord<CopyNumberTriState, Target> nextRecord = outputRecords.get(i);
        final HiddenStateSegmentRecord<CopyNumberTriState, Target> previousRecord = outputRecords.get(i - 1);
        final IndexRange nextRange = targets.indexRange(nextRecord.getSegment());
        final IndexRange previousRange = targets.indexRange(previousRecord.getSegment());
        Assert.assertTrue(nextRange.from >= previousRange.from);
    }
}
Also used : IndexRange(org.broadinstitute.hellbender.utils.IndexRange) Target(org.broadinstitute.hellbender.tools.exome.Target) CopyNumberTriState(org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriState)

Aggregations

IndexRange (org.broadinstitute.hellbender.utils.IndexRange)32 Test (org.testng.annotations.Test)11 BaseTest (org.broadinstitute.hellbender.utils.test.BaseTest)9 IntStream (java.util.stream.IntStream)7 Target (org.broadinstitute.hellbender.tools.exome.Target)7 CopyNumberTriState (org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriState)7 MathUtils (org.broadinstitute.hellbender.utils.MathUtils)7 Allele (htsjdk.variant.variantcontext.Allele)5 VariantContext (htsjdk.variant.variantcontext.VariantContext)5 Collectors (java.util.stream.Collectors)5 MathArrays (org.apache.commons.math3.util.MathArrays)5 SimpleInterval (org.broadinstitute.hellbender.utils.SimpleInterval)5 Utils (org.broadinstitute.hellbender.utils.Utils)5 VisibleForTesting (com.google.common.annotations.VisibleForTesting)4 File (java.io.File)4 java.util (java.util)4 Genotype (htsjdk.variant.variantcontext.Genotype)3 Arrays (java.util.Arrays)3 List (java.util.List)3 Map (java.util.Map)3