-
Notifications
You must be signed in to change notification settings - Fork 0
/
1D Pooling simulation 2.Rmd
154 lines (117 loc) · 4.36 KB
/
1D Pooling simulation 2.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
---
title: '1D Pooling: Simulation'
author: "Xianbin Cheng"
date: "December 17, 2018"
output: html_document
---
# Objective
1. We want to simulate one-dimensional pooling.
2. We want to analyze the cost.
# Method
1. Load necessary libraries and files.
```{r, message=FALSE, warning = FALSE}
library(tidyverse)
library(Hmisc)
library(mc2d)
source("1D_pooling_simulation.R")
source("STD_simulation.R")
```
2. Set input parameters
* `n` = the total number of kernels in a plate (48 by default).
* `n_pool` = the number of pools.
*
```{r}
#n = 48
#n_pos_vals = 1:(n-1)
thresh = 20
n_iter = 1000
by = "column"
```
3. Generate simulated data. Calculate sensitivity, specificity, and the total number of tests needed. The function used to generate random contamination levels is the same as the one used to generate data for STD simulation.
```{r, eval = FALSE}
### 48-well
# Pool by rows
set.seed(123)
pool_row_48 = tune_1d_n_pos(n_pos_vals = 1:(48-1), n_iter = n_iter, n = 48, thresh = thresh, by = "row")
saveRDS(object = pool_row_48, file = "Pool_row_48_1000.rds")
# Pool by columns
set.seed(123)
pool_col_48 = tune_1d_n_pos(n_pos_vals = 1:(48-1), n_iter = n_iter, n = 48, thresh = thresh, by = "column")
saveRDS(object = pool_col_48, file = "Pool_col_48_1000.rds")
### 96-well
# Pool by rows
set.seed(123)
pool_row_96 = tune_1d_n_pos(n_pos_vals = 1:(96-1), n_iter = n_iter, n = 96, thresh = thresh, by = "row")
saveRDS(object = pool_row_96, file = "Pool_row_96_1000.rds")
# Pool by columns
set.seed(123)
pool_col_96 = tune_1d_n_pos(n_pos_vals = 1:(96-1), n_iter = n_iter, n = 96, thresh = thresh, by = "column")
saveRDS(object = pool_col_96, file = "Pool_col_96_1000.rds")
```
```{r, echo = FALSE}
pool_row_48 = readRDS("Pool_row_48_1000.rds")
pool_col_48 = readRDS("Pool_col_48_1000.rds")
pool_row_96 = readRDS("Pool_row_96_1000.rds")
pool_col_96 = readRDS("Pool_col_96_1000.rds")
```
```{r}
### 48-well
# Pool by rows
metrics_row_48 = do.call(what = rbind.data.frame, args = pool_row_48) %>%
gather(data = ., key = "Type", value = "Value", - n_pos)
cost_row_48 = calc_1d_cost(data = pool_row_48, n = 48, by = "row")
# Pool by columns
metrics_col_48 = do.call(what = rbind.data.frame, args = pool_col_48) %>%
gather(data = ., key = "Type", value = "Value", - n_pos)
cost_col_48 = calc_1d_cost(data = pool_col_48, n = 48, by = "column")
### 96-well
# Pool by rows
metrics_row_96 = do.call(what = rbind.data.frame, args = pool_row_96) %>%
gather(data = ., key = "Type", value = "Value", - n_pos)
cost_row_96 = calc_1d_cost(data = pool_row_96, n = 96, by = "row")
# Pool by columns
metrics_col_96 = do.call(what = rbind.data.frame, args = pool_col_96) %>%
gather(data = ., key = "Type", value = "Value", - n_pos)
cost_col_96 = calc_1d_cost(data = pool_col_96, n = 96, by = "column")
```
# Results
1. Visualization of metrics with simulated data.
```{r}
### 48-well
# Pool by rows
draw_metrics(df = metrics_row_48, n = 48, method = "median_hilow")
# Pool by columns
draw_metrics(df = metrics_col_48, n = 48, method = "median_hilow")
### 96-well
# Pool by rows
draw_metrics(df = metrics_row_96, n = 96, method = "median_hilow")
# Pool by columns
draw_metrics(df = metrics_col_96, n = 96, method = "median_hilow")
```
2. Cost analysis.
1) The total number of tests needed:
```{r}
### 48-well
# Pool by rows
draw_cost(df = cost_row_48, n = 48, var = n_test_total, ylab = "The total number of tests needed")
# Pool by columns
draw_cost(df = cost_col_48, n = 48, var = n_test_total, ylab = "The total number of tests needed")
### 96-well
# Pool by rows
draw_cost(df = cost_row_96, n = 96, var = n_test_total, ylab = "The total number of tests needed")
# Pool by columns
draw_cost(df = cost_col_96, n = 96, var = n_test_total, ylab = "The total number of tests needed")
```
2) The number of pipetting:
```{r}
### 48-well
# Pool by rows
draw_cost(df = cost_row_48, n = 48, var = n_pipette, ylab = "The number of transfers")
# Pool by columns
draw_cost(df = cost_col_48, n = 48, var = n_pipette, ylab = "The number of transfers")
### 96-well
# Pool by rows
draw_cost(df = cost_row_96, n = 96, var = n_pipette, ylab = "The number of transfers")
# Pool by columns
draw_cost(df = cost_col_96, n = 96, var = n_pipette, ylab = "The number of transfers")
```