Purple: improved somatic fit calc purity value matching

hartwigmedical · Jan 15, 2025 · 3b8b2d0 · 3b8b2d0
1 parent 9cd9ba3
commit 3b8b2d0
Show file tree

Hide file tree

Showing 5 changed files with 54 additions and 18 deletions.
diff --git a/purple/src/main/java/com/hartwig/hmftools/purple/FittingConfig.java b/purple/src/main/java/com/hartwig/hmftools/purple/FittingConfig.java
@@ -119,6 +119,11 @@ public FittingConfig(final ConfigBuilder configBuilder, boolean targetedMode)
         GcRatioExponent = getConfigDecimal(configBuilder, GC_RATIO_EXPONENT, targetedMode ? TARGETED_GC_RATIO_EXPONENT_DEFAULT : 0);
     }
 
+    public boolean hasValidValues()
+    {
+        return MinPurity > 0 && MaxPurity >= MinPurity && MaxPurity <= 1.0 && MinPloidy > 0 && MaxPloidy >= MinPloidy && PurityIncrement > 0;
+    }
+
     public static void addConfig(final ConfigBuilder configBuilder)
     {
         configBuilder.addDecimal(MIN_PURITY, "Minimum purity", MIN_PURITY_DEFAULT);

diff --git a/purple/src/main/java/com/hartwig/hmftools/purple/PurpleApplication.java b/purple/src/main/java/com/hartwig/hmftools/purple/PurpleApplication.java
@@ -113,6 +113,12 @@ private PurpleApplication(final ConfigBuilder configBuilder) throws IOException
 
     public void run()
     {
+        if(!mConfig.Fitting.hasValidValues())
+        {
+            PPL_LOGGER.error("invalid configured fitting values");
+            System.exit(1);
+        }
+
         long startTimeMs = System.currentTimeMillis();
 
         try

diff --git a/purple/src/main/java/com/hartwig/hmftools/purple/fitting/FittedPurityFactory.java b/purple/src/main/java/com/hartwig/hmftools/purple/fitting/FittedPurityFactory.java
@@ -1,6 +1,5 @@
 package com.hartwig.hmftools.purple.fitting;
 
-import static com.hartwig.hmftools.common.utils.Doubles.lessOrEqual;
 import static com.hartwig.hmftools.common.utils.Doubles.positiveOrZero;
 
 import java.util.Collection;
@@ -107,17 +106,29 @@ public void fitPurity() throws ExecutionException, InterruptedException
     {
         FittingConfig config = mConfig.Fitting;
 
+        List<Double> purityValues = Lists.newArrayList();
+
+        double purityValue = config.MinPurity;
+
+        while(purityValue <= config.MaxPurity)
+        {
+            purityValues.add(purityValue);
+
+            purityValue = Doubles.round(purityValue + config.PurityIncrement, 2);
+        }
+
         if(mConfig.Threads <= 1)
         {
-            for(double purity = config.MinPurity; lessOrEqual(purity, config.MaxPurity); purity += config.PurityIncrement)
+            for(Double purity : purityValues)
             {
                 mFittedPurities.addAll(fitPurity(purity));
             }
         }
         else
         {
             List<Future<List<FittedPurity>>> futures = Lists.newArrayList();
-            for(double purity = config.MinPurity; lessOrEqual(purity, config.MaxPurity); purity += config.PurityIncrement)
+
+            for(Double purity : purityValues)
             {
                 futures.add(mExecutorService.submit(callableFitPurity(purity)));
             }
@@ -261,8 +272,9 @@ private static List<Double> sequence(double inclusiveMin, double exclusiveMax, d
         while(Doubles.lessThan(ploidy, exclusiveMax))
         {
             results.add(ploidy);
-            ploidy += increment;
+            ploidy = Doubles.round(ploidy + increment, 2);
         }
+
         return results;
     }
 }
diff --git a/purple/src/main/java/com/hartwig/hmftools/purple/fittingsnv/SomaticKernelDensityPeaks.java b/purple/src/main/java/com/hartwig/hmftools/purple/fittingsnv/SomaticKernelDensityPeaks.java
@@ -1,6 +1,5 @@
 package com.hartwig.hmftools.purple.fittingsnv;
 
-import static java.lang.Math.abs;
 import static java.lang.String.format;
 
 import static com.hartwig.hmftools.purple.PurpleUtils.PPL_LOGGER;
@@ -45,7 +44,7 @@ public final class SomaticKernelDensityPeaks
             {
                 if(peak.Count >= minPeak)
                 {
-                    FittedPurity matchedFittedPurity = findMatchedFittedPurity(impliedPurity, allCandidates, 0.001);
+                    FittedPurity matchedFittedPurity = findMatchedFittedPurity(impliedPurity, allCandidates);
                     if(matchedFittedPurity != null)
                     {
                         PPL_LOGGER.debug("somatic implied purity({})", impliedPurity);
@@ -70,7 +69,7 @@ public final class SomaticKernelDensityPeaks
                 double impliedPurity = peak.AlleleFrequency * 2;
                 if(peak.Count == maxPeak)
                 {
-                    FittedPurity matchedFittedPurity = findMatchedFittedPurity(impliedPurity, allCandidates, 0.001);
+                    FittedPurity matchedFittedPurity = findMatchedFittedPurity(impliedPurity, allCandidates);
                     if(matchedFittedPurity != null)
                     {
                         PPL_LOGGER.debug("somatic implied purity({})", impliedPurity);

diff --git a/purple/src/main/java/com/hartwig/hmftools/purple/fittingsnv/SomaticPurityFitter.java b/purple/src/main/java/com/hartwig/hmftools/purple/fittingsnv/SomaticPurityFitter.java
@@ -30,8 +30,6 @@
 import java.util.Optional;
 import java.util.StringJoiner;
 
-import javax.swing.DefaultListSelectionModel;
-
 import com.hartwig.hmftools.common.drivercatalog.panel.DriverGene;
 import com.hartwig.hmftools.common.drivercatalog.panel.DriverGenePanel;
 import com.hartwig.hmftools.common.genome.chromosome.HumanChromosome;
@@ -67,7 +65,8 @@ public class SomaticPurityFitter
     private final int mKdMinPeak;
     private final int mKdMinSomatics;
 
-    public SomaticPurityFitter(int minPeak, int minSomatics, int minReadCount, int maxReadCount, double minPurity, double maxPurity)
+    public SomaticPurityFitter(
+            int minPeak, int minSomatics, int minReadCount, int maxReadCount, double minPurity, double maxPurity)
     {
         mMinReadCount = minReadCount;
         mMaxReadCount = maxReadCount;
@@ -225,7 +224,7 @@ public FittedPurity fromSomatics(
 
         if(hotspotPurity > somaticPeakPurity)
         {
-            FittedPurity matchedFittedPurity = findMatchedFittedPurity(hotspotPurity, diploidCandidates, 0.005);
+            FittedPurity matchedFittedPurity = findMatchedFittedPurity(hotspotPurity, diploidCandidates);
 
             if(matchedFittedPurity != null)
                 return matchedFittedPurity;
@@ -237,7 +236,7 @@ public FittedPurity fromSomatics(
 
             if(reassessmentPurity > somaticPeakPurity)
             {
-                FittedPurity matchedFittedPurity = findMatchedFittedPurity(reassessmentPurity, diploidCandidates, 0.005);
+                FittedPurity matchedFittedPurity = findMatchedFittedPurity(reassessmentPurity, diploidCandidates);
 
                 if(matchedFittedPurity != null)
                     return matchedFittedPurity;
@@ -254,12 +253,27 @@ public static boolean highlyDiploidSomaticOrPanel(final FittedPurity normalPurit
             && normalPurityFit.ploidy() < SOMATIC_FIT_TUMOR_ONLY_PLOIDY_MAX;
     }
 
-    protected static FittedPurity findMatchedFittedPurity(double purity, final List<FittedPurity> allCandidates, final double epsilon)
+    protected static FittedPurity findMatchedFittedPurity(double purity, final List<FittedPurity> allCandidates)
     {
-        return allCandidates.stream()
-                .filter(x -> abs(x.purity() - purity) < epsilon)
-                .filter(x -> abs(x.ploidy() - 2.0) < epsilon) // assumes ploidy 2 as well since is looking a highly diploid solutions
-                .findFirst().orElse(null);
+        // find the closest purity with diploid ploidy
+        FittedPurity closestPurity = null;
+        double closestDiff = 0;
+
+        for(FittedPurity fittedPurity : allCandidates)
+        {
+            if(abs(fittedPurity.ploidy() - 2) > 0.005)
+                continue;
+
+            double diff = abs(fittedPurity.purity() - purity);
+
+            if(closestPurity == null || diff < closestDiff)
+            {
+                closestDiff = diff;
+                closestPurity = fittedPurity;
+            }
+        }
+
+        return closestPurity;
     }
 
     @Nullable
@@ -351,7 +365,7 @@ else if(variantVafs.size() == 2)
 
         PPL_LOGGER.info("somatic VAF-based purity({}) from {} variants", formatPurity(somaticPurity), variantVafs.size());
 
-        FittedPurity matchedFittedPurity = findMatchedFittedPurity(somaticPurity, allCandidates, 0.005);
+        FittedPurity matchedFittedPurity = findMatchedFittedPurity(somaticPurity, allCandidates);
 
         if(matchedFittedPurity != null)
             return matchedFittedPurity;