Examples with AllelicCount org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCount used on opensource projects

Example 31 with AllelicCount

use of org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCount in project gatk by broadinstitute.

the class AlleleFractionHMMUnitTest method obviousCallsTest.

// if we ignore ref bias, choose  well-separated minor fractions e.g. 0.1 and 0.5, and give really obvious
// allele counts like 10/90 (obviously 0.1) and 50/50 (obviously 0.5), the Viterbi algorithm should make the right calls
// this essentially tests whether we correctly defined the likelihood in terms of methods in the Allele Fraction model
@Test
public void obviousCallsTest() {
    // only the second state
    final List<Double> weights = Arrays.asList(0.5, 0.5);
    final List<Double> minorAlleleFractions = Arrays.asList(0.1, 0.5);
    final double memoryLength = 1e3;
    final AlleleFractionHMM model = new AlleleFractionHMM(minorAlleleFractions, weights, memoryLength, AllelicPanelOfNormals.EMPTY_PON, NO_BIAS_OR_OUTLIERS_PARAMS);
    final Random random = new Random(13);
    final int chainLength = 10000;
    final List<SimpleInterval> positions = new ArrayList<>();
    final List<AllelicCount> data = new ArrayList<>();
    int position = 1;
    for (int n = 0; n < chainLength; n++) {
        position += random.nextInt((int) (2 * memoryLength)) + (int) memoryLength;
        final SimpleInterval interval = new SimpleInterval("chr1", position, position);
        positions.add(interval);
        if (n < 2500) {
            // minor fraction = 0.1
            data.add(new AllelicCount(interval, 10, 90));
        } else if (n < 5000) {
            // minor fraction = 0.5
            data.add(new AllelicCount(interval, 50, 50));
        } else if (n < 7500) {
            // minor fraction = 0.1
            data.add(new AllelicCount(interval, 90, 10));
        } else {
            // minor fraction = 0.5
            data.add(new AllelicCount(interval, 40, 60));
        }
    }
    final List<Integer> states = ViterbiAlgorithm.apply(data, positions, model);
    for (int pos = 0; pos < chainLength; pos++) {
        if (pos < 2500) {
            Assert.assertEquals((int) states.get(pos), 0);
        } else if (pos < 5000) {
            Assert.assertEquals((int) states.get(pos), 1);
        } else if (pos < 7500) {
            Assert.assertEquals((int) states.get(pos), 0);
        } else {
            Assert.assertEquals((int) states.get(pos), 1);
        }
    }
}

Example 32 with AllelicCount

use of org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCount in project gatk-protected by broadinstitute.

the class GetBayesianHetCoverageIntegrationTest method testNormalJob.

