-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy path.Rhistory
executable file
·512 lines (512 loc) · 21.2 KB
/
.Rhistory
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
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
ggtitle(bquote(Distribution~of~Mixture~Model~Parameters*","~alpha~"="~.(round(alpha1,2))))
x
document()
install("../oncomix/")
install()
library(devtools)
build_vignettes()
?stats
build_vignettes()
scatterMixPlot <- function(mmParams, selIndThresh = 1, gene_labels = NULL){
mmParams = as.data.frame(mmParams)
one_over_alpha = diff(range(mmParams$deltaMu2))
alpha1 = 1/one_over_alpha
quants = c(0.01, 0.10, 0.50, 0.90, 0.99) #add in the quantiles
colors_red=RColorBrewer::brewer.pal(n=length(quants), name="Reds")
deltaMu2Quant <- stats::quantile(mmParams[,"deltaMu2"], quants)
deltaMu1Quant <- stats::quantile(1/(abs(mmParams[,"deltaMu1"]) + alpha1), quants)
deltaMu2 <- deltaMu1 <- NULL
x = ggplot(data = as.data.frame(mmParams), aes(x = deltaMu2, y = 1/(abs(deltaMu1) + alpha1))) +
theme_classic() +
geom_hline(yintercept = deltaMu1Quant, col = colors_red, size = c(1,1,1,1,1)) +
geom_vline(xintercept = deltaMu2Quant, col = colors_red, size = c(1,1,1,1,1)) +
geom_point(alpha= 0.5) +
xlab(expression(paste(Delta, mu[2]))) +
ylab(expression(paste(frac(1, paste(Delta, mu[1], " + ", alpha)))))
if(selIndThresh < 1){
mmParams.si = mmParams[mmParams$SI > selIndThresh,]
x = x + geom_point(data = as.data.frame(mmParams.si),
aes(x = deltaMu2, y = 1/(abs(deltaMu1)+alpha1)),
size = 10, alpha=0.1,
col=colors_red[length(colors_red)],
fill=colors_red[length(colors_red)]) +
ggtitle(bquote(Distribution~of~Mixture~Model~Parameters*","~alpha~"="~.(round(alpha1,2))*", SI >"~.(selIndThresh)))
} else if(!is.null(gene_labels)){
mmParams.si = mmParams[gene_labels,]
mmParams.si$gene_labels = gene_labels
x = x + geom_point(data = as.data.frame(mmParams.si),
aes(x = deltaMu2, y = 1/(abs(deltaMu1)+alpha1)),
size = 10, alpha=0.1,
col=colors_red[length(colors_red)],
fill=colors_red[length(colors_red)]) +
ggrepel::geom_text_repel(data = mmParams.si, aes(x = deltaMu2, y = 1/(abs(deltaMu1)+alpha1)),
label = rownames(mmParams.si)) +
ggtitle(bquote(Distribution~of~Mixture~Model~Parameters*","~alpha~"="~.(round(alpha1,2))))
} else{
x = x + ggtitle(bquote(Distribution~of~Mixture~Model~Parameters*","~alpha~"="~.(round(alpha1,2))))
}
return(x)
}
library(ggrepel)
scatterMixPlot(mmParams, selIndThresh = 1, gene_labels = NULL)
document()
deltaMu2Quant <- stats::quantile(mmParams[,"deltaMu2"], quants)
mmParams = as.data.frame(mmParams)
one_over_alpha = diff(range(mmParams$deltaMu2))
alpha1 = 1/one_over_alpha
quants = c(0.01, 0.10, 0.50, 0.90, 0.99) #add in the quantiles
colors_red=RColorBrewer::brewer.pal(n=length(quants), name="Reds")
deltaMu2Quant <- stats::quantile(mmParams[,"deltaMu2"], quants)
deltaMu2Quant
build_vignettes()
geom_text_repel
document()
build_vignette()
build_vignettes()
devtools::check()
document()
library(devtools)
document ()
build_vignettes()
getwd()
devtools::check()
document()
document()
devtools::check()
library(ggrepel)
setwd("../")
check("oncomix")
check("oncomix")
df3 = data.frame(cbind("expr" = c(rnorm(113, 3), c(rnorm(56,3), rnorm(57, mean =6)))), "type"=as.factor(c(rep(2, 113), rep(1,113))))
df3
check("oncomix") d1 = data.frame(cbind("expr" = c(rnorm(113, 3), c(rnorm(56,3), rnorm(57, mean =6)))), "type"=as.factor(c(rep(2, 113), rep(1,113))))
d1 = data.frame(cbind("expr" = c(rnorm(113, 3), c(rnorm(56,3), rnorm(57, mean =6)))), "type"=as.factor(c(rep(2, 113), rep(1,113))))
ggplot(d1, aes(x="expr", color="type", fill="type", group="type")) + theme_classic() + theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank(),legend.position = "none", plot.title = element_text(size = 12), axis.text=element_text(size=8), axis.title=element_text(size=8), axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())+ ggtitle("Theoretical") + xlim(-0.2,10) +
stat_function(fun = dnorm, colour = "#00BFC4", args = list(mean = oe_means[1], sd=1), size = 5) +
stat_function(fun = dnorm, colour = "#00BFC4", args = list(mean = oe_means[2], sd=1), size = 5)
oe_means=c(3,7)
ggplot(d1, aes(x="expr", color="type", fill="type", group="type")) + theme_classic() + theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank(),legend.position = "none", plot.title = element_text(size = 12), axis.text=element_text(size=8), axis.title=element_text(size=8), axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())+ ggtitle("Theoretical") + xlim(-0.2,10) +
stat_function(fun = dnorm, colour = "#00BFC4", args = list(mean = oe_means[1], sd=1), size = 5) +
stat_function(fun = dnorm, colour = "#00BFC4", args = list(mean = oe_means[2], sd=1), size = 5)
df3 = data.frame(cbind("expr" = c(rnorm(113, 3), c(rnorm(56,3), rnorm(57, mean =6)))), "type"=as.factor(c(rep(2, 113), rep(1,113))))
ggplot(df3, aes(x=expr, color=type, fill=type, group=type)) +
theme_classic() +
theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank(),
legend.position = "none", plot.title = element_text(size = 12), axis.text=element_text(size=8),
axis.title=element_text(size=8), axis.title.x=element_blank(), axis.text.x=element_blank(),
axis.ticks.x=element_blank()) +
ggtitle("Theoretical") + xlim(oe_means[1]-3.2,oe_means[2] + 3) +
stat_function(fun = "dnorm", colour = "#F8766D", args = list(mean = 3.1, sd=0.7), size=3) +
stat_function(fun = "dnorm", colour = "#F8766D", args = list(mean=2.9, sd=0.7), size = 3) +
stat_function(fun = "dnorm", colour = "#00BFC4", args = list(mean = oe_means[1], sd=1), size = 3) +
stat_function(fun = "dnorm", colour = "#00BFC4", args = list(mean = oe_means[2], sd=1), size = 3)
d1 = data.frame(cbind("expr" = c(rnorm(113, 3), c(rnorm(56,3), rnorm(57, mean =6)))), "type"=as.factor(c(rep(2, 113), rep(1,113))))
type <- expr <- NULL
ggplot(d1, aes(x=expr, color=type, fill=type, group=type)) + theme_classic() + theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank(),legend.position = "none", plot.title = element_text(size = 12), axis.text=element_text(size=8), axis.title=element_text(size=8), axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())+ ggtitle("Theoretical") + xlim(-0.2,10) +
stat_function(fun = dnorm, colour = "#00BFC4", args = list(mean = oe_means[1], sd=1), size = 5) +
stat_function(fun = dnorm, colour = "#00BFC4", args = list(mean = oe_means[2], sd=1), size = 5)
df3 = data.frame(cbind("expr" = c(rnorm(113, 3), c(rnorm(56,3), rnorm(57, mean =6)))), "type"=as.factor(c(rep(2, 113), rep(1,113))))
type <- expr <- NULL
ggplot(df3, aes(x=expr, color=type, fill=type, group=type)) +
theme_classic() +
theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank(),
legend.position = "none", plot.title = element_text(size = 12), axis.text=element_text(size=8),
axis.title=element_text(size=8), axis.title.x=element_blank(), axis.text.x=element_blank(),
axis.ticks.x=element_blank()) +
ggtitle("Theoretical") + xlim(oe_means[1]-3.2,oe_means[2] + 3) +
stat_function(fun = "dnorm", colour = "#F8766D", args = list(mean = 3.1, sd=0.7), size=3) +
stat_function(fun = "dnorm", colour = "#F8766D", args = list(mean=2.9, sd=0.7), size = 3) +
stat_function(fun = "dnorm", colour = "#00BFC4", args = list(mean = oe_means[1], sd=1), size = 3) +
stat_function(fun = "dnorm", colour = "#00BFC4", args = list(mean = oe_means[2], sd=1), size = 3)
check("oncomix/")
check("oncomix")
Sys.which(Sys.getenv("R_QPDF", "qpdf"))
Sys.getenv("PATH")
check("oncomix")
check("oncomix")
check("oncomix")
check("oncomix")
check("oncomix")
setwd("..")
library(devtools)
check("oncomix")
seq_len(50)
dfNml <- as.data.frame(matrix(data=rgamma(n=150, shape=2, rate=2),
nrow=10, ncol=15))
colnames(dfNml) <- paste0("patientN", seq_len(ncol(dfNml)))
rownames(dfNml) <- paste0("gene", seq_len(nrow(dfNml))
dfTumor <- as.data.frame(matrix(data=rgamma(n=150, shape=4, rate=3),
nrow=15, ncol=10))
colnames(dfTumor) <- paste0("patientT", seq_len(ncol(dfTumor))
rownames(dfTumor) <- paste0("gene", seq_len(nrow(dfTumor))
dfNml <- as.data.frame(matrix(data=rgamma(n=150, shape=2, rate=2),
nrow=10, ncol=15))
colnames(dfNml) <- paste0("patientN", seq_len(ncol(dfNml)))
rownames(dfNml) <- paste0("gene", seq_len(nrow(dfNml))
nrows <- 200; ncols <- 6
counts <- matrix(runif(nrows * ncols, 1, 1e4), nrows)
colData <- DataFrame(Treatment=rep(c("ChIP", "Input"), 3),
row.names=LETTERS[1:6])
se0 <- SummarizedExperiment(assays=SimpleList(counts=counts),
colData=colData)
library(SummarizedExperiment)
nrows <- 200; ncols <- 6
counts <- matrix(runif(nrows * ncols, 1, 1e4), nrows)
colData <- DataFrame(Treatment=rep(c("ChIP", "Input"), 3),
row.names=LETTERS[1:6])
se0 <- SummarizedExperiment(assays=SimpleList(counts=counts),
colData=colData)
SummarizedExperiment
class(SummarizedExperiment)
standardGeneric(SummarizedExperiment)
class(se0)
class(se0) == "SummarizedExperiment"
assay(se0)
assay(SummarizedExperiment())
assay(dfNml)
assay(se0)
class(dfNml)
mixModelParams <- function(dfNml, dfTumor) {
if(class(dfNml) == "SummarizedExperiment"){
dfNml <- assay(dfNml)
}
if(class(dfTumor) == "SummarizedExperiment"){
dfTumor <- assay(dfTumor)
}
paramsNormal <- apply(dfNml, 1, function(x) {
y <- mclust::Mclust(data=x, G=2, modelNames="E", verbose=FALSE)
z <- c(nMu=y$parameters$mean,
nVar=y$parameters$variance$sigmasq,
nPi1=y$parameters$pro[1])
return(z)})
paramsTumor <- apply(dfTumor, 1, function(x) {
y <- mclust::Mclust(data=x, G=2, modelNames="E", verbose=FALSE)
z <- c(tMu=y$parameters$mean,
tVar=y$parameters$variance$sigmasq,
tPi1=y$parameters$pro[1])
return(z)})
params <- rbind(paramsNormal, paramsTumor)
rownames(params)[c(1:2, 5:6)] = c("nMu1","nMu2","tMu1", "tMu2")
deltaMu2 <- params["tMu2",]-params["nMu2",]
deltaMu1 <- params["tMu1",]-params["nMu1",]
boundaryTumor <- (params["tMu2",] + params["tMu1",]) / 2
siCalc <- function(vectNml, boundTumor){ #calculate the sel. idx
x <- sum(boundTumor > vectNml) / length(vectNml)
return(x)
}
SI <- sapply(seq_len(dfNml),
function(i) siCalc(dfNml[,i], boundaryTumor[i]))
params <- rbind(params, deltaMu2, deltaMu1, SI)
mmParamsDf <- data.frame(t(params))
score <- NA
mmParamsDf$score <- mmParamsDf$SI*{
(mmParamsDf$deltaMu2-mmParamsDf$deltaMu1)-
(mmParamsDf$nVar+mmParamsDf$tVar)}
mmParamsDfS <- mmParamsDf[with(mmParamsDf, order(-score)), ]
return(mmParamsDfS)
}
mmParams <- mixModelParams(dfNml, dfTumor)
dfTumor <- as.data.frame(matrix(data=rgamma(n=150, shape=4, rate=3),
nrow=15, ncol=10))
colnames(dfTumor) <- paste0("patientT", seq_len(ncol(dfTumor))
rownames(dfTumor) <- paste0("gene", seq_len(nrow(dfTumor))
dfTumor <- as.data.frame(matrix(data=rgamma(n=150, shape=4, rate=3),
nrow=15, ncol=10))
colnames(dfTumor) <- paste0("patientT", seq_len(ncol(dfTumor))
rownames(dfTumor) <- paste0("gene", seq_len(nrow(dfTumor))
dfTumor <- as.data.frame(matrix(data=rgamma(n=150, shape=4, rate=3),
nrow=15, ncol=10))
colnames(dfTumor) <- paste0("patientT", seq_len(ncol(dfTumor)))
rownames(dfTumor) <- paste0("gene", seq_len(nrow(dfTumor)))
mixModelParams <- function(dfNml, dfTumor) {
if(class(dfNml) == "SummarizedExperiment"){
dfNml <- assay(dfNml)
}
if(class(dfTumor) == "SummarizedExperiment"){
dfTumor <- assay(dfTumor)
}
paramsNormal <- apply(dfNml, 1, function(x) {
y <- mclust::Mclust(data=x, G=2, modelNames="E", verbose=FALSE)
z <- c(nMu=y$parameters$mean,
nVar=y$parameters$variance$sigmasq,
nPi1=y$parameters$pro[1])
return(z)})
paramsTumor <- apply(dfTumor, 1, function(x) {
y <- mclust::Mclust(data=x, G=2, modelNames="E", verbose=FALSE)
z <- c(tMu=y$parameters$mean,
tVar=y$parameters$variance$sigmasq,
tPi1=y$parameters$pro[1])
return(z)})
params <- rbind(paramsNormal, paramsTumor)
rownames(params)[c(1:2, 5:6)] = c("nMu1","nMu2","tMu1", "tMu2")
deltaMu2 <- params["tMu2",]-params["nMu2",]
deltaMu1 <- params["tMu1",]-params["nMu1",]
boundaryTumor <- (params["tMu2",] + params["tMu1",]) / 2
siCalc <- function(vectNml, boundTumor){ #calculate the sel. idx
x <- sum(boundTumor > vectNml) / length(vectNml)
return(x)
}
SI <- sapply(seq_len(dfNml),
function(i) siCalc(dfNml[,i], boundaryTumor[i]))
params <- rbind(params, deltaMu2, deltaMu1, SI)
mmParamsDf <- data.frame(t(params))
score <- NA
mmParamsDf$score <- mmParamsDf$SI*{
(mmParamsDf$deltaMu2-mmParamsDf$deltaMu1)-
(mmParamsDf$nVar+mmParamsDf$tVar)}
mmParamsDfS <- mmParamsDf[with(mmParamsDf, order(-score)), ]
return(mmParamsDfS)
}
mmParams <- mixModelParams(dfNml, dfTumor)
library(mclust)
mixModelParams <- function(dfNml, dfTumor) {
if(class(dfNml) == "SummarizedExperiment"){
dfNml <- assay(dfNml)
}
if(class(dfTumor) == "SummarizedExperiment"){
dfTumor <- assay(dfTumor)
}
paramsNormal <- apply(dfNml, 1, function(x) {
y <- mclust::Mclust(data=x, G=2, modelNames="E", verbose=FALSE)
z <- c(nMu=y$parameters$mean,
nVar=y$parameters$variance$sigmasq,
nPi1=y$parameters$pro[1])
return(z)})
paramsTumor <- apply(dfTumor, 1, function(x) {
y <- mclust::Mclust(data=x, G=2, modelNames="E", verbose=FALSE)
z <- c(tMu=y$parameters$mean,
tVar=y$parameters$variance$sigmasq,
tPi1=y$parameters$pro[1])
return(z)})
params <- rbind(paramsNormal, paramsTumor)
rownames(params)[c(1:2, 5:6)] = c("nMu1","nMu2","tMu1", "tMu2")
deltaMu2 <- params["tMu2",]-params["nMu2",]
deltaMu1 <- params["tMu1",]-params["nMu1",]
boundaryTumor <- (params["tMu2",] + params["tMu1",]) / 2
siCalc <- function(vectNml, boundTumor){ #calculate the sel. idx
x <- sum(boundTumor > vectNml) / length(vectNml)
return(x)
}
SI <- sapply(seq_len(dfNml),
function(i) siCalc(dfNml[,i], boundaryTumor[i]))
params <- rbind(params, deltaMu2, deltaMu1, SI)
mmParamsDf <- data.frame(t(params))
score <- NA
mmParamsDf$score <- mmParamsDf$SI*{
(mmParamsDf$deltaMu2-mmParamsDf$deltaMu1)-
(mmParamsDf$nVar+mmParamsDf$tVar)}
mmParamsDfS <- mmParamsDf[with(mmParamsDf, order(-score)), ]
return(mmParamsDfS)
}
mmParams <- mixModelParams(dfNml, dfTumor)
paramsNormal <- apply(dfNml, 1, function(x) {
y <- mclust::Mclust(data=x, G=2, modelNames="E", verbose=FALSE)
z <- c(nMu=y$parameters$mean,
nVar=y$parameters$variance$sigmasq,
nPi1=y$parameters$pro[1])
return(z)})
paramsTumor <- apply(dfTumor, 1, function(x) {
y <- mclust::Mclust(data=x, G=2, modelNames="E", verbose=FALSE)
z <- c(tMu=y$parameters$mean,
tVar=y$parameters$variance$sigmasq,
tPi1=y$parameters$pro[1])
return(z)})
paramsTumor
params <- rbind(paramsNormal, paramsTumor)
dim(paramsNormal)
dim(paramsTumor)
dfNml <- as.data.frame(matrix(data=rgamma(n=150, shape=2, rate=2),
nrow=10, ncol=15))
colnames(dfNml) <- paste0("patientN", seq_len(ncol(dfNml)))
rownames(dfNml) <- paste0("gene", seq_len(nrow(dfNml)))
dfTumor <- as.data.frame(matrix(data=rgamma(n=150, shape=4, rate=3),
nrow=10, ncol=15))
colnames(dfTumor) <- paste0("patientT", seq_len(ncol(dfTumor)))
rownames(dfTumor) <- paste0("gene", seq_len(nrow(dfTumor)))
mmParams <- mixModelParams(dfNml, dfTumor)
mixModelParams <- function(dfNml, dfTumor) {
if(class(dfNml) == "SummarizedExperiment"){
dfNml <- assay(dfNml)
}
if(class(dfTumor) == "SummarizedExperiment"){
dfTumor <- assay(dfTumor)
}
paramsNormal <- apply(dfNml, 1, function(x) {
y <- mclust::Mclust(data=x, G=2, modelNames="E", verbose=FALSE)
z <- c(nMu=y$parameters$mean,
nVar=y$parameters$variance$sigmasq,
nPi1=y$parameters$pro[1])
return(z)})
paramsTumor <- apply(dfTumor, 1, function(x) {
y <- mclust::Mclust(data=x, G=2, modelNames="E", verbose=FALSE)
z <- c(tMu=y$parameters$mean,
tVar=y$parameters$variance$sigmasq,
tPi1=y$parameters$pro[1])
return(z)})
params <- rbind(paramsNormal, paramsTumor)
rownames(params)[c(1:2, 5:6)] = c("nMu1","nMu2","tMu1", "tMu2")
deltaMu2 <- params["tMu2",]-params["nMu2",]
deltaMu1 <- params["tMu1",]-params["nMu1",]
boundaryTumor <- (params["tMu2",] + params["tMu1",]) / 2
siCalc <- function(vectNml, boundTumor){ #calculate the sel. idx
x <- sum(boundTumor > vectNml) / length(vectNml)
return(x)
}
SI <- sapply(seq_len(nrow(dfNml)),
function(i) siCalc(dfNml[,i], boundaryTumor[i]))
params <- rbind(params, deltaMu2, deltaMu1, SI)
mmParamsDf <- data.frame(t(params))
score <- NA
mmParamsDf$score <- mmParamsDf$SI*{
(mmParamsDf$deltaMu2-mmParamsDf$deltaMu1)-
(mmParamsDf$nVar+mmParamsDf$tVar)}
mmParamsDfS <- mmParamsDf[with(mmParamsDf, order(-score)), ]
return(mmParamsDfS)
}
mmParams <- mixModelParams(dfNml, dfTumor)
mmParams
plotGeneHist <- function(mmParams, dfNml, dfTumor, isof){
tidyDf <- as.data.frame(cbind(as.numeric(c(dfTumor[isof,], dfNml[isof,])),
as.factor(c(rep("tumor",ncol(dfTumor)), rep("normal",ncol(dfNml))))),
stringsAsFactors=FALSE)
colnames(tidyDf) <- c("expr", "type")
expr <- type <- ..density.. <- NULL # Setting the variables to NULL first
p1 <- ggplot(tidyDf, aes(x=expr, color=as.factor(type),
fill=as.factor(type),
group=as.factor(type))) +
theme_classic() +
geom_histogram(data=subset(tidyDf,type == 1),fill="#F8766D",
alpha=0.2, aes(y=..density..), binwidth = .2) +
geom_histogram(data=subset(tidyDf,type == 2),fill="#00BFC4",
alpha=0.2, aes(y=..density..), binwidth = .2) +
geom_rug(alpha=0.3, show.legend=FALSE) +
theme(axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position="none",
plot.title=element_text(size=12),
axis.text=element_text(size=8),
axis.title=element_text(size=8)) +
ggtitle(paste0(isof, " : SI = ",
round(mmParams[isof,"SI"],4))) +
xlab(expression(Log[2] *"(TPM Reads)")) +
stat_function(fun="dnorm", colour="#F8766D",
args=list(mean=mmParams[isof,"nMu2"],
sd=sqrt(mmParams[isof,"nVar"]))) +
stat_function(fun="dnorm", colour="#F8766D",
args=list(mean=mmParams[isof,"nMu1"],
sd=sqrt(mmParams[isof,"nVar"]))) +
stat_function(fun="dnorm", colour="#00BFC4",
args=list(mean=mmParams[isof,"tMu2"],
sd=sqrt(mmParams[isof,"tVar"]))) +
stat_function(fun="dnorm", colour="#00BFC4",
args=list(mean=mmParams[isof,"tMu1"],
sd=sqrt(mmParams[isof,"tVar"])))
print(p1)
}
library(ggplot2)
library(ggplot2)
library(ggplot2)
install.packages("ggplot2")
library(ggplot2)
plotGeneHist <- function(mmParams, dfNml, dfTumor, isof){
tidyDf <- as.data.frame(cbind(as.numeric(c(dfTumor[isof,], dfNml[isof,])),
as.factor(c(rep("tumor",ncol(dfTumor)), rep("normal",ncol(dfNml))))),
stringsAsFactors=FALSE)
colnames(tidyDf) <- c("expr", "type")
expr <- type <- ..density.. <- NULL # Setting the variables to NULL first
p1 <- ggplot(tidyDf, aes(x=expr, color=as.factor(type),
fill=as.factor(type),
group=as.factor(type))) +
theme_classic() +
geom_histogram(data=subset(tidyDf,type == 1),fill="#F8766D",
alpha=0.2, aes(y=..density..), binwidth = .2) +
geom_histogram(data=subset(tidyDf,type == 2),fill="#00BFC4",
alpha=0.2, aes(y=..density..), binwidth = .2) +
geom_rug(alpha=0.3, show.legend=FALSE) +
theme(axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position="none",
plot.title=element_text(size=12),
axis.text=element_text(size=8),
axis.title=element_text(size=8)) +
ggtitle(paste0(isof, " : SI = ",
round(mmParams[isof,"SI"],4))) +
xlab(expression(Log[2] *"(TPM Reads)")) +
stat_function(fun="dnorm", colour="#F8766D",
args=list(mean=mmParams[isof,"nMu2"],
sd=sqrt(mmParams[isof,"nVar"]))) +
stat_function(fun="dnorm", colour="#F8766D",
args=list(mean=mmParams[isof,"nMu1"],
sd=sqrt(mmParams[isof,"nVar"]))) +
stat_function(fun="dnorm", colour="#00BFC4",
args=list(mean=mmParams[isof,"tMu2"],
sd=sqrt(mmParams[isof,"tVar"]))) +
stat_function(fun="dnorm", colour="#00BFC4",
args=list(mean=mmParams[isof,"tMu1"],
sd=sqrt(mmParams[isof,"tVar"])))
print(p1)
}
dfNml <- as.data.frame(matrix(data=rgamma(n=150, shape=2, rate=2),
nrow=10, ncol=15))
colnames(dfNml) <- paste0("patientN", seq_len(ncol(dfNml)))
rownames(dfNml) <- paste0("gene", seq_len(nrow(dfNml)))
dfTumor <- as.data.frame(matrix(data=rgamma(n=150, shape=4, rate=3),
nrow=10, ncol=15))
colnames(dfTumor) <- paste0("patientT", seq_len(ncol(dfTumor)))
rownames(dfTumor) <- paste0("gene", seq_len(nrow(dfTumor)))
mmParams <- mixModelParams(dfNml, dfTumor)
isof <- rownames(mmParams)[1]
plotGeneHist(mmParams, dfNml, dfTumor, isof)
mixModelParams <- function(dfNml, dfTumor) {
if(class(dfNml) == "SummarizedExperiment"){
dfNml <- assay(dfNml)
}
if(class(dfTumor) == "SummarizedExperiment"){
dfTumor <- assay(dfTumor)
}
paramsNormal <- apply(dfNml, 1, function(x) {
y <- mclust::Mclust(data=x, G=2, modelNames="E", verbose=FALSE)
z <- c(nMu=y$parameters$mean,
nVar=y$parameters$variance$sigmasq,
nPi1=y$parameters$pro[1])
return(z)})
paramsTumor <- apply(dfTumor, 1, function(x) {
y <- mclust::Mclust(data=x, G=2, modelNames="E", verbose=FALSE)
z <- c(tMu=y$parameters$mean,
tVar=y$parameters$variance$sigmasq,
tPi1=y$parameters$pro[1])
return(z)})
params <- rbind(paramsNormal, paramsTumor)
rownames(params)[c(1:2, 5:6)] = c("nMu1","nMu2","tMu1", "tMu2")
deltaMu2 <- params["tMu2",]-params["nMu2",]
deltaMu1 <- params["tMu1",]-params["nMu1",]
boundaryTumor <- (params["tMu2",] + params["tMu1",]) / 2
siCalc <- function(vectNml, boundTumor){ #calculate the sel. idx
x <- sum(boundTumor > vectNml) / length(vectNml)
return(x)
}
SI <- sapply(seq_len(nrow(dfNml)),
function(i) siCalc(dfNml[,i], boundaryTumor[i]))
params <- rbind(params, deltaMu2, deltaMu1, SI)
mmParamsDf <- data.frame(t(params))
score <- NA
mmParamsDf$score <- mmParamsDf$SI*{
(mmParamsDf$deltaMu2-mmParamsDf$deltaMu1)-
(mmParamsDf$nVar+mmParamsDf$tVar)}
mmParamsDfS <- mmParamsDf[with(mmParamsDf, order(-score)), ]
return(mmParamsDfS)
}
mmParams <- mixModelParams(dfNml, dfTumor)
isof <- rownames(mmParams)[1]
plotGeneHist(mmParams, dfNml, dfTumor, isof)
library(mclust)
mmParams <- mixModelParams(dfNml, dfTumor)
isof <- rownames(mmParams)[1]
plotGeneHist(mmParams, dfNml, dfTumor, isof)
library(ggplot2)
remove.packages(c("ggplot2", "data.table"))
install.packages('Rcpp', dependencies = TRUE)
install.packages('ggplot2', dependencies = TRUE)
library(ggplot2)