Browse code

More accurate CIs in callMutationBurden.

Markus Riester authored on 05/11/2017 23:27:56
Showing 4 changed files

... ...
@@ -2,8 +2,8 @@ Package: PureCN
2 2
 Type: Package
3 3
 Title: Copy number calling and SNV classification using
4 4
     targeted short read sequencing
5
-Version: 1.9.1
6
-Date: 2017-10-22
5
+Version: 1.9.2
6
+Date: 2017-11-05
7 7
 Authors@R: c(person("Markus", "Riester", role=c("aut", "cre"),
8 8
         email="[email protected]"),
9 9
     person("Angad P.", "Singh", role="aut"))
... ...
@@ -112,6 +112,7 @@ importFrom(graphics,text)
112 112
 importFrom(rtracklayer,import)
113 113
 importFrom(stats,C)
114 114
 importFrom(stats,aggregate)
115
+importFrom(stats,binom.test)
115 116
 importFrom(stats,complete.cases)
116 117
 importFrom(stats,cutree)
117 118
 importFrom(stats,dbeta)
... ...
@@ -126,7 +127,6 @@ importFrom(stats,loess)
126 127
 importFrom(stats,pbeta)
127 128
 importFrom(stats,prcomp)
128 129
 importFrom(stats,predict)
129
-importFrom(stats,qpois)
130 130
 importFrom(stats,runif)
131 131
 importFrom(stats,t.test)
132 132
 importFrom(stats,weighted.mean)
... ...
@@ -1,3 +1,13 @@
1
+Changes in version 1.10.0
2
+------------------------
3
+
4
+NEW FEATURES
5
+
6
+    o More accurate confidence intervals in callMutationBurden
7
+
8
+
9
+    
10
+
1 11
 Changes in version 1.8.0
2 12
 ------------------------
3 13
 
... ...
@@ -50,7 +50,7 @@
50 50
 #'     callable=callableBed, exclude=exclude, fun.countMutation=myVcfFilter)
51 51
 #' 
52 52
 #' @export callMutationBurden
53
-#' @importFrom stats qpois
53
+#' @importFrom stats binom.test
54 54
 callMutationBurden <- function(res, id = 1, remove.flagged = TRUE, 
55 55
     min.prior.somatic=0.1, min.cellfraction=0, 
56 56
     fun.countMutation=function(vcf) width(vcf)==1,
... ...
@@ -121,18 +121,22 @@ callMutationBurden <- function(res, id = 1, remove.flagged = TRUE,
121 121
     )
122 122
     
123 123
     if (!is.na(ret$callable.bases.ontarget) && ret$callable.bases.ontarget > 0) {
124
-        lambda <- max(1, ret$somatic.ontarget)
125 124
         delta <- ret$callable.bases.ontarget/1e+6
125
+        ci <- binom.test(x = ret$somatic.ontarget, 
126
+            n = ret$callable.bases.ontarget)$conf.int
126 127
         ret <- cbind(ret, data.frame(
127
-            somatic.rate.ontarget=ret$somatic.ontarget/delta,
128
-            somatic.rate.ontarget.95.lower=qpois(0.025, lambda=lambda)/delta,
129
-            somatic.rate.ontarget.95.upper=qpois(0.975, lambda=lambda)/delta
128
+            somatic.rate.ontarget = ret$somatic.ontarget/delta,
129
+            somatic.rate.ontarget.95.lower = ci[1] * 1e+6,
130
+            somatic.rate.ontarget.95.upper = ci[2] * 1e+6
130 131
         ))
131
-        lambda <- max(1, ret$private.germline.ontarget)
132
+        # if multiple CI methods were provided, then ret contains multiple rows
133
+        # now
134
+        ci <- binom.test(x = ret$private.germline.ontarget, 
135
+            n = ret$callable.bases.ontarget)$conf.int
132 136
         ret <- cbind(ret, data.frame(
133
-            private.germline.rate.ontarget=ret$private.germline.ontarget/delta,
134
-            private.germline.rate.ontarget.95.lower=qpois(0.025, lambda=lambda)/delta,
135
-            private.germline.rate.ontarget.95.upper=qpois(0.975, lambda=lambda)/delta
137
+            private.germline.rate.ontarget = ret$private.germline.ontarget/delta,
138
+            private.germline.rate.ontarget.95.lower = ci[1] * 1e+6,
139
+            private.germline.rate.ontarget.95.upper = ci[2] * 1e+6
136 140
         ))
137 141
     }  
138 142