@Test
public void testNormalJob() {
    final File normalOutputFile = createTempFile("normal-test", ".tsv");
    Pulldown normalHetPulldownExpected, normalHetPulldownResult;
    /* test 1: tight calling stringency */
    final String[] args_1 = { "-" + StandardArgumentDefinitions.REFERENCE_SHORT_NAME, REF_FILE.getAbsolutePath(), "-" + ExomeStandardArgumentDefinitions.SNP_FILE_SHORT_NAME, SNP_FILE.getAbsolutePath(), "-" + ExomeStandardArgumentDefinitions.NORMAL_BAM_FILE_SHORT_NAME, NORMAL_BAM_FILE.getAbsolutePath(), "-" + ExomeStandardArgumentDefinitions.NORMAL_ALLELIC_COUNTS_FILE_SHORT_NAME, normalOutputFile.getAbsolutePath(), "-" + GetBayesianHetCoverage.READ_DEPTH_THRESHOLD_SHORT_NAME, Integer.toString(10), "-" + GetBayesianHetCoverage.HET_CALLING_STRINGENCY_SHORT_NAME, Double.toString(15.0) };
    runCommandLine(args_1);
    normalHetPulldownExpected = new Pulldown(normalHeader);
    normalHetPulldownExpected.add(new AllelicCount(new SimpleInterval("1", 14630, 14630), 9, 8, Nucleotide.T, Nucleotide.G, 17, 39.99));
    normalHetPulldownResult = new Pulldown(normalOutputFile, normalHeader);
    Assert.assertEquals(normalHetPulldownExpected, normalHetPulldownResult);
    /* test 2: loose calling stringency */
    final String[] args_2 = { "-" + StandardArgumentDefinitions.REFERENCE_SHORT_NAME, REF_FILE.getAbsolutePath(), "-" + ExomeStandardArgumentDefinitions.SNP_FILE_SHORT_NAME, SNP_FILE.getAbsolutePath(), "-" + ExomeStandardArgumentDefinitions.NORMAL_BAM_FILE_SHORT_NAME, NORMAL_BAM_FILE.getAbsolutePath(), "-" + ExomeStandardArgumentDefinitions.NORMAL_ALLELIC_COUNTS_FILE_SHORT_NAME, normalOutputFile.getAbsolutePath(), "-" + GetBayesianHetCoverage.READ_DEPTH_THRESHOLD_SHORT_NAME, Integer.toString(10), "-" + GetBayesianHetCoverage.HET_CALLING_STRINGENCY_SHORT_NAME, Double.toString(2.0) };
    runCommandLine(args_2);
    normalHetPulldownExpected = new Pulldown(normalHeader);
    normalHetPulldownExpected.add(new AllelicCount(new SimpleInterval("1", 11522, 11522), 7, 4, Nucleotide.G, Nucleotide.A, 11, 18.60));
    normalHetPulldownExpected.add(new AllelicCount(new SimpleInterval("1", 12098, 12098), 8, 6, Nucleotide.G, Nucleotide.T, 14, 29.29));
    normalHetPulldownExpected.add(new AllelicCount(new SimpleInterval("1", 14630, 14630), 9, 8, Nucleotide.T, Nucleotide.G, 17, 39.99));
    normalHetPulldownExpected.add(new AllelicCount(new SimpleInterval("2", 14689, 14689), 6, 9, Nucleotide.T, Nucleotide.G, 15, 28.60));
    normalHetPulldownExpected.add(new AllelicCount(new SimpleInterval("2", 14982, 14982), 6, 5, Nucleotide.G, Nucleotide.C, 11, 24.99));
    normalHetPulldownResult = new Pulldown(normalOutputFile, normalHeader);
    Assert.assertEquals(normalHetPulldownExpected, normalHetPulldownResult);
}

Example 33 with AllelicCount

use of org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCount in project gatk-protected by broadinstitute.

the class GetHetCoverageIntegrationTest method initHeaders.

@BeforeClass
public void initHeaders() throws IOException {
    try (final SamReader normalBamReader = SamReaderFactory.makeDefault().open(NORMAL_BAM_FILE);
        final SamReader tumorBamReader = SamReaderFactory.makeDefault().open(TUMOR_BAM_FILE)) {
        normalHeader = normalBamReader.getFileHeader();
        tumorHeader = tumorBamReader.getFileHeader();
        normalHetPulldownExpected = new Pulldown(normalHeader);
        normalHetPulldownExpected.add(new AllelicCount(new SimpleInterval("1", 11522, 11522), 7, 4));
        normalHetPulldownExpected.add(new AllelicCount(new SimpleInterval("1", 12098, 12098), 8, 6));
        normalHetPulldownExpected.add(new AllelicCount(new SimpleInterval("1", 14630, 14630), 9, 8));
        normalHetPulldownExpected.add(new AllelicCount(new SimpleInterval("2", 14689, 14689), 6, 9));
        normalHetPulldownExpected.add(new AllelicCount(new SimpleInterval("2", 14982, 14982), 6, 5));
        tumorHetPulldownExpected = new Pulldown(tumorHeader);
        tumorHetPulldownExpected.add(new AllelicCount(new SimpleInterval("1", 11522, 11522), 7, 4));
        tumorHetPulldownExpected.add(new AllelicCount(new SimpleInterval("1", 12098, 12098), 8, 6));
        tumorHetPulldownExpected.add(new AllelicCount(new SimpleInterval("1", 14630, 14630), 9, 8));
        tumorHetPulldownExpected.add(new AllelicCount(new SimpleInterval("2", 14689, 14689), 6, 9));
        tumorHetPulldownExpected.add(new AllelicCount(new SimpleInterval("2", 14982, 14982), 6, 5));
    }
}

Example 34 with AllelicCount

use of org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCount in project gatk-protected by broadinstitute.

the class GetBayesianHetCoverageIntegrationTest method testTumorJob.

