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
## # ℹ 394,811 more rows

You 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,
        cores = 1, 
        
        # This is ONLY for speeding up the Vignette execution
        draws_after_tail = 1
    )

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 × 13]>                  0                        0

The 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.4.0 beta (2024-04-15 r86425)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [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       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ppcseq_1.12.0      rstan_2.32.6       StanHeaders_2.32.7 dplyr_1.1.4       
## [5] knitr_1.46         BiocStyle_2.32.0  
## 
## loaded via a namespace (and not attached):
##  [1] benchmarkmeData_1.0.4 gtable_0.3.5          tensorA_0.36.2.1     
##  [4] xfun_0.43             bslib_0.7.0           ggplot2_3.5.1        
##  [7] QuickJSR_1.1.3        inline_0.3.19         lattice_0.22-6       
## [10] vctrs_0.6.5           tools_4.4.0           generics_0.1.3       
## [13] stats4_4.4.0          curl_5.2.1            parallel_4.4.0       
## [16] tibble_3.2.1          fansi_1.0.6           highr_0.10           
## [19] pkgconfig_2.0.3       Matrix_1.7-0          checkmate_2.3.1      
## [22] distributional_0.4.0  RcppParallel_5.1.7    lifecycle_1.0.4      
## [25] farver_2.1.1          stringr_1.5.1         compiler_4.4.0       
## [28] tinytex_0.50          statmod_1.5.0         munsell_0.5.1        
## [31] codetools_0.2-20      htmltools_0.5.8.1     sass_0.4.9           
## [34] yaml_2.3.8            pillar_1.9.0          jquerylib_0.1.4      
## [37] tidyr_1.3.1           arrayhelpers_1.1-0    cachem_1.0.8         
## [40] limma_3.60.0          magick_2.8.3          iterators_1.0.14     
## [43] foreach_1.5.2         abind_1.4-5           posterior_1.5.0      
## [46] tidybayes_3.0.6       tidyselect_1.2.1      locfit_1.5-9.9       
## [49] digest_0.6.35         svUnit_1.0.6          stringi_1.8.3        
## [52] purrr_1.0.2           bookdown_0.39         labeling_0.4.3       
## [55] fastmap_1.1.1         grid_4.4.0            colorspace_2.1-0     
## [58] cli_3.6.2             magrittr_2.0.3        loo_2.7.0            
## [61] pkgbuild_1.4.4        utf8_1.2.4            withr_3.0.0          
## [64] edgeR_4.2.0           scales_1.3.0          backports_1.4.1      
## [67] httr_1.4.7            rmarkdown_2.26        matrixStats_1.3.0    
## [70] benchmarkme_1.0.8     gridExtra_2.3         coda_0.19-4.1        
## [73] evaluate_0.23         doParallel_1.0.17     V8_4.4.2             
## [76] ggdist_3.3.2          rstantools_2.4.0      rlang_1.1.3          
## [79] Rcpp_1.0.12           glue_1.7.0            BiocManager_1.30.22  
## [82] jsonlite_1.8.8        R6_2.5.1

Probabilistic Outlier Identification for RNA Sequencing Generalized Linear Models

2024-05-01

Package

1 Introduction

2 Installation and use