library("IHW")
library("ggplot2")
pfin <- read.table('pfins.tsv', sep='\t', header=TRUE, row.names = 1)
ihwRes <- ihw(p ~ var, data = pfin, alpha = 0.1)
pfin$padj <- adj_pvalues(ihwRes)
pfin$BH <- p.adjust(pfin$p, method = "BH")
head(pfin)
## p var padj BH
## 1 0.08518485 0.04687183 0.7391697 0.7391697
## 2 0.29819847 0.03476625 0.8619217 0.8619217
## 3 0.21340210 0.03651571 0.8320083 0.8320083
## 4 0.52053214 0.11713521 0.9242342 0.9242342
## 5 0.50901341 0.34885276 0.9219028 0.9219028
## 6 0.43508199 0.09018223 0.9014622 0.9014622
plot(ihwRes)
ggplot(aes(x = rank(var), y = -log10(p)), data=pfin) + geom_hex(bins = 100) + ylab(expression(-log[10]~p))
ggplot(pfin, aes(x = p)) + geom_histogram(binwidth = 0.005, boundary = 0)
There seems to be a small enrichment of small p-values; we can also verify this using the Schweder-Spjøtvoll plot, see the orange line.
metap::schweder(pfin$p, drawline = c("ls", "bh"),
ls.col = "red", ls.lty = 1, ls.control = list(frac = 1/3),
bh.col = "blue")
Also the BH method indicates that there are a few discoveries to be made, but only with FDR of 40% or higher.
ggplot(dplyr::filter(pfin, BH < 1), aes(x = BH)) +
geom_histogram(binwidth = 0.01, boundary = 0)
ggplot(pfin, aes(x = var)) + geom_histogram(binwidth = 0.01)
ggplot(pfin, aes(x = p)) + geom_histogram(binwidth = 0.040, boundary = 0) +
facet_wrap(~ groups_by_filter(pfin$var, 8), nrow = 2)
There seems to be a weak trend that the enrichment of small p-values is strongest for small var. Perhaps this could also be a normalization / data preprocessing artefact?
var is not a strong covariate. You’ll have to search for others.
devtools::session_info()
## ─ Session info ────────────────────────────────────────────────────────────
## setting value
## version R version 3.5.1 (2018-07-02)
## os macOS 10.14.1
## system x86_64, darwin15.6.0
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz Europe/Berlin
## date 2018-11-18
##
## ─ Packages ────────────────────────────────────────────────────────────────
## package * version date lib source
## assertthat 0.2.0 2017-04-11 [1] CRAN (R 3.5.0)
## backports 1.1.2 2017-12-13 [1] CRAN (R 3.5.0)
## base64enc 0.1-3 2015-07-28 [1] CRAN (R 3.5.0)
## bibtex 0.4.2 2017-06-30 [1] CRAN (R 3.5.0)
## bindr 0.1.1 2018-03-13 [1] CRAN (R 3.5.0)
## bindrcpp * 0.2.2 2018-03-29 [1] CRAN (R 3.5.0)
## BiocGenerics 0.28.0 2018-10-30 [1] Bioconductor
## BiocManager 1.30.3 2018-10-10 [1] CRAN (R 3.5.0)
## BiocStyle * 2.10.0 2018-10-30 [1] Bioconductor
## bookdown 0.7 2018-02-18 [1] CRAN (R 3.5.0)
## callr 3.0.0 2018-08-24 [1] CRAN (R 3.5.0)
## cli 1.0.1 2018-09-25 [1] CRAN (R 3.5.0)
## codetools 0.2-15 2016-10-05 [1] CRAN (R 3.5.0)
## colorspace 1.3-2 2016-12-14 [1] CRAN (R 3.5.0)
## crayon 1.3.4 2017-09-16 [1] CRAN (R 3.5.0)
## desc 1.2.0 2018-05-01 [1] CRAN (R 3.5.0)
## devtools 2.0.1 2018-10-26 [1] CRAN (R 3.5.1)
## digest 0.6.18 2018-10-10 [1] CRAN (R 3.5.0)
## dplyr 0.7.8 2018-11-10 [1] CRAN (R 3.5.0)
## evaluate 0.12 2018-10-09 [1] CRAN (R 3.5.0)
## fdrtool 1.2.15 2015-07-08 [1] CRAN (R 3.5.0)
## fs 1.2.6 2018-08-23 [1] CRAN (R 3.5.0)
## gbRd 0.4-11 2012-10-01 [1] CRAN (R 3.5.0)
## ggplot2 * 3.1.0 2018-10-25 [1] CRAN (R 3.5.0)
## glue 1.3.0 2018-07-17 [1] CRAN (R 3.5.0)
## gtable 0.2.0 2016-02-26 [1] CRAN (R 3.5.0)
## hexbin * 1.27.2 2018-01-15 [1] CRAN (R 3.5.0)
## htmltools 0.3.6 2017-04-28 [1] CRAN (R 3.5.0)
## IHW * 1.10.0 2018-10-30 [1] Bioconductor
## knitr 1.20 2018-02-20 [1] CRAN (R 3.5.0)
## labeling 0.3 2014-08-23 [1] CRAN (R 3.5.0)
## lattice 0.20-38 2018-11-04 [1] CRAN (R 3.5.0)
## lazyeval 0.2.1 2017-10-29 [1] CRAN (R 3.5.0)
## lpsymphony 1.10.0 2018-10-30 [1] Bioconductor (R 3.5.1)
## magrittr 1.5 2014-11-22 [1] CRAN (R 3.5.0)
## memoise 1.1.0 2017-04-21 [1] CRAN (R 3.5.0)
## metap 1.0 2018-07-25 [1] CRAN (R 3.5.0)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 3.5.0)
## pillar 1.3.0 2018-07-14 [1] CRAN (R 3.5.0)
## pkgbuild 1.0.2 2018-10-16 [1] CRAN (R 3.5.0)
## pkgconfig 2.0.2 2018-08-16 [1] CRAN (R 3.5.0)
## pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.5.0)
## plyr 1.8.4 2016-06-08 [1] CRAN (R 3.5.0)
## prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.5.0)
## processx 3.2.0 2018-08-16 [1] CRAN (R 3.5.0)
## ps 1.2.1 2018-11-06 [1] CRAN (R 3.5.0)
## purrr 0.2.5 2018-05-29 [1] CRAN (R 3.5.0)
## R6 2.3.0 2018-10-04 [1] CRAN (R 3.5.0)
## Rcpp 1.0.0 2018-11-07 [1] CRAN (R 3.5.0)
## Rdpack 0.10-1 2018-10-04 [1] CRAN (R 3.5.0)
## remotes 2.0.2 2018-10-30 [1] CRAN (R 3.5.0)
## rlang 0.3.0.1 2018-10-25 [1] CRAN (R 3.5.0)
## rmarkdown 1.10 2018-06-11 [1] CRAN (R 3.5.0)
## rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.5.0)
## scales 1.0.0 2018-08-09 [1] CRAN (R 3.5.1)
## sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.5.0)
## slam 0.1-43 2018-04-23 [1] CRAN (R 3.5.0)
## stringi 1.2.4 2018-07-20 [1] CRAN (R 3.5.0)
## stringr 1.3.1 2018-05-10 [1] CRAN (R 3.5.0)
## testthat 2.0.1 2018-10-13 [1] CRAN (R 3.5.1)
## tibble 1.4.2 2018-01-22 [1] CRAN (R 3.5.0)
## tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.5.0)
## usethis 1.4.0 2018-08-14 [1] CRAN (R 3.5.0)
## withr 2.1.2 2018-03-15 [1] CRAN (R 3.5.0)
## xfun 0.4 2018-10-23 [1] CRAN (R 3.5.0)
## yaml 2.2.0 2018-07-25 [1] CRAN (R 3.5.0)
##
## [1] /Library/Frameworks/R.framework/Versions/3.5/Resources/library