Changes to MBEM, fix tests

mikessh · Nov 12, 2016 · b0bd81b · b0bd81b
1 parent 17a26bf
commit b0bd81b
Show file tree

Hide file tree

Showing 7 changed files with 20 additions and 13 deletions.
diff --git a/src/main/java/com/antigenomics/mageri/core/assemble/DummyMinorCaller.java b/src/main/java/com/antigenomics/mageri/core/assemble/DummyMinorCaller.java
@@ -81,7 +81,7 @@ public int getTotalMigs() {
     }
 
     @Override
-    public int getGeometricMeanMigSize() {
+    public double getGeometricMeanMigSize() {
         return (int) Math.exp(logMigSize.get() / getTotalMigs());
     }
 }
diff --git a/src/main/java/com/antigenomics/mageri/core/assemble/MinorCaller.java b/src/main/java/com/antigenomics/mageri/core/assemble/MinorCaller.java
@@ -37,5 +37,5 @@ public abstract class MinorCaller<M extends MinorCaller> extends PipelineBlock {
 
     public abstract int getTotalMigs();
 
-    public abstract int getGeometricMeanMigSize();
+    public abstract double getGeometricMeanMigSize();
 }
diff --git a/src/main/java/com/antigenomics/mageri/core/assemble/PoissonTestMinorCaller.java b/src/main/java/com/antigenomics/mageri/core/assemble/PoissonTestMinorCaller.java
@@ -158,8 +158,8 @@ public int getTotalMigs() {
     }
 
     @Override
-    public int getGeometricMeanMigSize() {
-        return (int) Math.exp(logMigSize.get() / getTotalMigs());
+    public double getGeometricMeanMigSize() {
+        return Math.exp(logMigSize.get() / getTotalMigs());
     }
 
     private int getM(int from, int to) {

diff --git a/src/main/java/com/antigenomics/mageri/core/variant/model/MinorBasedErrorModel.java b/src/main/java/com/antigenomics/mageri/core/variant/model/MinorBasedErrorModel.java
@@ -87,8 +87,10 @@ public static double computePropagateProb(double efficiency, double order) {
     }
 
     public static double computeBaseErrorRateEstimate(double minorRate, double fdr,
-                                                      double readFractionEstForCalledMinors) {
-        return minorRate * (1.0 - fdr) * readFractionEstForCalledMinors;
+                                                      double readFractionEstForCalledMinors,
+                                                      double geomMeanMigSize) {
+        minorRate *= (1.0 - fdr);
+        return minorRate * (readFractionEstForCalledMinors + 1.0 / geomMeanMigSize);
     }
 
     @Override
@@ -118,10 +120,11 @@ public ErrorRateEstimate computeErrorRate(int pos, int from, int to) {
         // Share of minors that are misidentified sequencing errors
         double fdr = minorCaller.computeFdr(from, to);
 
-        double errorRateBase = computeBaseErrorRateEstimate(minorRate, fdr, readFractionForCalledMinors);
+        double errorRateBase = computeBaseErrorRateEstimate(minorRate, fdr,
+                readFractionForCalledMinors, minorCaller.getGeometricMeanMigSize());
 
         // Share of minors to expect in the absence of sequencing errors due to sampling
-        double recall = Math.exp(-errorRateBase * minorCaller.getGeometricMeanMigSize());
+        double recall = 1.0 - Math.exp(-errorRateBase * minorCaller.getGeometricMeanMigSize());
 
         // Adjust for probability of error propagation
         double majorErrorRate = errorRateBase / cycles / lambda * (1.0 + lambda) * propagateProb;

diff --git a/src/test/java/com/antigenomics/mageri/core/assemble/AssemblerMinorTest.java b/src/test/java/com/antigenomics/mageri/core/assemble/AssemblerMinorTest.java
@@ -180,7 +180,11 @@ public void exactTestIndel() {
                 new DummyMinorCaller());
         SConsensus consensus = assembler.assemble(mig);
 
-        Assert.assertEquals("Incorrect consensus sequence", cons, consensus.getConsensusSQPair().getSequence().toString());
+        System.out.println(cons);
+        System.out.println(consensus.getConsensusSQPair().getSequence());
+
+        Assert.assertEquals("Incorrect consensus sequence", cons,
+                consensus.getConsensusSQPair().getSequence().toString());
 
         System.out.println("Expected minors:");
         for (int minor : minors) {

diff --git a/src/test/java/com/antigenomics/mageri/core/assemble/MinorCallerTest.java b/src/test/java/com/antigenomics/mageri/core/assemble/MinorCallerTest.java
@@ -126,7 +126,7 @@ private MinorCaller simulationTest(int nCycles, double lambda,
 
         double minorRate = (calledTruePCRMinors + calledFalsePCRMinors) / (double) (nMigs * readLength);
         double errorRateEst = MinorBasedErrorModel.computeBaseErrorRateEstimate(minorRate,
-                minorCaller.computeFdr(0, 0), minorCaller.getReadFractionForCalledMinors(0, 0)) /
+                minorCaller.computeFdr(0, 0), minorCaller.getReadFractionForCalledMinors(0, 0), migSize) /
                 nCycles / lambda * (1.0 + lambda);
         System.out.println("Error rate est = " + errorRateEst);
         Assert.assertTrue("No more than order of magnitude difference between expected PCR " +

diff --git a/src/test/java/com/antigenomics/mageri/core/variant/model/MinorBasedErrorModelTest.java b/src/test/java/com/antigenomics/mageri/core/variant/model/MinorBasedErrorModelTest.java
@@ -40,11 +40,11 @@ public class MinorBasedErrorModelTest {
     @Test
     @Category(ComplexRandomTests.class)
     public void test() {
-        int nMigs = 30000, migSize = 100;
+        int nMigs = 20000, migSize = 100;
 
         RandomReferenceGenerator randomReferenceGenerator = new RandomReferenceGenerator();
-        randomReferenceGenerator.setReferenceSizeMin(100);
-        randomReferenceGenerator.setReferenceSizeMax(100);
+        randomReferenceGenerator.setReferenceSizeMin(50);
+        randomReferenceGenerator.setReferenceSizeMax(50);
 
         ReferenceLibrary referenceLibrary = randomReferenceGenerator.nextReferenceLibrary(1);
         Reference reference = referenceLibrary.getAt(0);