From 7cba23beac756768400d5f90d234dc58cb45e632 Mon Sep 17 00:00:00 2001
From: ctokheim <collintokheim@gmail.com>
Date: Sun, 29 May 2016 19:15:05 -0400
Subject: [PATCH] Change total mutation count calculation

---
 R/ratioMetric.R | 78 ++++++++++++++++++++++++++-----------------------
 1 file changed, 42 insertions(+), 36 deletions(-)

diff --git a/R/ratioMetric.R b/R/ratioMetric.R
index e4c07ce..bf1de58 100644
--- a/R/ratioMetric.R
+++ b/R/ratioMetric.R
@@ -33,24 +33,27 @@ ratiometric.binom.power <- function(p, N, mu,
   for(i in N){
     # step one, find the # of mutations where
     # it is expected to occur at least 90% of the time
-    j <- 1
-    while(j){
-      prob <- pbinom(j-1, round(Leff*i), muEffect)
-      if(prob >= .1){
-        mutEffect <- j
-        break
-      }
-      j <- j+1
-    }
-    j <- 1
-    while(j){
-      prob <- pbinom(j-1, round((1-nonsil.frac)*1/nonsil.frac*Leff*i), mu)
-      if(prob >= .1){
-        mutEffect <- mutEffect + j
-        break
-      }
-      j <- j+1
-    }
+    #j <- 1
+    #while(j){
+      #prob <- pbinom(j-1, round(Leff*i), muEffect)
+      #if(prob >= .1){
+        #mutEffect <- j
+        #break
+      #}
+      #j <- j+1
+    #}
+    #j <- 1
+    #while(j){
+      #prob <- pbinom(j-1, round((1-nonsil.frac)*1/nonsil.frac*Leff*i), mu)
+      #if(prob >= .1){
+        #mutEffect <- mutEffect + j
+        #break
+      #}
+      #j <- j+1
+    #}
+    totRiskBases = 1/nonsil.frac*Leff*i
+    nonsilRiskBases = Leff*i
+    mutEffect <- ceiling(nonsilRiskBases*muEffect + (totRiskBases-nonsilRiskBases)*mu)
     
     # step two, find critical threshold
     j <- 1
@@ -109,24 +112,27 @@ ratiometric.bbd.power <- function(my.alpha, my.beta,
   for(i in N){
     # step one, find the # of mutations where
     # it is expected to occur at least 90% of the time
-    j <- 1
-    while(j){
-      prob <- pbinom(j-1, round(Leff*i), muEffect)
-      if(prob >= .1){
-        mutEffect <- j
-        break
-      }
-      j <- j+1
-    }
-    j <- 1
-    while(j){
-      prob <- pbinom(j-1, round((1-nonsil.frac)*1/nonsil.frac*Leff*i), mu)
-      if(prob >= .1){
-        mutEffect <- mutEffect + j
-        break
-      }
-      j <- j+1
-    }
+    #j <- 1
+    #while(j){
+      #prob <- pbinom(j-1, round(Leff*i), muEffect)
+      #if(prob >= .1){
+        #mutEffect <- j
+        #break
+      #}
+      #j <- j+1
+    #}
+    #j <- 1
+    #while(j){
+      #prob <- pbinom(j-1, round((1-nonsil.frac)*1/nonsil.frac*Leff*i), mu)
+      #if(prob >= .1){
+        #mutEffect <- mutEffect + j
+        #break
+      #}
+      #j <- j+1
+    #}
+    totRiskBases = 1/nonsil.frac*Leff*i
+    nonsilRiskBases = Leff*i
+    mutEffect <- ceiling(nonsilRiskBases*muEffect + (totRiskBases-nonsilRiskBases)*mu)
     
     # step two, find critical threshold
     j <- 1