Permutation p-values
We use Myc and Pten as examples to show how permutation p-values are calculated.
library("ggplot2")
load("DataSummary/Pooled_MA_Zscore.RData")
attach(Pooled_MA_Zscore)
source("/sibcb2/bioinformatics/Script/PathwayEnrichment.R")
gene = c("Myc", "Pten")
control_sgRNA = rownames(annT)[annT$group == "control"]
control_list = list()
for(i in 1:5000)
control_list[[i]] = sample(control_sgRNA, 6)
Log2FC = sort(FcM[, "Vitro_vs_Plasmid"])
randomT = data.frame(controlZS = SimpleRankTest(Log2FC, control_list))
test_zscore = zM[gene, "Vitro_vs_Plasmid"]
test_pvalue = pM[gene, "Vitro_vs_Plasmid"]
p <- ggplot(randomT, aes(controlZS) )
p <- p + theme_bw() + labs(x = "Z-score for control sgRNAs", y = "Density", title = "")
p <- p + geom_line(stat="density", size = 0.6) + theme(legend.title=element_blank()) + xlim(-5, 5)
p <- p + geom_vline(xintercept = test_zscore, size = 1, linetype = "longdash", color = c("blue", "red"))
p <- p + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + geom_text(x = -3, y = 0.3, label = paste("Myc\np =", test_pvalue[1]), color = "blue")
p <- p + geom_text(x = 4, y = 0.3, label = paste("Pten\np =", test_pvalue[2]), color = "red")
print(p)
sortRank = 1:length(Log2FC)
test_zscore = signif(test_zscore, 4)
MycRank = sortRank[names(Log2FC) %in% rownames(annT)[annT$Plate == "Myc"]]
PtenRank = sortRank[names(Log2FC) %in% rownames(annT)[annT$Plate == "Pten"]]
p <- ggplot(data = data.frame(sortRank), aes(x = sortRank))
p <- p + theme_bw() + labs(x = "", y = "", title = paste("Myc, Z-score =", test_zscore[1]))
p <- p + xlim(1, max(sortRank))
p <- p + geom_vline(data = data.frame(MycRank), aes(xintercept = MycRank), size = 1.5, color = "blue")
p <- p + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + theme(axis.ticks = element_blank(), axis.text.x = element_blank())
print(p)
p <- ggplot(data = data.frame(sortRank), aes(x = sortRank))
p <- p + theme_bw() + labs(x = "", y = "", title = paste("Pten, Z-score =", test_zscore[2]))
p <- p + xlim(1, max(sortRank))
p <- p + geom_vline(data = data.frame(PtenRank), aes(xintercept = PtenRank), size = 1.5, color = "red")
p <- p + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + theme(axis.ticks = element_blank(), axis.text.x = element_blank())
print(p)
sessionInfo
R version 4.1.0 (2021-05-18)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server 7.6 (Maipo)
Matrix products: default
BLAS/LAPACK: /sibcb2/bioinformatics/software/Miniconda3/lib/libopenblasp-r0.3.15.so
locale:
[1] C
attached base packages:
[1] grid stats graphics grDevices utils datasets
[7] methods base
other attached packages:
[1] stringr_1.5.0 dplyr_1.1.2 scales_1.2.0
[4] RColorBrewer_1.1-3 pheatmap_1.0.12 VennDiagram_1.6.20
[7] futile.logger_1.4.3 ggrepel_0.9.1 ggplot2_3.4.2
[10] rmarkdown_2.20
loaded via a namespace (and not attached):
[1] Rcpp_1.0.9 highr_0.9 jquerylib_0.1.4
[4] bslib_0.2.5.1 pillar_1.9.0 compiler_4.1.0
[7] formatR_1.11 futile.options_1.0.1 tools_4.1.0
[10] downlit_0.4.2 digest_0.6.29 jsonlite_1.8.7
[13] memoise_2.0.1 evaluate_0.20 lifecycle_1.0.3
[16] tibble_3.2.1 gtable_0.3.0 pkgconfig_2.0.3
[19] rlang_1.1.1 rstudioapi_0.15.0 cli_3.6.1
[22] distill_1.5 yaml_2.3.5 xfun_0.37
[25] fastmap_1.1.0 withr_2.5.0 knitr_1.42
[28] sass_0.4.0 generics_0.1.3 vctrs_0.6.3
[31] tidyselect_1.2.1 glue_1.6.2 R6_2.5.1
[34] fansi_1.0.3 bookdown_0.33 farver_2.1.1
[37] lambda.r_1.2.4 magrittr_2.0.3 htmltools_0.5.2
[40] colorspace_2.0-3 labeling_0.4.2 utf8_1.2.2
[43] stringi_1.7.8 munsell_0.5.0 cachem_1.0.6