Major refactoring. Use preset error model by default. Change preproce…

…ssor default parameters.
mikessh · Nov 27, 2016 · 07698d4 · 07698d4
1 parent 282171b
commit 07698d4
Show file tree

Hide file tree

Showing 16 changed files with 323 additions and 179 deletions.
diff --git a/src/main/java/com/antigenomics/mageri/core/assemble/AssemblerParameters.java b/src/main/java/com/antigenomics/mageri/core/assemble/AssemblerParameters.java
@@ -30,13 +30,13 @@ public final class AssemblerParameters implements ParameterSet {
     public static AssemblerParameters DEFAULT = new AssemblerParameters(
             4, 8, 4, 0,
             0.3, 0.0, 0.3, 0.0,
-            1e-3,
+            1e-2,
             false, true, true, false);
 
     public static AssemblerParameters TORRENT454 = new AssemblerParameters(
             4, 8, 4, 2,
             0.7, 0.3, 0.3, 0.5,
-            1e-3,
+            1e-2,
             true, true, false, false);
 
     public AssemblerParameters(int offsetRange, int anchorRegion, int maxMMs, int maxConsequentMMs,

diff --git a/src/main/java/com/antigenomics/mageri/core/input/PreprocessorParameters.java b/src/main/java/com/antigenomics/mageri/core/input/PreprocessorParameters.java
@@ -27,14 +27,14 @@ public class PreprocessorParameters implements ParameterSet {
     private final double minUmiMismatchRatio;
 
     public static final PreprocessorParameters DEFAULT = new PreprocessorParameters(
-            QualityDefaults.PH33_BAD_QUAL, QualityDefaults.PH33_GOOD_QUAL,
+            QualityDefaults.PH33_BAD_QUAL, (byte) 25,
             true,
-            20.0, false, 5);
+            10.0, true, 5);
 
     public static final PreprocessorParameters IGNORE_QUAL = new PreprocessorParameters(
             (byte) 0, (byte) 0,
             true,
-            20.0, false, 5);
+            10.0, true, 5);
 
     public PreprocessorParameters(byte umiQualThreshold, byte goodQualityThreshold,
                                   boolean trimAdapters,

diff --git a/src/main/java/com/antigenomics/mageri/core/variant/VariantCaller.java b/src/main/java/com/antigenomics/mageri/core/variant/VariantCaller.java
@@ -92,19 +92,18 @@ public VariantCaller(ConsensusAligner consensusAligner, MinorCaller minorCaller,
                         int majorCount = mutationsTable.getMajorMigCount(pos, to),
                                 coverage = mutationsTable.getMigCoverage(pos);
 
-                        ErrorRateEstimate errorRateEstimate = errorModel.computeErrorRate(mutation);
-                        double score = getQScore(majorCount, coverage,
-                                errorRateEstimate.getErrorRate(), variantCallerParameters);
+                        VariantQuality variantQuality = errorModel.computeQuality(majorCount,
+                                coverage, mutation);
 
                         NucleotideSequenceBuilder nsb = new NucleotideSequenceBuilder(1);
                         nsb.setCode(0, mutationsTable.getAncestralBase(pos));
 
                         variant = new Variant(reference,
                                 mutation, majorCount,
                                 mutationsTable.getMigCoverage(pos),
-                                score, mutationsTable.getMeanCqs(pos, to),
+                                variantQuality.getScore(), mutationsTable.getMeanCqs(pos, to),
                                 nsb.create(), mutationsTable.hasReferenceBase(pos),
-                                errorRateEstimate);
+                                variantQuality.getErrorRateEstimate());
                     } else if (variantCallerParameters.isNoIndels()) {
                         continue;
                     } else {
@@ -167,58 +166,6 @@ public VariantCaller(ConsensusAligner consensusAligner, MinorCaller minorCaller,
         Collections.sort(variants);
     }
 
-    static double getBinomialQScore(int majorCount, int total, double errorRate) {
-        BinomialDistribution binomialDistribution = new BinomialDistributionImpl(total,
-                errorRate);
-
-        double score;
-        try {
-            score = -10 * Math.log10(1.0 - binomialDistribution.cumulativeProbability(majorCount) +
-                    0.5 * binomialDistribution.probability(majorCount));
-        } catch (MathException e) {
-            e.printStackTrace();
-            return VcfUtil.UNDEF_QUAL;
-        }
-
-        return score;
-    }
-
-    static double getNegBinomialQScore(int majorCount, int total, double errorRate,
-                                       VariantCallerParameters variantCallerParameters) {
-        return getNegBinomialQScore(majorCount, total, errorRate,
-                variantCallerParameters.getCompoundQScoreMu(),
-                variantCallerParameters.getCompoundQScoreSD());
-    }
-
-    static double getNegBinomialQScore(int majorCount, int total, double errorRate,
-                                       double modelMu, double modelSD) {
-        double r = 1.0 / (Math.exp(modelSD * modelSD) - 1),
-                p = 1 / (1 + r / Math.exp(modelMu + modelSD * modelSD / 2) / errorRate / total);
-
-        if (r <= 0 || p >= 1 || p <= 0) {
-            return VcfUtil.UNDEF_QUAL;
-        }
-
-        return -10 * Math.log10(1.0 - AuxiliaryStats.negativeBinomialCdf(majorCount, r, p));
-    }
-
-    static double getQScore(int majorCount, int total, double errorRate,
-                            VariantCallerParameters variantCallerParameters) {
-        if (majorCount == 0) {
-            return 0;
-        }
-
-        if (Double.isNaN(errorRate) || errorRate == 0) {
-            return VcfUtil.UNDEF_QUAL;
-        }
-
-        double score = variantCallerParameters.getErrorModelType() == ErrorModelType.MinorBased ?
-                getNegBinomialQScore(majorCount, total, errorRate, variantCallerParameters) :
-                getBinomialQScore(majorCount, total, errorRate);
-
-        return Double.isInfinite(score) ? VcfUtil.MAX_QUAL : score;
-    }
-
     public ReferenceLibrary getReferenceLibrary() {
         return referenceLibrary;
     }