diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/PedigreeAnnotation.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/PedigreeAnnotation.java index 2621f2ebbc9..b64741f5191 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/PedigreeAnnotation.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/PedigreeAnnotation.java @@ -16,7 +16,7 @@ * A common interface for handling annotations that require pedigree file information either in the form of explicitly * selected founderIDs or in the form of an imported pedigreeFile. * - * In order to use the the behavior, simply extend Pedigree annotation and access its constructors, then call + * In order to use the behavior, simply extend Pedigree annotation and access its constructors, then call * getFounderGenotypes() to extract only the genotypes corresponding to requested founder samples or that appear as founders * in the provided pedigree file. If no founderIDs or pedigreeFiles are present, then it defaults to returning all genotypes. * @@ -103,6 +103,26 @@ void validateArguments(Collection founderIds, GATKPath pedigreeFile) { } } + /** + * Warning generator for when a pedigree file is required and founderIDs cannot be used for this annotation. + * @param founderIds + * @param pedigreeFile + */ + protected String validateArgumentsWhenPedigreeRequired(Collection founderIds, GATKPath pedigreeFile) { + if (pedigreeFile == null) { + if ((founderIds != null && !founderIds.isEmpty())) { + return "PossibleDenovo annotation will not be calculated, must provide a valid PED file (-ped). Founder-id arguments cannot be used for this annotation"; + } else { + return "PossibleDenovo Annotation will not be calculated, must provide a valid PED file (-ped) from the command line."; + } + } else { + if ((founderIds != null && !founderIds.isEmpty())) { + return "PossibleDenovo annotation does not take founder-id arguments, trio information will be extracted only from the provided PED file"; + } + } + return null; + } + /** * This is a getter for the pedigree file, which could be null. * @@ -111,4 +131,19 @@ void validateArguments(Collection founderIds, GATKPath pedigreeFile) { public @Nullable GATKPath getPedigreeFile() { return pedigreeFile; } + + /** + * Helper function to check if the variant context has GQs for the trio + * @param vc variant context + * @param trio trio to check for GQs in the variant context + */ + protected static boolean contextHasTrioGQs(final VariantContext vc, final Trio trio) { + final String mom = trio.getMaternalID(); + final String dad = trio.getPaternalID(); + final String kid = trio.getChildID(); + + return (!mom.isEmpty() && vc.hasGenotype(mom) && vc.getGenotype(mom).hasGQ()) + && (!dad.isEmpty() && vc.hasGenotype(dad) && vc.getGenotype(dad).hasGQ()) + && (!kid.isEmpty() && vc.hasGenotype(kid) && vc.getGenotype(kid).hasGQ()); + } } diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/PossibleDeNovo.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/PossibleDeNovo.java index c03fc8549f6..6c453d8b5e5 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/PossibleDeNovo.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/PossibleDeNovo.java @@ -5,6 +5,7 @@ import htsjdk.variant.variantcontext.VariantContext; import org.apache.logging.log4j.LogManager; import org.apache.logging.log4j.Logger; +import org.broadinstitute.barclay.argparser.Argument; import org.broadinstitute.barclay.help.DocumentedFeature; import org.broadinstitute.hellbender.engine.GATKPath; import org.broadinstitute.hellbender.engine.ReferenceContext; @@ -35,15 +36,30 @@ */ @DocumentedFeature(groupName=HelpConstants.DOC_CAT_ANNOTATORS, groupSummary=HelpConstants.DOC_CAT_ANNOTATORS_SUMMARY, summary="Existence of a de novo mutation in at least one of the given families (hiConfDeNovo, loConfDeNovo)") public final class PossibleDeNovo extends PedigreeAnnotation implements InfoFieldAnnotation { + public static final String DENOVO_PARENT_GQ_THRESHOLD = "denovo-parent-gq-threshold"; + public static final String DENOVO_DEPTH_THRESHOLD = "denovo-depth-threshold"; protected final Logger warning = LogManager.getLogger(this.getClass()); private final MendelianViolation mendelianViolation; private Set trios; + @Argument(fullName = DENOVO_PARENT_GQ_THRESHOLD, doc = "Minimum genotype quality for parents to be considered for de novo calculation (separate from GQ thershold for full trio).", optional = true) + public int parentGQThreshold = hi_GQ_threshold; + + @Argument(fullName = DENOVO_DEPTH_THRESHOLD, doc = "Minimum depth (DP) for all trio members to be considered for de novo calculation.", optional = true) + public int depthThreshold = 0; + @VisibleForTesting - public PossibleDeNovo(final Set trios, final double minGenotypeQualityP) { + public PossibleDeNovo(final Set trios, final double minGenotypeQualityP, final int parentGQThreshold, final int depthThreshold) { super((Set) null); this.trios = Collections.unmodifiableSet(new LinkedHashSet<>(trios)); mendelianViolation = new MendelianViolation(minGenotypeQualityP); + this.parentGQThreshold = parentGQThreshold; + this.depthThreshold = depthThreshold; + } + + @VisibleForTesting + public PossibleDeNovo(final Set trios, final double minGenotypeQualityP) { + this(trios, minGenotypeQualityP, hi_GQ_threshold, 0); } public PossibleDeNovo(final GATKPath pedigreeFile){ @@ -58,19 +74,13 @@ public PossibleDeNovo(){ @Override void validateArguments(Collection founderIds, GATKPath pedigreeFile) { - if (pedigreeFile == null) { - if ((founderIds != null && !founderIds.isEmpty())) { - warning.warn("PossibleDenovo annotation will not be calculated, must provide a valid PED file (-ped). Founder-id arguments cannot be used for this annotation"); - } else { - warning.warn("PossibleDenovo Annotation will not be calculated, must provide a valid PED file (-ped) from the command line."); - } - } else { - if ((founderIds != null && !founderIds.isEmpty())) { - warning.warn("PossibleDenovo annotation does not take founder-id arguments, trio information will be extracted only from the provided PED file"); - } + String warningString = validateArgumentsWhenPedigreeRequired(founderIds, pedigreeFile); + if (warningString != null) { + warning.warn(warningString); } } + // Static thresholds for the denovo calculation public final static double DEFAULT_MIN_GENOTYPE_QUALITY_P = 0; // TODO should this be exposed as a command line argument? private static final int hi_GQ_threshold = 20; //WARNING - If you change this value, update the description in GATKVCFHeaderLines @@ -105,7 +115,7 @@ public Map annotate(final ReferenceContext ref, final List lowConfDeNovoChildren = new ArrayList<>(); for (final Trio trio : trioSet) { if (vc.isBiallelic() && - PossibleDeNovo.contextHasTrioGQs(vc, trio) && + contextHasTrioGQs(vc, trio) && mendelianViolation.isViolation(trio.getMother(), trio.getFather(), trio.getChild(), vc) && mendelianViolation.getParentsRefRefChildHet() > 0) { @@ -113,7 +123,12 @@ public Map annotate(final ReferenceContext ref, final int momGQ = vc.getGenotype(trio.getMaternalID()).getGQ(); final int dadGQ = vc.getGenotype(trio.getPaternalID()).getGQ(); - if (childGQ >= hi_GQ_threshold && momGQ >= hi_GQ_threshold && dadGQ >= hi_GQ_threshold) { + final int childDP = vc.getGenotype(trio.getChildID()).getDP(); + final int momDP = vc.getGenotype(trio.getMaternalID()).getDP(); + final int dadDP = vc.getGenotype(trio.getPaternalID()).getDP(); + + if (childGQ >= hi_GQ_threshold && momGQ >= parentGQThreshold && dadGQ >= parentGQThreshold && + childDP >= depthThreshold && momDP >= depthThreshold && dadDP >= depthThreshold) { highConfDeNovoChildren.add(trio.getChildID()); } else if (childGQ >= lo_GQ_threshold && momGQ > 0 && dadGQ > 0) { lowConfDeNovoChildren.add(trio.getChildID()); @@ -135,14 +150,4 @@ public Map annotate(final ReferenceContext ref, return attributeMap; } - private static boolean contextHasTrioGQs(final VariantContext vc, final Trio trio) { - final String mom = trio.getMaternalID(); - final String dad = trio.getPaternalID(); - final String kid = trio.getChildID(); - - return (!mom.isEmpty() && vc.hasGenotype(mom) && vc.getGenotype(mom).hasGQ()) - && (!dad.isEmpty() && vc.hasGenotype(dad) && vc.getGenotype(dad).hasGQ()) - && (!kid.isEmpty() && vc.hasGenotype(kid) && vc.getGenotype(kid).hasGQ()); - } - } diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/TransmittedSingleton.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/TransmittedSingleton.java new file mode 100644 index 00000000000..06932f5f2d0 --- /dev/null +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/TransmittedSingleton.java @@ -0,0 +1,144 @@ +package org.broadinstitute.hellbender.tools.walkers.annotator; + +import com.google.common.annotations.VisibleForTesting; +import htsjdk.variant.variantcontext.Allele; +import htsjdk.variant.variantcontext.VariantContext; +import htsjdk.variant.vcf.VCFConstants; +import org.apache.logging.log4j.LogManager; +import org.apache.logging.log4j.Logger; +import org.broadinstitute.barclay.argparser.ExperimentalFeature; +import org.broadinstitute.barclay.help.DocumentedFeature; +import org.broadinstitute.hellbender.engine.GATKPath; +import org.broadinstitute.hellbender.engine.ReferenceContext; +import org.broadinstitute.hellbender.utils.Utils; +import org.broadinstitute.hellbender.utils.genotyper.AlleleLikelihoods; +import org.broadinstitute.hellbender.utils.help.HelpConstants; +import org.broadinstitute.hellbender.utils.logging.OneShotLogger; +import org.broadinstitute.hellbender.utils.read.GATKRead; +import org.broadinstitute.hellbender.utils.samples.Trio; +import org.broadinstitute.hellbender.utils.variant.GATKVCFConstants; + +import java.util.*; + + +/** + * Existence of a transmitted or non-transmitted singleton in at least one of the given families + * + *

This annotation uses the genotype information from individuals in family trios to identify transmitted and non-transmitted singletons and the sample(s) in which they occur. + * Transmitted singletons occur at sites in a cohort where the allele count is two and these two alleles occur in one parent and the child of a trio. A non-transmitted singleton + * are sites with an allele count of one and this one allele occurs in a parent, but not the child of a trio. In both cases the other parent must have a high quality hom ref call. + * + *

Caveats

+ *
    + *
  • The calculation assumes that the organism is diploid.
    • + *
  • This annotation requires a valid pedigree file.
    • + *
  • This annotation is only valid for trios, not quads or other more complicated family structures.
    • + *
  • Only reports possible singletons for families where each of the three samples has high GQ (>20) and high depth (>20)
    • + *
  • Only reports possible singletons at sites with a Call Rate greater than 90% (meaning less than 10% of the samples at the given site were no-calls)
    • + *
+ */ +@DocumentedFeature(groupName= HelpConstants.DOC_CAT_ANNOTATORS, groupSummary=HelpConstants.DOC_CAT_ANNOTATORS_SUMMARY, summary="Existence of a transmitted (or non-transmitted) singleton in at least one of the given families (transmittedSingleton, nonTransmittedSingleton)") +public final class TransmittedSingleton extends PedigreeAnnotation implements InfoFieldAnnotation { + protected final Logger warning = LogManager.getLogger(this.getClass()); + private final OneShotLogger oneShotLogger = new OneShotLogger(this.getClass()); + private Set trios; + private final int HI_GQ_THRESHOLD = 20; + private final int HI_DP_THRESHOLD = 20; + private final double CALL_RATE_THRESHOLD = 0.90; + + @VisibleForTesting + public TransmittedSingleton(final Set trios) { + super((Set) null); + this.trios = Collections.unmodifiableSet(new LinkedHashSet<>(trios)); + } + + public TransmittedSingleton(final GATKPath pedigreeFile){ + super(pedigreeFile); + } + + public TransmittedSingleton(){ + super((Set) null); + } + + private Set initializeAndGetTrios() { + if (trios == null) { + trios = getTrios(); + } + final long numOfTrios = trios.stream().map(trio -> trio.getMother().getFamilyID()).count(); + if (numOfTrios == 0) { + oneShotLogger.warn("Submitted pedigree has no trios. TransmittedSingleton annotation will not be calculated."); + } else if (numOfTrios != trios.stream().map(trio -> trio.getMother().getFamilyID()).distinct().count()) { + oneShotLogger.warn("Submitted pedigree has non-trio families. TransmittedSingleton annotation is only valid for trios. Non-trio families (such as quads) will be ignored."); + } + return trios; + } + + @Override + void validateArguments(Collection founderIds, GATKPath pedigreeFile) { + String warningString = validateArgumentsWhenPedigreeRequired(founderIds, pedigreeFile); + if (warningString != null) { + warning.warn(warningString); + } + } + + @Override + public List getKeyNames() { + return Arrays.asList(GATKVCFConstants.TRANSMITTED_SINGLETON, GATKVCFConstants.NON_TRANSMITTED_SINGLETON); + } + @Override + public Map annotate(final ReferenceContext ref, + final VariantContext vc, + final AlleleLikelihoods likelihoods) { + Utils.nonNull(vc); + Set trioSet = initializeAndGetTrios(); + if (!vc.isBiallelic() || trioSet.isEmpty()) { + return Collections.emptyMap(); + } + long highQualCalls = vc.getGenotypes().stream().filter(gt -> gt.getGQ() > HI_GQ_THRESHOLD).count(); + if ((double) highQualCalls / vc.getNSamples() < CALL_RATE_THRESHOLD) { + return Collections.emptyMap(); + } + final List transmittedSingletonParent = new ArrayList<>(); + final List nonTransmittedSingletonParent = new ArrayList<>(); + for (final Trio trio : trioSet) { + if (vc.isBiallelic() && + contextHasTrioGQs(vc, trio)) { + + final boolean oneParentHasAllele = (vc.getGenotype(trio.getMaternalID()).isHet() && vc.getGenotype(trio.getPaternalID()).isHomRef()) || (vc.getGenotype(trio.getMaternalID()).isHomRef() && vc.getGenotype(trio.getPaternalID()).isHet()); + final String matchingParentId = vc.getGenotype(trio.getMaternalID()).isHet() ? trio.getMaternalID() : trio.getPaternalID(); + + final boolean momIsHighGQ = vc.getGenotype(trio.getMaternalID()).getGQ() >= HI_GQ_THRESHOLD; + final boolean dadIsHighGQ = vc.getGenotype(trio.getPaternalID()).getGQ() >= HI_GQ_THRESHOLD; + + final boolean childIsHighGQHet = vc.getGenotype(trio.getChildID()).isHet() && vc.getGenotype(trio.getChildID()).getGQ() >= HI_GQ_THRESHOLD; + final boolean childIsHighGQHomRef = vc.getGenotype(trio.getChildID()).isHomRef() && vc.getGenotype(trio.getChildID()).getGQ() >= HI_GQ_THRESHOLD; + + final boolean childIsHighDepth = vc.getGenotype(trio.getChildID()).getDP() >= HI_DP_THRESHOLD; + final boolean momIsHighDepth = vc.getGenotype(trio.getChildID()).getDP() >= HI_DP_THRESHOLD; + final boolean dadIsHighDepth = vc.getGenotype(trio.getChildID()).getDP() >= HI_DP_THRESHOLD; + + //TODO: This only works for trios (not quads or other more complicated family structures that would make the AC>2) + if (childIsHighDepth && momIsHighDepth && dadIsHighDepth && + vc.getAttributeAsInt(VCFConstants.ALLELE_COUNT_KEY, 0) == 2 && + childIsHighGQHet && oneParentHasAllele && momIsHighGQ && dadIsHighGQ) { + transmittedSingletonParent.add(matchingParentId); + } + //TODO: This only works for trios (not quads or other more complicated family structures that would make the AC>1) + if (childIsHighDepth && momIsHighDepth && dadIsHighDepth && + vc.getAttributeAsInt(VCFConstants.ALLELE_COUNT_KEY, 0) == 1 && + childIsHighGQHomRef && momIsHighGQ && dadIsHighGQ) { + nonTransmittedSingletonParent.add(matchingParentId); + } + } + } + + final Map attributeMap = new LinkedHashMap<>(1); + if (!transmittedSingletonParent.isEmpty()) { + attributeMap.put(GATKVCFConstants.TRANSMITTED_SINGLETON, transmittedSingletonParent); + } + if (!nonTransmittedSingletonParent.isEmpty()) { + attributeMap.put(GATKVCFConstants.NON_TRANSMITTED_SINGLETON, nonTransmittedSingletonParent); + } + return attributeMap; + } +} diff --git a/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFConstants.java b/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFConstants.java index ab6ee049e30..5870d1ad745 100644 --- a/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFConstants.java +++ b/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFConstants.java @@ -53,6 +53,8 @@ public final class GATKVCFConstants { public static final String GQ_STDEV_KEY = "GQ_STDDEV"; public static final String HAPLOTYPE_SCORE_KEY = "HaplotypeScore"; public static final String HI_CONF_DENOVO_KEY = "hiConfDeNovo"; + public static final String TRANSMITTED_SINGLETON = "transmittedSingleton"; + public static final String NON_TRANSMITTED_SINGLETON = "nonTransmittedSingleton"; public static final String INTERVAL_GC_CONTENT_KEY = "IGC"; public static final String INBREEDING_COEFFICIENT_KEY = "InbreedingCoeff"; public static final String AS_INBREEDING_COEFFICIENT_KEY = "AS_InbreedingCoeff"; diff --git a/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFHeaderLines.java b/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFHeaderLines.java index 0aa128eb777..21d48aeb104 100644 --- a/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFHeaderLines.java +++ b/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFHeaderLines.java @@ -180,6 +180,8 @@ public static VCFFormatHeaderLine getEquivalentFormatHeaderLine(final String inf addInfoLine(new VCFInfoHeaderLine(AS_VQS_SENS_KEY, VCFHeaderLineCount.A, VCFHeaderLineType.String, "For each alt allele, the calibration sensitivity threshold of being a true variant versus being false under the trained gaussian mixture model")); addInfoLine(new VCFInfoHeaderLine(AS_YNG_STATUS_KEY, VCFHeaderLineCount.A, VCFHeaderLineType.String, "For each alt allele, the yay/nay/grey status (yay are known good alleles, nay are known false positives, grey are unknown)")); + addInfoLine(new VCFInfoHeaderLine(TRANSMITTED_SINGLETON, 1, VCFHeaderLineType.String, "Possible transmitted singleton (site with AC=2 from parent and child). Parent ID is listed.")); + addInfoLine(new VCFInfoHeaderLine(NON_TRANSMITTED_SINGLETON, 1, VCFHeaderLineType.String, "Possible non transmitted singleton (site with AC=1 in just one parent). Parent ID is listed.")); addInfoLine(new VCFInfoHeaderLine(HI_CONF_DENOVO_KEY, 1, VCFHeaderLineType.String, "High confidence possible de novo mutation (GQ >= 20 for all trio members)=[comma-delimited list of child samples]")); addInfoLine(new VCFInfoHeaderLine(LO_CONF_DENOVO_KEY, 1, VCFHeaderLineType.String, "Low confidence possible de novo mutation (GQ >= 10 for child, GQ > 0 for parents)=[comma-delimited list of child samples]")); addInfoLine(new VCFInfoHeaderLine(QUAL_BY_DEPTH_KEY, 1, VCFHeaderLineType.Float, "Variant Confidence/Quality by Depth")); diff --git a/src/test/java/org/broadinstitute/hellbender/cmdline/GATKPlugin/GATKAnnotationPluginDescriptorUnitTest.java b/src/test/java/org/broadinstitute/hellbender/cmdline/GATKPlugin/GATKAnnotationPluginDescriptorUnitTest.java index 8920a445117..015e8b83c11 100644 --- a/src/test/java/org/broadinstitute/hellbender/cmdline/GATKPlugin/GATKAnnotationPluginDescriptorUnitTest.java +++ b/src/test/java/org/broadinstitute/hellbender/cmdline/GATKPlugin/GATKAnnotationPluginDescriptorUnitTest.java @@ -693,4 +693,17 @@ public void testConfigFileControlsAnnotationPackages() throws IOException { runToolInNewJVM("TestAnnotationsTool", args); Assert.assertEquals(Files.readAllLines(output.toPath()), Collections.singletonList(TestAnnotation.class.getSimpleName())); } + + @Test + public void testPossibleDeNovoAnnotationArguments() { + CommandLineParser clp = new CommandLineArgumentParser( + new Object(), + Collections.singletonList(new GATKAnnotationPluginDescriptor( Collections.singletonList(new PossibleDeNovo()), null)), + Collections.emptySet()); + String[] args = {"--" + PossibleDeNovo.DENOVO_DEPTH_THRESHOLD,"20","--" + PossibleDeNovo.DENOVO_PARENT_GQ_THRESHOLD,"40"}; + clp.parseArguments(nullMessageStream, args); + PossibleDeNovo annot = (PossibleDeNovo) instantiateAnnotations(clp).get(0); + Assert.assertEquals(annot.parentGQThreshold, 40); + Assert.assertEquals(annot.depthThreshold, 20); + } } \ No newline at end of file diff --git a/src/test/java/org/broadinstitute/hellbender/tools/walkers/annotator/PossibleDeNovoUnitTest.java b/src/test/java/org/broadinstitute/hellbender/tools/walkers/annotator/PossibleDeNovoUnitTest.java index 0cfb9b24b45..44e926fe864 100644 --- a/src/test/java/org/broadinstitute/hellbender/tools/walkers/annotator/PossibleDeNovoUnitTest.java +++ b/src/test/java/org/broadinstitute/hellbender/tools/walkers/annotator/PossibleDeNovoUnitTest.java @@ -16,6 +16,7 @@ public final class PossibleDeNovoUnitTest extends GATKBaseTest { private static final int GQ30 = 30; + private static final int GQ45 = 45; @DataProvider(name = "deNovo") public Object[][] deNovo() { @@ -30,105 +31,112 @@ public Object[][] deNovo() { final List het = Arrays.asList(refAllele, altAllele); final List notCall = Arrays.asList(noCallAllele, noCallAllele); - final Genotype g00 = new GenotypeBuilder("2", homRef).DP(10).PL(new int[]{0,100,100}).AD(new int[]{10, 0}).GQ(GQ30).make(); - final Genotype g01 = new GenotypeBuilder("1", het).DP(10).PL(new int[]{100, 0, 100}).AD(new int[]{5, 5}).GQ(GQ30).make(); - final Genotype g11 = new GenotypeBuilder("3", homVar).DP(10).PL(new int[]{100, 100, 0}).AD(new int[]{0, 10}).GQ(GQ30).make(); + final Genotype g00gq30 = new GenotypeBuilder("2", homRef).DP(10).PL(new int[]{0,100,100}).AD(new int[]{10, 0}).GQ(GQ30).make(); + final Genotype g01gq30 = new GenotypeBuilder("1", het).DP(10).PL(new int[]{100, 0, 100}).AD(new int[]{5, 5}).GQ(GQ30).make(); + final Genotype g11gq30 = new GenotypeBuilder("3", homVar).DP(10).PL(new int[]{100, 100, 0}).AD(new int[]{0, 10}).GQ(GQ30).make(); final Genotype gNo = new GenotypeBuilder("4", notCall).make(); - tests.add(new Object[]{g00, g00, g00, false}); - tests.add(new Object[]{g00, g00, g01, true}); - tests.add(new Object[]{g00, g00, g11, false}); - tests.add(new Object[]{g00, g00, gNo, false}); - - tests.add(new Object[]{g00, g01, g00, false}); - tests.add(new Object[]{g00, g01, g01, false}); - tests.add(new Object[]{g00, g01, g11, false}); - tests.add(new Object[]{g00, g01, gNo, false}); - - tests.add(new Object[]{g00, g11, g00, false}); - tests.add(new Object[]{g00, g11, g01, false}); - tests.add(new Object[]{g00, g11, g11, false}); - tests.add(new Object[]{g00, g11, gNo, false}); - - tests.add(new Object[]{g00, gNo, g00, false}); - tests.add(new Object[]{g00, gNo, g01, false}); - tests.add(new Object[]{g00, gNo, g11, false}); - tests.add(new Object[]{g00, gNo, gNo, false}); - - tests.add(new Object[]{g01, g00, g00, false}); - tests.add(new Object[]{g01, g00, g01, false}); - tests.add(new Object[]{g01, g00, g11, false}); - tests.add(new Object[]{g01, g00, gNo, false}); - - tests.add(new Object[]{g01, g01, g00, false}); - tests.add(new Object[]{g01, g01, g01, false}); - tests.add(new Object[]{g01, g01, g11, false}); - tests.add(new Object[]{g01, g01, gNo, false}); - - tests.add(new Object[]{g01, g11, g00, false}); - tests.add(new Object[]{g01, g11, g01, false}); - tests.add(new Object[]{g01, g11, g11, false}); - tests.add(new Object[]{g01, g11, gNo, false}); - - tests.add(new Object[]{g01, gNo, g00, false}); - tests.add(new Object[]{g01, gNo, g01, false}); - tests.add(new Object[]{g01, gNo, g11, false}); - tests.add(new Object[]{g01, gNo, gNo, false}); - - - tests.add(new Object[]{g11, g00, g00, false}); - tests.add(new Object[]{g11, g00, g01, false}); - tests.add(new Object[]{g11, g00, g11, false}); - tests.add(new Object[]{g11, g00, gNo, false}); - - tests.add(new Object[]{g11, g01, g00, false}); - tests.add(new Object[]{g11, g01, g01, false}); - tests.add(new Object[]{g11, g01, g11, false}); - tests.add(new Object[]{g11, g01, gNo, false}); - - tests.add(new Object[]{g11, g11, g00, false}); - tests.add(new Object[]{g11, g11, g01, false}); - tests.add(new Object[]{g11, g11, g11, false}); - tests.add(new Object[]{g11, g11, gNo, false}); - - tests.add(new Object[]{g11, gNo, g00, false}); - tests.add(new Object[]{g11, gNo, g01, false}); - tests.add(new Object[]{g11, gNo, g11, false}); - tests.add(new Object[]{g11, gNo, gNo, false}); - - - tests.add(new Object[]{gNo, g00, g00, false}); - tests.add(new Object[]{gNo, g00, g01, false}); - tests.add(new Object[]{gNo, g00, g11, false}); - tests.add(new Object[]{gNo, g00, gNo, false}); - - tests.add(new Object[]{gNo, g01, g00, false}); - tests.add(new Object[]{gNo, g01, g01, false}); - tests.add(new Object[]{gNo, g01, g11, false}); - tests.add(new Object[]{gNo, g01, gNo, false}); - - tests.add(new Object[]{gNo, g11, g00, false}); - tests.add(new Object[]{gNo, g11, g01, false}); - tests.add(new Object[]{gNo, g11, g11, false}); - tests.add(new Object[]{gNo, g11, gNo, false}); - - tests.add(new Object[]{gNo, gNo, g00, false}); - tests.add(new Object[]{gNo, gNo, g01, false}); - tests.add(new Object[]{gNo, gNo, g11, false}); - tests.add(new Object[]{gNo, gNo, gNo, false}); + tests.add(new Object[]{g00gq30, g00gq30, g00gq30, false, false}); + //test case where there is a possible DeNovo, but depth and parent GQs are lower than set thresholds (20 and 40 respectively) + tests.add(new Object[]{g00gq30, g00gq30, g01gq30, true, false}); + tests.add(new Object[]{g00gq30, g00gq30, g11gq30, false, false}); + tests.add(new Object[]{g00gq30, g00gq30, gNo, false, false}); + + tests.add(new Object[]{g00gq30, g01gq30, g00gq30, false, false}); + tests.add(new Object[]{g00gq30, g01gq30, g01gq30, false, false}); + tests.add(new Object[]{g00gq30, g01gq30, g11gq30, false, false}); + tests.add(new Object[]{g00gq30, g01gq30, gNo, false, false}); + + tests.add(new Object[]{g00gq30, g11gq30, g00gq30, false, false}); + tests.add(new Object[]{g00gq30, g11gq30, g01gq30, false, false}); + tests.add(new Object[]{g00gq30, g11gq30, g11gq30, false, false}); + tests.add(new Object[]{g00gq30, g11gq30, gNo, false, false}); + + tests.add(new Object[]{g00gq30, gNo, g00gq30, false, false}); + tests.add(new Object[]{g00gq30, gNo, g01gq30, false, false}); + tests.add(new Object[]{g00gq30, gNo, g11gq30, false, false}); + tests.add(new Object[]{g00gq30, gNo, gNo, false, false}); + + tests.add(new Object[]{g01gq30, g00gq30, g00gq30, false, false}); + tests.add(new Object[]{g01gq30, g00gq30, g01gq30, false, false}); + tests.add(new Object[]{g01gq30, g00gq30, g11gq30, false, false}); + tests.add(new Object[]{g01gq30, g00gq30, gNo, false, false}); + + tests.add(new Object[]{g01gq30, g01gq30, g00gq30, false, false}); + tests.add(new Object[]{g01gq30, g01gq30, g01gq30, false, false}); + tests.add(new Object[]{g01gq30, g01gq30, g11gq30, false, false}); + tests.add(new Object[]{g01gq30, g01gq30, gNo, false, false}); + + tests.add(new Object[]{g01gq30, g11gq30, g00gq30, false, false}); + tests.add(new Object[]{g01gq30, g11gq30, g01gq30, false, false}); + tests.add(new Object[]{g01gq30, g11gq30, g11gq30, false, false}); + tests.add(new Object[]{g01gq30, g11gq30, gNo, false, false}); + + tests.add(new Object[]{g01gq30, gNo, g00gq30, false, false}); + tests.add(new Object[]{g01gq30, gNo, g01gq30, false, false}); + tests.add(new Object[]{g01gq30, gNo, g11gq30, false, false}); + tests.add(new Object[]{g01gq30, gNo, gNo, false, false}); + + + tests.add(new Object[]{g11gq30, g00gq30, g00gq30, false, false}); + tests.add(new Object[]{g11gq30, g00gq30, g01gq30, false, false}); + tests.add(new Object[]{g11gq30, g00gq30, g11gq30, false, false}); + tests.add(new Object[]{g11gq30, g00gq30, gNo, false, false}); + + tests.add(new Object[]{g11gq30, g01gq30, g00gq30, false, false}); + tests.add(new Object[]{g11gq30, g01gq30, g01gq30, false, false}); + tests.add(new Object[]{g11gq30, g01gq30, g11gq30, false, false}); + tests.add(new Object[]{g11gq30, g01gq30, gNo, false, false}); + + tests.add(new Object[]{g11gq30, g11gq30, g00gq30, false, false}); + tests.add(new Object[]{g11gq30, g11gq30, g01gq30, false, false}); + tests.add(new Object[]{g11gq30, g11gq30, g11gq30, false, false}); + tests.add(new Object[]{g11gq30, g11gq30, gNo, false, false}); + + tests.add(new Object[]{g11gq30, gNo, g00gq30, false, false}); + tests.add(new Object[]{g11gq30, gNo, g01gq30, false, false}); + tests.add(new Object[]{g11gq30, gNo, g11gq30, false, false}); + tests.add(new Object[]{g11gq30, gNo, gNo, false, false}); + + + tests.add(new Object[]{gNo, g00gq30, g00gq30, false, false}); + tests.add(new Object[]{gNo, g00gq30, g01gq30, false, false}); + tests.add(new Object[]{gNo, g00gq30, g11gq30, false, false}); + tests.add(new Object[]{gNo, g00gq30, gNo, false, false}); + + tests.add(new Object[]{gNo, g01gq30, g00gq30, false, false}); + tests.add(new Object[]{gNo, g01gq30, g01gq30, false, false}); + tests.add(new Object[]{gNo, g01gq30, g11gq30, false, false}); + tests.add(new Object[]{gNo, g01gq30, gNo, false, false}); + + tests.add(new Object[]{gNo, g11gq30, g00gq30, false, false}); + tests.add(new Object[]{gNo, g11gq30, g01gq30, false, false}); + tests.add(new Object[]{gNo, g11gq30, g11gq30, false, false}); + tests.add(new Object[]{gNo, g11gq30, gNo, false, false}); + + tests.add(new Object[]{gNo, gNo, g00gq30, false, false}); + tests.add(new Object[]{gNo, gNo, g01gq30, false, false}); + tests.add(new Object[]{gNo, gNo, g11gq30, false, false}); + tests.add(new Object[]{gNo, gNo, gNo, false, false}); final Genotype g00NoPl = new GenotypeBuilder("2", homRef).DP(10).GQ(GQ30).make(); final Genotype g01NoPl = new GenotypeBuilder("1", het).DP(10).GQ(GQ30).make(); final Genotype g11NoPl = new GenotypeBuilder("3", homVar).DP(10).GQ(GQ30).make(); - tests.add(new Object[]{g00NoPl, g01NoPl, g11NoPl, false}); - tests.add(new Object[]{g00NoPl, g00NoPl, g01NoPl, true}); + tests.add(new Object[]{g00NoPl, g01NoPl, g11NoPl, false, false}); + tests.add(new Object[]{g00NoPl, g00NoPl, g01NoPl, true, false}); + + final Genotype g00gq45 = new GenotypeBuilder("2", homRef).DP(50).PL(new int[]{0,100,100}).AD(new int[]{10, 0}).GQ(GQ45).make(); + final Genotype g01gq45 = new GenotypeBuilder("1", het).DP(50).PL(new int[]{100, 0, 100}).AD(new int[]{5, 5}).GQ(GQ45).make(); + + //test case when there is a possible deNovo (het child with hom ref parents) and GQ is high enough to pass parent GQ threshold of 40 and depth threshold of 20 + tests.add(new Object[]{g00gq45, g00gq45, g01gq45, true, true}); return tests.toArray(new Object[][]{}); } @Test(dataProvider = "deNovo") - public void testUsingVC(final Genotype gMom, final Genotype gDad, final Genotype gChild, final boolean expectedDeNovo) { + public void testUsingVC(final Genotype gMom, final Genotype gDad, final Genotype gChild, final boolean expectedDeNovo, final boolean expectedDeNovoWithExtraThresholds) { final Sample sMom = new Sample("mom", "fam", null, null, Sex.FEMALE); final Sample sDad = new Sample("dad", sMom.getFamilyID(), null, null, Sex.MALE); @@ -157,11 +165,18 @@ public void testUsingVC(final Genotype gMom, final Genotype gDad, final Genotype final Trio trio = new Trio(sMom, sDad, sChild); final InfoFieldAnnotation ann = new PossibleDeNovo(Collections.singleton(trio), 0); + final InfoFieldAnnotation hiQualAnn = new PossibleDeNovo(Collections.singleton(trio), 0, 40, 20); final Map result = ann.annotate(null, vc, null); + final Map hiQualResult = hiQualAnn.annotate(null, vc, null); if (expectedDeNovo) { - Assert.assertEquals(result.get(GATKVCFConstants.HI_CONF_DENOVO_KEY), Arrays.asList(sChild.getID())); + Assert.assertEquals(result.get(GATKVCFConstants.HI_CONF_DENOVO_KEY), Collections.singletonList(sChild.getID())); Assert.assertFalse(result.containsKey(GATKVCFConstants.LO_CONF_DENOVO_KEY)); - } else { + } + if (expectedDeNovoWithExtraThresholds) { + Assert.assertEquals(hiQualResult.get(GATKVCFConstants.HI_CONF_DENOVO_KEY), Collections.singletonList(sChild.getID())); + Assert.assertFalse(result.containsKey(GATKVCFConstants.LO_CONF_DENOVO_KEY)); + } + if (!expectedDeNovo && !expectedDeNovoWithExtraThresholds) { Assert.assertTrue(result.isEmpty()); } diff --git a/src/test/java/org/broadinstitute/hellbender/tools/walkers/annotator/TransmittedSingletonUnitTest.java b/src/test/java/org/broadinstitute/hellbender/tools/walkers/annotator/TransmittedSingletonUnitTest.java new file mode 100644 index 00000000000..dd03ed0bf87 --- /dev/null +++ b/src/test/java/org/broadinstitute/hellbender/tools/walkers/annotator/TransmittedSingletonUnitTest.java @@ -0,0 +1,234 @@ +package org.broadinstitute.hellbender.tools.walkers.annotator; + +import htsjdk.variant.variantcontext.*; +import htsjdk.variant.vcf.VCFConstants; +import org.broadinstitute.hellbender.utils.samples.Sample; +import org.broadinstitute.hellbender.utils.samples.Sex; +import org.broadinstitute.hellbender.utils.samples.Trio; +import org.broadinstitute.hellbender.utils.variant.GATKVCFConstants; +import org.broadinstitute.hellbender.utils.variant.GATKVCFHeaderLines; +import org.testng.Assert; +import org.testng.annotations.DataProvider; +import org.testng.annotations.Test; + +import java.util.*; + +public class TransmittedSingletonUnitTest { + + private static final int GQ30 = 30; + public static final String MOM = "mom"; + public static final String DAD = "dad"; + public static final String CHILD = "child"; + + @DataProvider(name = "transmittedSingletons") + public Object[][] transmittedSingletons() { + final List tests = new ArrayList<>(); + + final Allele refAllele = Allele.create("A", true); + final Allele altAllele = Allele.create("T"); + final Allele noCallAllele = Allele.NO_CALL; + + final List homRef = Arrays.asList(refAllele, refAllele); + final List homVar = Arrays.asList(altAllele, altAllele); + final List het = Arrays.asList(refAllele, altAllele); + final List notCall = Arrays.asList(noCallAllele, noCallAllele); + + final Genotype g00gq30 = new GenotypeBuilder("2", homRef).DP(30).PL(new int[]{0,100,100}).AD(new int[]{10, 0}).GQ(GQ30).make(); + final Genotype g01gq30 = new GenotypeBuilder("1", het).DP(30).PL(new int[]{100, 0, 100}).AD(new int[]{5, 5}).GQ(GQ30).make(); + final Genotype g11gq30 = new GenotypeBuilder("3", homVar).DP(30).PL(new int[]{100, 100, 0}).AD(new int[]{0, 10}).GQ(GQ30).make(); + final Genotype gNo = new GenotypeBuilder("4", notCall).make(); + + tests.add(new Object[]{g00gq30, g00gq30, g00gq30, false, false, null}); + tests.add(new Object[]{g00gq30, g00gq30, g01gq30, false, false, null}); + tests.add(new Object[]{g00gq30, g00gq30, g11gq30, false, false, null}); + tests.add(new Object[]{g00gq30, g00gq30, gNo, false, false, null}); + + tests.add(new Object[]{g00gq30, g01gq30, g00gq30, false, true, DAD}); + tests.add(new Object[]{g00gq30, g01gq30, g01gq30, true, false, DAD}); + tests.add(new Object[]{g00gq30, g01gq30, g11gq30, false, false, null}); + tests.add(new Object[]{g00gq30, g01gq30, gNo, false, false, null}); + + tests.add(new Object[]{g00gq30, g11gq30, g00gq30, false, false, null}); + tests.add(new Object[]{g00gq30, g11gq30, g01gq30, false, false, null}); + tests.add(new Object[]{g00gq30, g11gq30, g11gq30, false, false, null}); + tests.add(new Object[]{g00gq30, g11gq30, gNo, false, false, null}); + + tests.add(new Object[]{g00gq30, gNo, g00gq30, false, false, null}); + tests.add(new Object[]{g00gq30, gNo, g01gq30, false, false, null}); + tests.add(new Object[]{g00gq30, gNo, g11gq30, false, false, null}); + tests.add(new Object[]{g00gq30, gNo, gNo, false, false, null}); + + tests.add(new Object[]{g01gq30, g00gq30, g00gq30, false, true, MOM}); + tests.add(new Object[]{g01gq30, g00gq30, g01gq30, true, false, MOM}); + tests.add(new Object[]{g01gq30, g00gq30, g11gq30, false, false, null}); + tests.add(new Object[]{g01gq30, g00gq30, gNo, false, false, null}); + + tests.add(new Object[]{g01gq30, g01gq30, g00gq30, false, false, null}); + tests.add(new Object[]{g01gq30, g01gq30, g01gq30, false, false, null}); + tests.add(new Object[]{g01gq30, g01gq30, g11gq30, false, false, null}); + tests.add(new Object[]{g01gq30, g01gq30, gNo, false, false, null}); + + tests.add(new Object[]{g01gq30, g11gq30, g00gq30, false, false, null}); + tests.add(new Object[]{g01gq30, g11gq30, g01gq30, false, false, null}); + tests.add(new Object[]{g01gq30, g11gq30, g11gq30, false, false, null}); + tests.add(new Object[]{g01gq30, g11gq30, gNo, false, false, null}); + + tests.add(new Object[]{g01gq30, gNo, g00gq30, false, false, null}); + tests.add(new Object[]{g01gq30, gNo, g01gq30, false, false, null}); + tests.add(new Object[]{g01gq30, gNo, g11gq30, false, false, null}); + tests.add(new Object[]{g01gq30, gNo, gNo, false, false, null}); + + + tests.add(new Object[]{g11gq30, g00gq30, g00gq30, false, false, null}); + tests.add(new Object[]{g11gq30, g00gq30, g01gq30, false, false, null}); + tests.add(new Object[]{g11gq30, g00gq30, g11gq30, false, false, null}); + tests.add(new Object[]{g11gq30, g00gq30, gNo, false, false, null}); + + tests.add(new Object[]{g11gq30, g01gq30, g00gq30, false, false, null}); + tests.add(new Object[]{g11gq30, g01gq30, g01gq30, false, false, null}); + tests.add(new Object[]{g11gq30, g01gq30, g11gq30, false, false, null}); + tests.add(new Object[]{g11gq30, g01gq30, gNo, false, false, null}); + + tests.add(new Object[]{g11gq30, g11gq30, g00gq30, false, false, null}); + tests.add(new Object[]{g11gq30, g11gq30, g01gq30, false, false, null}); + tests.add(new Object[]{g11gq30, g11gq30, g11gq30, false, false, null}); + tests.add(new Object[]{g11gq30, g11gq30, gNo, false, false, null}); + + tests.add(new Object[]{g11gq30, gNo, g00gq30, false, false, null}); + tests.add(new Object[]{g11gq30, gNo, g01gq30, false, false, null}); + tests.add(new Object[]{g11gq30, gNo, g11gq30, false, false, null}); + tests.add(new Object[]{g11gq30, gNo, gNo, false, false, null}); + + + tests.add(new Object[]{gNo, g00gq30, g00gq30, false, false, null}); + tests.add(new Object[]{gNo, g00gq30, g01gq30, false, false, null}); + tests.add(new Object[]{gNo, g00gq30, g11gq30, false, false, null}); + tests.add(new Object[]{gNo, g00gq30, gNo, false, false, null}); + + tests.add(new Object[]{gNo, g01gq30, g00gq30, false, false, null}); + tests.add(new Object[]{gNo, g01gq30, g01gq30, false, false, null}); + tests.add(new Object[]{gNo, g01gq30, g11gq30, false, false, null}); + tests.add(new Object[]{gNo, g01gq30, gNo, false, false, null}); + + tests.add(new Object[]{gNo, g11gq30, g00gq30, false, false, null}); + tests.add(new Object[]{gNo, g11gq30, g01gq30, false, false, null}); + tests.add(new Object[]{gNo, g11gq30, g11gq30, false, false, null}); + tests.add(new Object[]{gNo, g11gq30, gNo, false, false, null}); + + tests.add(new Object[]{gNo, gNo, g00gq30, false, false, null}); + tests.add(new Object[]{gNo, gNo, g01gq30, false, false, null}); + tests.add(new Object[]{gNo, gNo, g11gq30, false, false, null}); + tests.add(new Object[]{gNo, gNo, gNo, false, false, null}); + + final Genotype g00NoPl = new GenotypeBuilder("2", homRef).DP(30).GQ(GQ30).make(); + final Genotype g01NoPl = new GenotypeBuilder("1", het).DP(30).GQ(GQ30).make(); + final Genotype g11NoPl = new GenotypeBuilder("3", homVar).DP(30).GQ(GQ30).make(); + + tests.add(new Object[]{g00NoPl, g01NoPl, g11NoPl, false, false, null}); + tests.add(new Object[]{g00NoPl, g00NoPl, g01NoPl, false, false, null}); + tests.add(new Object[]{g00NoPl, g01NoPl, g01NoPl, true, false, DAD}); + tests.add(new Object[]{g00NoPl, g01NoPl, g00NoPl, false, true, DAD}); + + return tests.toArray(new Object[][]{}); + } + + @Test(dataProvider = "transmittedSingletons") + public void testUsingVC(final Genotype gMom, final Genotype gDad, final Genotype gChild, final boolean expectedTransmittedSingleton, + final boolean expectedNonTransmittedSingleton, final String matchingParentSample) { + + final Sample sMom = new Sample(MOM, "fam", null, null, Sex.FEMALE); + final Sample sDad = new Sample(DAD, sMom.getFamilyID(), null, null, Sex.MALE); + final Sample sChild = new Sample(CHILD, sMom.getFamilyID(), sDad.getID(), sMom.getID(), Sex.MALE); + + final Allele refAllele = Allele.create("A", true); + final Allele altAllele = Allele.create("T"); + final Allele noCallAllele = Allele.NO_CALL; + + final ArrayList genotypes = new ArrayList<>(Arrays.asList(gMom, gDad, gChild)); + final Map sampleNameToOffset= new LinkedHashMap<>(3); + sampleNameToOffset.put(MOM, 0); + sampleNameToOffset.put(DAD, 1); + sampleNameToOffset.put(CHILD, 2); + final List sampleNamesInOrder=Arrays.asList(MOM, DAD, CHILD); + final GenotypesContext context = GenotypesContext.create(genotypes, sampleNameToOffset, sampleNamesInOrder); + + final Collection alleles= new LinkedHashSet<>(2); + alleles.addAll(gMom.getAlleles()); + alleles.addAll(gDad.getAlleles()); + alleles.addAll(gChild.getAlleles()); + alleles.add(refAllele); + alleles.add(altAllele); + alleles.remove(noCallAllele); + final int ac = (int) context.stream().mapToLong(g-> g.getAlleles().stream().filter(a->a.isCalled()&&a.isNonReference()).count()).sum(); + final VariantContext vc = new VariantContextBuilder("test", "20", 10, 10, alleles).genotypes(context).make(); + vc.getCommonInfo().putAttribute(VCFConstants.ALLELE_COUNT_KEY, ac); + + final Trio trio = new Trio(sMom, sDad, sChild); + final InfoFieldAnnotation ann = new TransmittedSingleton(Collections.singleton(trio)); + final Map result = ann.annotate(null, vc, null); + if (expectedTransmittedSingleton) { + Assert.assertEquals(result.get(GATKVCFConstants.TRANSMITTED_SINGLETON), Collections.singletonList(matchingParentSample)); + } else if (expectedNonTransmittedSingleton) { + Assert.assertEquals(result.get(GATKVCFConstants.NON_TRANSMITTED_SINGLETON), Collections.singletonList(matchingParentSample)); + } else { + Assert.assertTrue(result.isEmpty()); + } + + Assert.assertEquals(ann.getDescriptions(), Arrays.asList(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.TRANSMITTED_SINGLETON), + GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.NON_TRANSMITTED_SINGLETON))); + Assert.assertEquals(ann.getKeyNames(), Arrays.asList(GATKVCFConstants.TRANSMITTED_SINGLETON, GATKVCFConstants.NON_TRANSMITTED_SINGLETON)); + + final Map resultNoTrio = new TransmittedSingleton(Collections.emptySet()).annotate(null, vc, null); + Assert.assertTrue(resultNoTrio.isEmpty()); + } + + // tests quad family which won't work for transmitted singletons + @Test + public void testNonTrioFamilies() { + final Sample sMom = new Sample(MOM, "fam", null, null, Sex.FEMALE); + final Sample sDad = new Sample(DAD, sMom.getFamilyID(), null, null, Sex.MALE); + final Sample sSib1 = new Sample(CHILD, sMom.getFamilyID(), sDad.getID(), sMom.getID(), Sex.MALE); + final Sample sSib2 = new Sample(CHILD + "2", sMom.getFamilyID(), sDad.getID(), sMom.getID(), Sex.FEMALE); + + final Allele refAllele = Allele.create("A", true); + final Allele altAllele = Allele.create("T"); + + final List homRef = Arrays.asList(refAllele, refAllele); + final List het = Arrays.asList(refAllele, altAllele); + + final Genotype gMom = new GenotypeBuilder("2", homRef).DP(30).PL(new int[]{0,100,100}).AD(new int[]{10, 0}).GQ(GQ30).make(); + final Genotype gDad = new GenotypeBuilder("1", het).DP(30).PL(new int[]{100, 0, 100}).AD(new int[]{5, 5}).GQ(GQ30).make(); + final Genotype gSib1 = new GenotypeBuilder("3", het).DP(30).PL(new int[]{100, 100, 0}).AD(new int[]{0, 10}).GQ(GQ30).make(); + final Genotype gSib2 = new GenotypeBuilder("4", het).DP(30).PL(new int[]{100, 0, 100}).AD(new int[]{5, 5}).GQ(GQ30).make(); + + final ArrayList genotypes = new ArrayList<>(Arrays.asList(gMom, gDad, gSib1, gSib2)); + final Map sampleNameToOffset= new LinkedHashMap<>(3); + sampleNameToOffset.put(MOM, 0); + sampleNameToOffset.put(DAD, 1); + sampleNameToOffset.put(CHILD, 2); + sampleNameToOffset.put(CHILD + "2", 3); + final List sampleNamesInOrder=Arrays.asList(MOM, DAD, CHILD, CHILD+"2"); + final GenotypesContext context = GenotypesContext.create(genotypes, sampleNameToOffset, sampleNamesInOrder); + + final Collection alleles= new LinkedHashSet<>(2); + alleles.addAll(gMom.getAlleles()); + alleles.addAll(gDad.getAlleles()); + alleles.addAll(gSib1.getAlleles()); + alleles.addAll(gSib2.getAlleles()); + alleles.add(refAllele); + alleles.add(altAllele); + final int ac = (int) context.stream().mapToLong(g-> g.getAlleles().stream().filter(a->a.isCalled()&&a.isNonReference()).count()).sum(); + final VariantContext vc = new VariantContextBuilder("test", "20", 10, 10, alleles).genotypes(context).make(); + vc.getCommonInfo().putAttribute(VCFConstants.ALLELE_COUNT_KEY, ac); + + final Trio trio1 = new Trio(sMom, sDad, sSib1); + final Trio trio2 = new Trio(sMom, sDad, sSib2); + + final InfoFieldAnnotation ann = new TransmittedSingleton(new HashSet<>(Arrays.asList(trio1, trio2))); + final Map result = ann.annotate(null, vc, null); + + //quads are not implemented yet so we skip over them + Assert.assertEquals(result.size(), 0); + } + +}