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
|
get_Poisson_conf_intervals = function(seq_range, quantile_vec = c(0.05,0.95), plot=T) {
num_rows = length(seq_range)
num_cols = length(quantile_vec)
m = matrix(nrow=num_rows, ncol=num_cols)
colnames(m) = quantile_vec
rownames(m) = seq_range
row_count = 0
for (i in seq_range) {
q = qpois(quantile_vec, i)
row_count = row_count + 1
m[row_count,] = q
}
results = list()
results$mat = m
if (plot) {
percents = plot_conf_intervals(m)
results$pct = percents
}
return(results)
}
get_NB_conf_intervals = function(seq_range, dispersion, quantile_vec = c(0.05,0.95), plot=T) {
size = 1/dispersion #according to the mu-definition of size in nbinom of R, where var = mean + (1/size)mean^2
num_rows = length(seq_range)
num_cols = length(quantile_vec)
m = matrix(nrow=num_rows, ncol=num_cols)
colnames(m) = quantile_vec
rownames(m) = seq_range
row_count = 0
for (i in seq_range) {
q = qnbinom(quantile_vec, mu=i, size=size)
row_count = row_count + 1
m[row_count,] = q
}
results = list()
results$mat = m
if (plot) {
percents = plot_conf_intervals(m)
results$pct = percents
}
return(results)
}
plot_conf_intervals = function(m) {
par(mfrow=c(1,2))
c_names = colnames(m)
r_vals = as.numeric(rownames(m))
max_val = max(m)
plot(r_vals, r_vals, xlab="known read counts", ylab="Poisson read counts dist", t='l')
# plot the confidence levels
line_colors = rainbow(length(c_names))
for(i in 1:length(c_names)) {
points(r_vals,m[,i], t='l', col=line_colors[i])
}
# plot the max percentage of value for 95% conf level.
percents=c()
for (i in 1:length(r_vals)) {
max_delta = max(abs(m[i,]-r_vals[i]))
percent = max_delta/r_vals[i]*100
percents[i]=percent;
}
plot(r_vals, percents, ylim=c(0,100), xlab="read counts", ylab="percent of value for 95% conf interval")
percents
}
|
