Skip to content

Commit

Permalink
Fix CGP tests -- GQ0 posterior is still hom-ref, but we can put that …
Browse files Browse the repository at this point in the history 
…in the refactoring PR
  • Loading branch information
@ldgauthier
ldgauthier committed Apr 5, 2024
1 parent 2353ddc commit 0fd8ac2
Show file tree
Hide file tree
Showing 5 changed files with 74 additions and 32 deletions.
6 changes: 5 additions & 1 deletion ...org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtils.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
import org.broadinstitute.hellbender.exceptions.UserException;
import org.broadinstitute.hellbender.tools.walkers.ReferenceConfidenceVariantContextMerger;
import org.broadinstitute.hellbender.tools.walkers.genotyper.GenotypeAssignmentMethod;
import org.broadinstitute.hellbender.tools.walkers.genotyper.GenotypesCache;
import org.broadinstitute.hellbender.utils.MathUtils;
import org.broadinstitute.hellbender.utils.SimpleInterval;
import org.broadinstitute.hellbender.utils.Utils;
Expand Down Expand Up @@ -138,8 +139,11 @@ else if (a.isSymbolic()) { //use the greater of the priors for non-ref

builder.attribute(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY,
Utils.listFromPrimitives(GenotypeLikelihoods.fromLog10Likelihoods(posteriors.get(genoIdx)).getAsPLs()));
//TODO: refactor me -- this is to make up for the GGVCFs changes fixing GQ0 hom-refs
if (!builder.makeWithShallowCopy().hasPL()) {
builder.log10PError(GenotypeLikelihoods.getGQLog10FromLikelihoods(MathUtils.maxElementIndex(posteriors.get(genoIdx)), posteriors.get(genoIdx)));
final int maxPosteriorIndex = MathUtils.maxElementIndex(posteriors.get(genoIdx));
builder.log10PError(GenotypeLikelihoods.getGQLog10FromLikelihoods(maxPosteriorIndex, posteriors.get(genoIdx)));
builder.alleles(GenotypesCache.get(vc1.getGenotype(genoIdx).getPloidy(), maxPosteriorIndex).asAlleleList(vc1.getAlleles()));
}
}
newContext.add(builder.make());
Expand Down
23 changes: 4 additions & 19 deletions src/test/java/org/broadinstitute/hellbender/tools/walkers/GenotypeGVCFsIntegrationTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,21 +77,6 @@ public class GenotypeGVCFsIntegrationTest extends CommandLineProgramTest {

private static final String ALLELE_SPECIFIC_DIRECTORY = toolsTestDir + "walkers/annotator/allelespecific";

private static <T> void assertForEachElementInLists(final List<T> actual, final List<T> expected, final BiConsumer<T, T> assertion) {
Assert.assertEquals(actual.size(), expected.size(), "different number of elements in lists:\n"
+ actual.stream().map(Object::toString).collect(Collectors.joining("\n","actual:\n","\n"))
+ expected.stream().map(Object::toString).collect(Collectors.joining("\n","expected:\n","\n")));
for (int i = 0; i < actual.size(); i++) {
assertion.accept(actual.get(i), expected.get(i));
}
}

private static <T> void assertCountForEachElementInList(final List<T> actual, Integer count, BiConsumer<T, Integer> assertion) {
for (int i = 0; i < actual.size(); i++) {
assertion.accept(actual.get(i), count);
}
}

@DataProvider(name = "gvcfsToGenotype")
public Object[][] gvcfsToGenotype() {
return new Object[][]{
Expand Down Expand Up @@ -419,7 +404,7 @@ private void runAndCheckGenomicsDBOutput(final ArgumentsBuilder args, final File
// runs in 0.06 minutes with no input intervals specified
final List<VariantContext> expectedVC = VariantContextTestUtils.getVariantContexts(expected);
final List<VariantContext> actualVC = VariantContextTestUtils.getVariantContexts(output);
assertForEachElementInLists(actualVC, expectedVC, VariantContextTestUtils::assertVariantContextsHaveSameGenotypes);
VariantContextTestUtils.assertForEachElementInLists(actualVC, expectedVC, VariantContextTestUtils::assertVariantContextsHaveSameGenotypes);
}

@Test(dataProvider = "getGVCFsForGenomicsDBOverMultipleIntervals")
Expand Down Expand Up @@ -534,7 +519,7 @@ private void runGenotypeGVCFSAndAssertComparison(String input, File expected, Li
final List<VariantContext> expectedVC = VariantContextTestUtils.getVariantContexts(expected);
final List<VariantContext> actualVC = VariantContextTestUtils.getVariantContexts(output);
try {
assertForEachElementInLists(actualVC, expectedVC, assertion);
VariantContextTestUtils.assertForEachElementInLists(actualVC, expectedVC, assertion);
} catch (final AssertionError error) {
throw error;
}
Expand All @@ -548,7 +533,7 @@ private void runGenotypeGVCFSAndAssertCount(final String input, final List<Strin

final List<VariantContext> actualVC = VariantContextTestUtils.getVariantContexts(output);
Assert.assertTrue(actualVC.size() > 0);
assertCountForEachElementInList(actualVC, count, conditionOnCount);
VariantContextTestUtils.assertCountForEachElementInList(actualVC, count, conditionOnCount);
}

@Test
Expand Down Expand Up @@ -647,7 +632,7 @@ public void testGenotypingForSomaticGVCFs() throws IOException{
final File expectedFile = getTestFile("threeSamples.MT.vcf");
final List<VariantContext> expected = VariantContextTestUtils.getVariantContexts(expectedFile);
final VCFHeader header = VCFHeaderReader.readHeaderFrom(new SeekablePathStream(IOUtils.getPath(expectedFile.getAbsolutePath())));
assertForEachElementInLists(results, expected,
VariantContextTestUtils.assertForEachElementInLists(results, expected,
(a, e) -> VariantContextTestUtils.assertVariantContextsAreEqualAlleleOrderIndependent(a, e, ATTRIBUTES_TO_IGNORE, ATTRIBUTES_WITH_JITTER, header));

}
Expand Down
59 changes: 48 additions & 11 deletions ...ute/hellbender/tools/walkers/variantutils/CalculateGenotypePosteriorsIntegrationTest.java
View file
Original file line number Diff line number Diff line change
@@ -1,17 +1,24 @@
package org.broadinstitute.hellbender.tools.walkers.variantutils;

import htsjdk.samtools.util.BlockCompressedInputStream;
import htsjdk.variant.variantcontext.VariantContext;
import org.broadinstitute.hellbender.CommandLineProgramTest;
import org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions;
import org.broadinstitute.hellbender.exceptions.UserException;
import org.broadinstitute.hellbender.testutils.ArgumentsBuilder;
import org.broadinstitute.hellbender.testutils.IntegrationTestSpec;
import org.broadinstitute.hellbender.testutils.VariantContextTestUtils;
import org.broadinstitute.hellbender.tools.walkers.annotator.RMSMappingQuality;
import org.broadinstitute.hellbender.utils.Utils;
import org.testng.Assert;
import org.testng.annotations.DataProvider;
import org.testng.annotations.Test;

import java.io.*;
import java.util.Collections;
import java.util.List;

import static org.broadinstitute.hellbender.testutils.VariantContextTestUtils.assertForEachElementInLists;

public final class CalculateGenotypePosteriorsIntegrationTest extends CommandLineProgramTest {

Expand Down Expand Up @@ -59,16 +66,25 @@ public void testNumRefWithPanel() throws IOException {
//use the first 20 variants to save time; they have a nice range of AC from 4 to over 4000
//three variants have GTs with zero depth and PLs=[1,0,0], which should get PPs
public void testUsingDiscoveredAF() throws IOException {
final IntegrationTestSpec spec = new IntegrationTestSpec(
" -O %s" +
" -R " + b37_reference_20_21 + //NOTE: we need a reference for -L
//" -L 20:10,000,000-10,001,432" +
" -L 20:10001432" +
" -V " + largeDir + "1000G.phase3.broad.withGenotypes.chr20.10100000.vcf" +
" --" + StandardArgumentDefinitions.ADD_OUTPUT_VCF_COMMANDLINE +" false",
Collections.singletonList(largeDir + "CalculateGenotypePosteriors/expectedCGP_testUsingDiscoveredAF.vcf")
);
spec.executeTest("testUsingDiscoveredAF", this);
final File output = createTempFile("cgp_output", ".vcf");
final ArgumentsBuilder args = new ArgumentsBuilder();
args.addReference(new File(b37_reference_20_21)) //note we need a ref for -L arg
.add("V", largeDir + "1000G.phase3.broad.withGenotypes.chr20.10100000.vcf")
.add("L", "20:10,000,000-10,001,432")
.add(StandardArgumentDefinitions.ADD_OUTPUT_VCF_COMMANDLINE, false)
.addOutput(output);
runCommandLine(args);

Assert.assertTrue(output.exists());

File expected = new File(largeDir + "CalculateGenotypePosteriors/expectedCGP_testUsingDiscoveredAF.vcf");
final List<VariantContext> expectedVC = VariantContextTestUtils.getVariantContexts(expected);
final List<VariantContext> actualVC = VariantContextTestUtils.getVariantContexts(output);
try {
assertForEachElementInLists(actualVC, expectedVC, VariantContextTestUtils::assertVariantContextsHaveSameGenotypes);
} catch (final AssertionError error) {
throw error;
}
}

@Test
Expand All @@ -84,7 +100,28 @@ public void testMissingPriors() throws IOException {
" --" + StandardArgumentDefinitions.ADD_OUTPUT_VCF_COMMANDLINE +" false",
Collections.singletonList(largeDir + "CalculateGenotypePosteriors/expectedCGP_testMissingPriors.vcf")
);
spec.executeTest("testMissingPriors", this);
//spec.executeTest("testMissingPriors", this);

final File output = createTempFile("cgp_output", ".vcf");
final ArgumentsBuilder args = new ArgumentsBuilder();
args.addReference(new File(b37_reference_20_21)) //note we need a ref for -L arg
.addFlag("discovered-allele-count-priors-off")
.add("V", largeDir + "1000G.phase3.broad.withGenotypes.chr20.10100000.vcf")
.add("L", "20:10,000,000-10,001,432")
.add(StandardArgumentDefinitions.ADD_OUTPUT_VCF_COMMANDLINE, false)
.addOutput(output);
runCommandLine(args);

Assert.assertTrue(output.exists());

File expected = new File(largeDir + "CalculateGenotypePosteriors/expectedCGP_testMissingPriors.vcf");
final List<VariantContext> expectedVC = VariantContextTestUtils.getVariantContexts(expected);
final List<VariantContext> actualVC = VariantContextTestUtils.getVariantContexts(output);
try {
assertForEachElementInLists(actualVC, expectedVC, VariantContextTestUtils::assertVariantContextsHaveSameGenotypes);
} catch (final AssertionError error) {
throw error;
}
}

@Test
Expand Down
2 changes: 1 addition & 1 deletion src/test/resources/large/CalculateGenotypePosteriors/expectedCGP_testMissingPriors.vcf
View file

Large diffs are not rendered by default.

16 changes: 16 additions & 0 deletions src/testUtils/java/org/broadinstitute/hellbender/testutils/VariantContextTestUtils.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@
import java.io.IOException;
import java.io.PrintStream;
import java.util.*;
import java.util.function.BiConsumer;
import java.util.stream.Collectors;
import java.util.stream.Stream;
import java.util.stream.StreamSupport;
Expand Down Expand Up @@ -781,4 +782,19 @@ public static String keyForVariant(final VariantContext variant) {
return String.format("%s:%d-%d %s, %s", variant.getContig(), variant.getStart(), end, variant.getReference(),
variant.getAlternateAlleles().stream().map(Allele::getDisplayString).sorted().collect(Collectors.toList()));
}

public static <T> void assertForEachElementInLists(final List<T> actual, final List<T> expected, final BiConsumer<T, T> assertion) {
Assert.assertEquals(actual.size(), expected.size(), "different number of elements in lists:\n"
+ actual.stream().map(Object::toString).collect(Collectors.joining("\n","actual:\n","\n"))
+ expected.stream().map(Object::toString).collect(Collectors.joining("\n","expected:\n","\n")));
for (int i = 0; i < actual.size(); i++) {
assertion.accept(actual.get(i), expected.get(i));
}
}

public static <T> void assertCountForEachElementInList(final List<T> actual, Integer count, BiConsumer<T, Integer> assertion) {
for (int i = 0; i < actual.size(); i++) {
assertion.accept(actual.get(i), count);
}
}
}

0 comments on commit 0fd8ac2

Please sign in to comment.