Fixed a bug in AlleleFiltering that ignored more than a single sample (…

…#8841)
broadinstitute · Jun 13, 2024 · ab98a5d · ab98a5d
1 parent 2a420e4
commit ab98a5d
Show file tree

Hide file tree

Showing 2 changed files with 82 additions and 12 deletions.
diff --git a/...n/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AlleleFilteringHC.java b/...n/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AlleleFilteringHC.java
@@ -11,7 +11,10 @@
 import org.broadinstitute.hellbender.utils.read.GATKRead;
 
 import java.io.OutputStreamWriter;
+import java.util.ArrayList;
 import java.util.Arrays;
+import java.util.Collections;
+import java.util.List;
 
 /**
  * Filtering haplotypes that contribute weak alleles to the genotyping. This is a version that determines if allele is weak using
@@ -57,14 +60,16 @@ int getAlleleLikelihoodVsInverse(final AlleleLikelihoods<GATKRead, Allele> allel
 
         final GenotypingLikelihoods<Allele> genotypingLikelihoods = genotypesModel.calculateLikelihoods(alleleList,
                 genotypingData, null, 0, null);
-        AFCalculationResult af = afCalc.fastCalculateDiploidBasedOnGLs(genotypingLikelihoods, genotypingEngine.getPloidyModel().totalPloidy());
-        final double log10Confidence = af.log10ProbOnlyRefAlleleExists();
-        final double phredScaledConfidence = (10.0 * log10Confidence) + 0.0;
 
-        final int[] asPL = genotypingLikelihoods.sampleLikelihoods(0).getAsPLs();
-
-        logger.debug(() -> String.format("GAL:: %s: %d %d %d", allele.toString(), asPL[0], asPL[1], asPL[2]));
-        return Math.min(asPL[1]-asPL[0], asPL[2]-asPL[0]);
+        List<Integer> perSamplePLs = new ArrayList<>();
+        for (int i = 0; i < genotypingLikelihoods.numberOfSamples(); i++) {
+            final int[] pls = genotypingLikelihoods.sampleLikelihoods(i).getAsPLs();
+            perSamplePLs.add(Math.min(pls[1] - pls[0], pls[2] - pls[0]));
+            final int finalI = i;
+            logger.debug(() -> String.format("GAL (%s):: %s: %d %d %d",
+                    genotypingLikelihoods.getSample(finalI), allele.toString(), pls[0], pls[1], pls[2]));
+        }
+        return Collections.min(perSamplePLs);
     }
 
 }
diff --git a/.../org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AlleleFilteringUnitTest.java b/.../org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AlleleFilteringUnitTest.java
@@ -1,12 +1,7 @@
 package org.broadinstitute.hellbender.tools.walkers.haplotypecaller;
 
-import com.google.common.collect.ImmutableList;
 import htsjdk.samtools.TextCigarCodec;
 import htsjdk.variant.variantcontext.Allele;
-import org.broadinstitute.hellbender.tools.walkers.haplotypecaller.AlleleFiltering;
-import org.broadinstitute.hellbender.tools.walkers.haplotypecaller.AlleleFilteringHC;
-import org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCallerArgumentCollection;
-import org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCallerGenotypingEngine;
 import org.broadinstitute.hellbender.utils.SimpleInterval;
 import org.broadinstitute.hellbender.utils.genotyper.*;
 import org.broadinstitute.hellbender.utils.haplotype.Event;
@@ -67,6 +62,76 @@ public void testNoNeedToFilter(){
         AlleleLikelihoods<GATKRead, Haplotype> filtered_lks = alleleFiltering.filterAlleles(lks, 0, new HashSet<>());
         Assert.assertEquals(filtered_lks.alleles(), lks.alleles());
     }
+    @Test
+    public void testNoNeedToFilterTwoSamples(){
+        List<Haplotype> haplotypeList = new ArrayList<>();
+        final byte[] fullReferenceWithPadding = "CATGCATG".getBytes();
+        Haplotype haplotype = new Haplotype(fullReferenceWithPadding, true, 0, TextCigarCodec.decode("8M"));
+        haplotype.setGenomeLocation(new SimpleInterval("chr", 100, 108));
+        haplotype.setEventMap(EventMap.fromHaplotype(haplotype, fullReferenceWithPadding, 0));
+        haplotypeList.add(haplotype);
+
+        haplotype = new Haplotype("CAGGCATG".getBytes(), false, 0, TextCigarCodec.decode("8M"));
+        haplotype.setGenomeLocation(new SimpleInterval("chr", 100, 108));
+        haplotypeList.add(haplotype);
+        haplotype.setEventMap(EventMap.fromHaplotype(haplotype, fullReferenceWithPadding, 0));
+
+        haplotype = new Haplotype("CATTCATG".getBytes(), false, 0, TextCigarCodec.decode("8M"));
+        haplotype.setGenomeLocation(new SimpleInterval("chr", 100, 108));
+        haplotypeList.add(haplotype);
+        haplotype.setEventMap(EventMap.fromHaplotype(haplotype, fullReferenceWithPadding, 0));
+
+        AlleleList<Haplotype> haplotypes = new IndexedAlleleList<>(haplotypeList);
+        SampleList samples = new IndexedSampleList(Arrays.asList(new String[]{"sm1", "sm2"}));
+
+        List<GATKRead> readList = new ArrayList<>(30);
+        Map<String, List<GATKRead>>ebs = new HashMap<>();
+        ebs.put("sm1", readList);
+
+        for (int i = 0 ; i < 30; i++) {
+            readList.add(ArtificialReadUtils.createArtificialRead("20M"));
+        }
+
+        List<GATKRead> readList2 = new ArrayList<>(30);
+        ebs.put("sm2", readList2);
+        for (int i = 0 ; i < 30; i++) {
+            readList2.add(ArtificialReadUtils.createArtificialRead("20M"));
+        }
+
+        double[][] values = {{0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3},
+                {3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0},
+                {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
+        };
+        AlleleLikelihoods<GATKRead, Haplotype> lks = new AlleleLikelihoods<>(samples, haplotypes, ebs);
+        LikelihoodMatrix<GATKRead, Haplotype> lkm = lks.sampleMatrix(0);
+        for (int i = 0; i < lks.numberOfAlleles(); i++){
+            for (int j = 0 ; j < lkm.evidenceCount(); j++) {
+                lkm.set(i,j,values[i][j]);
+            }
+        }
+
+
+        double[][] values2 = {{0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3},
+                {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+                {3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0}
+        };
+        LikelihoodMatrix<GATKRead, Haplotype> lkm2 = lks.sampleMatrix(1);
+        for (int i = 0; i < lks.numberOfAlleles(); i++){
+            for (int j = 0 ; j < lkm2.evidenceCount(); j++) {
+                lkm2.set(i,j,values2[i][j]);
+            }
+        }
+
+        HaplotypeCallerArgumentCollection hcArgs = new HaplotypeCallerArgumentCollection();
+        HaplotypeCallerGenotypingEngine genotypingEngine = new HaplotypeCallerGenotypingEngine(hcArgs, samples, ! hcArgs.doNotRunPhysicalPhasing, false);
+
+
+        AlleleFiltering alleleFiltering = new AlleleFilteringHC(hcArgs, null, genotypingEngine);
+        AlleleLikelihoods<GATKRead, Haplotype> filtered_lks = alleleFiltering.filterAlleles(lks, 0, new HashSet<>());
+        Assert.assertEquals(filtered_lks.alleles(), lks.alleles());
+    }
+
+
 
     @Test
     public void testFilterCloseMis(){