Probabilistic Outlier Identification for RNA Sequencing Generalized Linear Models
2021-10-26
Package
ppcseq 1.2.0
1 Introduction
Relative transcript abundance has proven to be a valuable tool for understanding the function of genes in biological systems. For the differential analysis of transcript abundance using RNA sequencing data, the negative binomial model is by far the most frequently adopted. However, common methods that are based on a negative binomial model are not robust to extreme outliers, which we found to be abundant in public datasets. So far, no rigorous and probabilistic methods for detection of outliers have been developed for RNA sequencing data, leaving the identification mostly to visual inspection. Recent advances in Bayesian computation allow large-scale comparison of observed data against its theoretical distribution given in a statistical model. Here we propose ppcseq, a key quality-control tool for identifying transcripts that include outlier data points in differential expression analysis, which do not follow a negative binomial distribution. Applying ppcseq to analyse several publicly available datasets using popular tools, we show that from 3 to 10 percent of differentially abundant transcripts across algorithms and datasets had statistics inflated by the presence of outliers.
2 Installation and use
The input data set is a tidy representation of a differential gene transcript abundance analysis
library(dplyr)
library(ppcseq)To install:
Before install, for linux systems, in order to exploit multi-threading, from R write:
fileConn<-file("~/.R/Makevars")
writeLines(c( "CXX14FLAGS += -O3","CXX14FLAGS += -DSTAN_THREADS", "CXX14FLAGS += -pthread"), fileConn)
close(fileConn)Multi-threading allows the sampling or variational bayes to share the computation on multiple cores.
Then, install with
if(!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ppcseq")You can get the test dataset with
data("counts")
counts ## # A tibble: 394,821 × 9
## sample symbol logCPM LR PValue FDR value W Label
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <chr>
## 1 10922PP SLC16A12 1.39 41.1 1.46e-10 0.00000274 160 -0.129 High
## 2 10935PP SLC16A12 1.39 41.1 1.46e-10 0.00000274 150 -0.127 High
## 3 10973PP SLC16A12 1.39 41.1 1.46e-10 0.00000274 146 -0.426 High
## 4 10976PP SLC16A12 1.39 41.1 1.46e-10 0.00000274 347 -0.0164 High
## 5 10985PP SLC16A12 1.39 41.1 1.46e-10 0.00000274 175 -0.135 High
## 6 11026PP SLC16A12 1.39 41.1 1.46e-10 0.00000274 244 0.125 High
## 7 11045PP SLC16A12 1.39 41.1 1.46e-10 0.00000274 399 -0.0892 High
## 8 11082PP SLC16A12 1.39 41.1 1.46e-10 0.00000274 100 0.261 Neoadjuvant
## 9 11086PP SLC16A12 1.39 41.1 1.46e-10 0.00000274 37 -0.132 Neoadjuvant
## 10 11103PP SLC16A12 1.39 41.1 1.46e-10 0.00000274 73 0.146 Neoadjuvant
## # … with 394,811 more rowsYou can identify anrtefactual calls from your differential transcribt anundance analysis, due to outliers.
# Import libraries
if(Sys.info()[['sysname']] == "Linux")
counts.ppc =
counts %>%
mutate(is_significant = FDR < 0.0001) %>%
identify_outliers(
formula = ~ Label,
.sample = sample,
.transcript = symbol,
.abundance = value,
.significance = PValue,
.do_check = is_significant,
percent_false_positive_genes = 5,
approximate_posterior_inference = FALSE
)The new posterior predictive check has been added to the original data frame
if(Sys.info()[['sysname']] == "Linux")
counts.ppc ## # A tibble: 1 × 4
## symbol `sample wise data` `ppc samples failed` `tot deleterious outliers`
## <chr> <list> <int> <int>
## 1 SLC16A12 <tibble [21 × 14]> 0 0The new data frame contains plots for each gene
We can visualise the top five differentially transcribed genes
if(Sys.info()[['sysname']] == "Linux")
counts.ppc_plots =
counts.ppc %>%
plot_credible_intervals() if(Sys.info()[['sysname']] == "Linux")
counts.ppc_plots %>%
pull(plot) %>%
.[seq_len(1)]## [[1]]
sessionInfo()## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.14-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.14-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ppcseq_1.2.0 dplyr_1.0.7 knitr_1.36 BiocStyle_2.22.0
##
## loaded via a namespace (and not attached):
## [1] matrixStats_0.61.0 doParallel_1.0.16 httr_1.4.2
## [4] rstan_2.21.2 tensorA_0.36.2 tools_4.1.1
## [7] backports_1.2.1 bslib_0.3.1 utf8_1.2.2
## [10] R6_2.5.1 DBI_1.1.1 colorspace_2.0-2
## [13] ggdist_3.0.0 withr_2.4.2 tidyselect_1.1.1
## [16] gridExtra_2.3 prettyunits_1.1.1 processx_3.5.2
## [19] curl_4.3.2 compiler_4.1.1 cli_3.0.1
## [22] arrayhelpers_1.1-0 labeling_0.4.2 bookdown_0.24
## [25] posterior_1.1.0 sass_0.4.0 scales_1.1.1
## [28] checkmate_2.0.0 callr_3.7.0 stringr_1.4.0
## [31] digest_0.6.28 StanHeaders_2.21.0-7 rmarkdown_2.11
## [34] benchmarkmeData_1.0.4 pkgconfig_2.0.3 htmltools_0.5.2
## [37] parallelly_1.28.1 fastmap_1.1.0 limma_3.50.0
## [40] highr_0.9 rlang_0.4.12 jquerylib_0.1.4
## [43] farver_2.1.0 generics_0.1.1 svUnit_1.0.6
## [46] jsonlite_1.7.2 distributional_0.2.2 inline_0.3.19
## [49] magrittr_2.0.1 loo_2.4.1 Matrix_1.3-4
## [52] Rcpp_1.0.7 munsell_0.5.0 fansi_0.5.0
## [55] abind_1.4-5 lifecycle_1.0.1 furrr_0.2.3
## [58] stringi_1.7.5 yaml_2.2.1 edgeR_3.36.0
## [61] pkgbuild_1.2.0 grid_4.1.1 parallel_4.1.1
## [64] listenv_0.8.0 crayon_1.4.1 lattice_0.20-45
## [67] locfit_1.5-9.4 magick_2.7.3 ps_1.6.0
## [70] pillar_1.6.4 codetools_0.2-18 stats4_4.1.1
## [73] rstantools_2.1.1 glue_1.4.2 evaluate_0.14
## [76] V8_3.4.2 BiocManager_1.30.16 RcppParallel_5.1.4
## [79] vctrs_0.3.8 foreach_1.5.1 gtable_0.3.0
## [82] purrr_0.3.4 tidyr_1.1.4 benchmarkme_1.0.7
## [85] future_1.22.1 assertthat_0.2.1 ggplot2_3.3.5
## [88] xfun_0.27 coda_0.19-4 tibble_3.1.5
## [91] iterators_1.0.13 tidybayes_3.0.1 globals_0.14.0
## [94] ellipsis_0.3.2