@Test
public void testTumorJob() {
    final File tumorOutputFile = createTempFile("tumor-test", ".tsv");
    Pulldown tumorHetPulldownExpected, tumorHetPulldownResult;
    /* test 1: tight calling stringency */
    final String[] args_1 = { "-" + StandardArgumentDefinitions.REFERENCE_SHORT_NAME, REF_FILE.getAbsolutePath(), "-" + ExomeStandardArgumentDefinitions.SNP_FILE_SHORT_NAME, SNP_FILE.getAbsolutePath(), "-" + ExomeStandardArgumentDefinitions.TUMOR_BAM_FILE_SHORT_NAME, TUMOR_BAM_FILE.getAbsolutePath(), "-" + ExomeStandardArgumentDefinitions.TUMOR_ALLELIC_COUNTS_FILE_SHORT_NAME, tumorOutputFile.getAbsolutePath(), "-" + GetBayesianHetCoverage.READ_DEPTH_THRESHOLD_SHORT_NAME, Integer.toString(10), "-" + GetBayesianHetCoverage.HET_CALLING_STRINGENCY_SHORT_NAME, Double.toString(10.0) };
    runCommandLine(args_1);
    tumorHetPulldownExpected = new Pulldown(tumorHeader);
    tumorHetPulldownExpected.add(new AllelicCount(new SimpleInterval("1", 12098, 12098), 8, 6, Nucleotide.G, Nucleotide.T, 14, 24.85));
    tumorHetPulldownExpected.add(new AllelicCount(new SimpleInterval("1", 14630, 14630), 9, 8, Nucleotide.T, Nucleotide.G, 17, 34.03));
    tumorHetPulldownExpected.add(new AllelicCount(new SimpleInterval("2", 14689, 14689), 6, 9, Nucleotide.T, Nucleotide.G, 15, 24.23));
    tumorHetPulldownResult = new Pulldown(tumorOutputFile, tumorHeader);
    Assert.assertEquals(tumorHetPulldownExpected, tumorHetPulldownResult);
    /* test 2: tight calling stringency */
    final String[] args_2 = { "-" + StandardArgumentDefinitions.REFERENCE_SHORT_NAME, REF_FILE.getAbsolutePath(), "-" + ExomeStandardArgumentDefinitions.SNP_FILE_SHORT_NAME, SNP_FILE.getAbsolutePath(), "-" + ExomeStandardArgumentDefinitions.TUMOR_BAM_FILE_SHORT_NAME, TUMOR_BAM_FILE.getAbsolutePath(), "-" + ExomeStandardArgumentDefinitions.TUMOR_ALLELIC_COUNTS_FILE_SHORT_NAME, tumorOutputFile.getAbsolutePath(), "-" + GetBayesianHetCoverage.READ_DEPTH_THRESHOLD_SHORT_NAME, Integer.toString(10), "-" + GetBayesianHetCoverage.HET_CALLING_STRINGENCY_SHORT_NAME, Double.toString(5.0) };
    runCommandLine(args_2);
    tumorHetPulldownExpected = new Pulldown(tumorHeader);
    tumorHetPulldownExpected.add(new AllelicCount(new SimpleInterval("1", 11522, 11522), 7, 4, Nucleotide.G, Nucleotide.A, 11, 15.72));
    tumorHetPulldownExpected.add(new AllelicCount(new SimpleInterval("1", 12098, 12098), 8, 6, Nucleotide.G, Nucleotide.T, 14, 24.85));
    tumorHetPulldownExpected.add(new AllelicCount(new SimpleInterval("1", 14630, 14630), 9, 8, Nucleotide.T, Nucleotide.G, 17, 34.03));
    tumorHetPulldownExpected.add(new AllelicCount(new SimpleInterval("2", 14689, 14689), 6, 9, Nucleotide.T, Nucleotide.G, 15, 24.23));
    tumorHetPulldownExpected.add(new AllelicCount(new SimpleInterval("2", 14982, 14982), 6, 5, Nucleotide.G, Nucleotide.C, 11, 21.23));
    tumorHetPulldownResult = new Pulldown(tumorOutputFile, tumorHeader);
    Assert.assertEquals(tumorHetPulldownExpected, tumorHetPulldownResult);
}

Example 35 with AllelicCount

use of org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCount in project gatk by broadinstitute.

the class BayesianHetPulldownCalculatorUnitTest method testGetHetPulldown.

