Make error informative for non-diploid family likelihoods

This commit is contained in:
Ron Levine 2017-07-18 23:07:27 -04:00
parent e406bb76ea
commit 57588a5661
3 changed files with 26 additions and 3 deletions

View File

@ -185,6 +185,8 @@ import java.util.*;
* -o output.withPosteriors.vcf * -o output.withPosteriors.vcf
* </pre> * </pre>
* *
* <h3>Caveat</h3>
* <p>If applying family priors, only diploid family genotypes are supported</p>
*/ */
@DocumentedGATKFeature( groupName = HelpConstants.DOCS_CAT_VARDISC, extraDocs = {CommandLineGATK.class} ) @DocumentedGATKFeature( groupName = HelpConstants.DOCS_CAT_VARDISC, extraDocs = {CommandLineGATK.class} )
public class CalculateGenotypePosteriors extends RodWalker<Integer,Integer> { public class CalculateGenotypePosteriors extends RodWalker<Integer,Integer> {

View File

@ -51,6 +51,7 @@
package org.broadinstitute.gatk.tools.walkers.variantutils; package org.broadinstitute.gatk.tools.walkers.variantutils;
import htsjdk.variant.vcf.VCFConstants;
import org.apache.log4j.Logger; import org.apache.log4j.Logger;
import org.broadinstitute.gatk.engine.samples.Sample; import org.broadinstitute.gatk.engine.samples.Sample;
import org.broadinstitute.gatk.utils.MathUtils; import org.broadinstitute.gatk.utils.MathUtils;
@ -263,7 +264,7 @@ public class FamilyLikelihoodsUtils {
updateFamilyGenotypes(vc, mother, father, child, trioGenotypes); updateFamilyGenotypes(vc, mother, father, child, trioGenotypes);
//replace uses sample names to match genotypes, so order doesn't matter //replace uses sample names to match genotypes, so order doesn't matter
if (trioGenotypes.size() > 0) { if (!trioGenotypes.isEmpty()) {
genotypesContext.replace(trioGenotypes.get(0)); genotypesContext.replace(trioGenotypes.get(0));
genotypesContext.replace(trioGenotypes.get(1)); genotypesContext.replace(trioGenotypes.get(1));
genotypesContext.replace(trioGenotypes.get(2)); genotypesContext.replace(trioGenotypes.get(2));
@ -440,7 +441,7 @@ public class FamilyLikelihoodsUtils {
//In case of null, unavailable or no call, all likelihoods are log10(1/3) //In case of null, unavailable or no call, all likelihoods are log10(1/3)
else if(genotype == null || !hasCalledGT(genotype.getType()) || genotype.getLikelihoods() == null){ else if(genotype == null || !hasCalledGT(genotype.getType()) || genotype.getLikelihoods() == null){
likelihoods = new double[3]; likelihoods = new double[NUM_CALLED_GENOTYPETYPES];
likelihoods[0] = LOG10_OF_ONE_THIRD; likelihoods[0] = LOG10_OF_ONE_THIRD;
likelihoods[1] = LOG10_OF_ONE_THIRD; likelihoods[1] = LOG10_OF_ONE_THIRD;
likelihoods[2] = LOG10_OF_ONE_THIRD; likelihoods[2] = LOG10_OF_ONE_THIRD;
@ -450,6 +451,13 @@ public class FamilyLikelihoodsUtils {
else else
likelihoods = GeneralUtils.normalizeFromLog10(genotype.getLikelihoods().getAsVector(),true,true); likelihoods = GeneralUtils.normalizeFromLog10(genotype.getLikelihoods().getAsVector(),true,true);
if (likelihoods.length != NUM_CALLED_GENOTYPETYPES) {
final String key = genotype.hasExtendedAttribute(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY) ?
GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY : VCFConstants.GENOTYPE_PL_KEY;
throw new UserException(genotype + " has " + likelihoods.length + " " + key + " values, should be " + NUM_CALLED_GENOTYPETYPES +
" since only the diploid case is supported when applying family priors.");
}
likelihoodsMap.put(GenotypeType.HOM_REF,likelihoods[genotypeTypeToValue(GenotypeType.HOM_REF)]); likelihoodsMap.put(GenotypeType.HOM_REF,likelihoods[genotypeTypeToValue(GenotypeType.HOM_REF)]);
likelihoodsMap.put(GenotypeType.HET,likelihoods[genotypeTypeToValue(GenotypeType.HET)]); likelihoodsMap.put(GenotypeType.HET,likelihoods[genotypeTypeToValue(GenotypeType.HET)]);
likelihoodsMap.put(GenotypeType.HOM_VAR, likelihoods[genotypeTypeToValue(GenotypeType.HOM_VAR)]); likelihoodsMap.put(GenotypeType.HOM_VAR, likelihoods[genotypeTypeToValue(GenotypeType.HOM_VAR)]);

View File

@ -52,6 +52,7 @@
package org.broadinstitute.gatk.tools.walkers.variantutils; package org.broadinstitute.gatk.tools.walkers.variantutils;
import org.broadinstitute.gatk.engine.walkers.WalkerTest; import org.broadinstitute.gatk.engine.walkers.WalkerTest;
import org.broadinstitute.gatk.utils.exceptions.UserException;
import org.testng.annotations.Test; import org.testng.annotations.Test;
import java.util.Arrays; import java.util.Arrays;
@ -61,6 +62,7 @@ public class CalculateGenotypePosteriorsIntegrationTest extends WalkerTest {
private static String CEUtrioFamilyFile = privateTestDir + "CEUtrio.ped"; private static String CEUtrioFamilyFile = privateTestDir + "CEUtrio.ped";
private static String CEUtrioTest = privateTestDir + "CEUtrioTest.vcf"; private static String CEUtrioTest = privateTestDir + "CEUtrioTest.vcf";
private static String CEUtrioPopPriorsTest = privateTestDir + "CEUtrioPopPriorsTest.vcf"; private static String CEUtrioPopPriorsTest = privateTestDir + "CEUtrioPopPriorsTest.vcf";
private static String CEUtrioMixedPloidyTest = privateTestDir + "CEUtrioMixedPloidy.vcf";
private static String threeMemberNonTrioFamilyFile = privateTestDir + "threeMemberNonTrio.ped"; private static String threeMemberNonTrioFamilyFile = privateTestDir + "threeMemberNonTrio.ped";
private static String getThreeMemberNonTrioTest = privateTestDir + "threeMemberNonTrioTest.vcf"; private static String getThreeMemberNonTrioTest = privateTestDir + "threeMemberNonTrioTest.vcf";
@ -135,5 +137,16 @@ public class CalculateGenotypePosteriorsIntegrationTest extends WalkerTest {
executeTest("testFamilyPriors", spec); executeTest("testFamilyPriors", spec);
} }
@Test(enabled = true)
public void testFamilyPriorsMixedPloidy() {
WalkerTestSpec spec = new WalkerTestSpec(
"-T CalculateGenotypePosteriors --no_cmdline_in_header " +
" -o %s" +
" -R " + b37KGReference +
" -ped " + CEUtrioFamilyFile +
" -V " + CEUtrioMixedPloidyTest,
1,
UserException.class);
executeTest("testFamilyPriorsMixedPloidy", spec);
}
} }