@Test
public void testGetHetPulldown() {
    Pulldown testPulldown, expectedPulldown;
    File tempFile;
    /* write Pulldown to file while checking the results */
    try {
        /* test 1: normal, loose threshold */
        testPulldown = calculator.getHetPulldown(NORMAL_BAM_FILE, 2);
        tempFile = File.createTempFile("testPulldownNormalLoose", ".txt");
        testPulldown.write(tempFile, AllelicCountTableColumn.AllelicCountTableVerbosity.FULL);
        expectedPulldown = new Pulldown(normalHeader);
        expectedPulldown.add(new AllelicCount(new SimpleInterval("1", 11522, 11522), 7, 4, Nucleotide.G, Nucleotide.A, 11, 18.38));
        expectedPulldown.add(new AllelicCount(new SimpleInterval("1", 12098, 12098), 8, 6, Nucleotide.G, Nucleotide.T, 14, 28.84));
        expectedPulldown.add(new AllelicCount(new SimpleInterval("1", 14630, 14630), 9, 8, Nucleotide.T, Nucleotide.G, 17, 39.39));
        expectedPulldown.add(new AllelicCount(new SimpleInterval("2", 14689, 14689), 6, 9, Nucleotide.T, Nucleotide.G, 15, 28.23));
        expectedPulldown.add(new AllelicCount(new SimpleInterval("2", 14982, 14982), 6, 5, Nucleotide.G, Nucleotide.C, 11, 24.54));
        Assert.assertEquals(new Pulldown(tempFile, normalHeader), expectedPulldown);
        /* test 2: normal, tight threshold */
        testPulldown = calculator.getHetPulldown(NORMAL_BAM_FILE, 12);
        tempFile = File.createTempFile("testPulldownNormalTight", ".txt");
        testPulldown.write(tempFile, AllelicCountTableColumn.AllelicCountTableVerbosity.FULL);
        expectedPulldown = new Pulldown(normalHeader);
        expectedPulldown.add(new AllelicCount(new SimpleInterval("1", 12098, 12098), 8, 6, Nucleotide.G, Nucleotide.T, 14, 28.84));
        expectedPulldown.add(new AllelicCount(new SimpleInterval("1", 14630, 14630), 9, 8, Nucleotide.T, Nucleotide.G, 17, 39.39));
        expectedPulldown.add(new AllelicCount(new SimpleInterval("2", 14689, 14689), 6, 9, Nucleotide.T, Nucleotide.G, 15, 28.23));
        Assert.assertEquals(new Pulldown(tempFile, normalHeader), expectedPulldown);
        /* test 3: tumor, loose threshold */
        testPulldown = calculator.getHetPulldown(TUMOR_BAM_FILE, 2);
        tempFile = File.createTempFile("testPulldownTumorLoose", ".txt");
        testPulldown.write(tempFile, AllelicCountTableColumn.AllelicCountTableVerbosity.FULL);
        expectedPulldown = new Pulldown(tumorHeader);
        expectedPulldown.add(new AllelicCount(new SimpleInterval("1", 11522, 11522), 7, 4, Nucleotide.G, Nucleotide.A, 11, 15.63));
        expectedPulldown.add(new AllelicCount(new SimpleInterval("1", 12098, 12098), 8, 6, Nucleotide.G, Nucleotide.T, 14, 24.72));
        expectedPulldown.add(new AllelicCount(new SimpleInterval("1", 14630, 14630), 9, 8, Nucleotide.T, Nucleotide.G, 17, 33.89));
        expectedPulldown.add(new AllelicCount(new SimpleInterval("2", 14689, 14689), 6, 9, Nucleotide.T, Nucleotide.G, 15, 24.11));
        expectedPulldown.add(new AllelicCount(new SimpleInterval("2", 14982, 14982), 6, 5, Nucleotide.G, Nucleotide.C, 11, 21.10));
        Assert.assertEquals(new Pulldown(tempFile, normalHeader), expectedPulldown);
        /* test 4: tumor, tight threshold */
        testPulldown = calculator.getHetPulldown(TUMOR_BAM_FILE, 12);
        tempFile = File.createTempFile("testPulldownTumorTight", ".txt");
        testPulldown.write(tempFile, AllelicCountTableColumn.AllelicCountTableVerbosity.FULL);
        expectedPulldown = new Pulldown(tumorHeader);
        expectedPulldown.add(new AllelicCount(new SimpleInterval("1", 14630, 14630), 9, 8, Nucleotide.T, Nucleotide.G, 17, 33.89));
        Assert.assertEquals(new Pulldown(tempFile, normalHeader), expectedPulldown);
    } catch (final IOException e) {
        throw new UserException.CouldNotCreateOutputFile("Could not write pulldown to to file.", e);
    }
}