Contents

1 Installation

The package can be installed from bioconductor

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("scDDboost")

Issue can be reported at “https://github.com/wiscstatman/scDDboost/issues”

2 Introduction

scDDboost scores evidence of a gene being differentially distributed(DD) across two conditions for single cell RNA-seq data. Higher resolution brings several chanllenges for analyzing the data, specifically, the distribution of gene expression tends to have high prevalence of zero and multi-modes. To account for those characteristics and utilizing some biological intuition, we view the expression values sampled from a pool of cells mixed by distinct cellular subtypes blind to condition label. Consequently, the distributional change can be fully determined by the the change of subtype proportions. One tricky part is that not any change of proportions will lead to a distributional change. Given that some genes could be equivalent expressed across several subtypes, even the individual subytpe proportion may differ between conditions but as long as the aggregated proportions over those subtypes remain the same between conditions, it will not introduce different distribution. For example

Proportions of subtypes 1 and 2 changed between the 2 conditions. The gene is not DD if subtype 1 and 2 have the same expression level

For subtype 1 and 2 have different expression level, there is different distribution

3 Posterior probability of a gene being DD

pdd is the core function developed to quantify the posterior probabilities of DD for input genes. Let’s look at an example,

suppressMessages(library(scDDboost))

Next, we load the toy simulated example a object that we will use for identifying and classifying DD genes.

data(sim_dat)

Verify that this object is a member of the SingleCellExperiment class and that it contains 200 cells and 1000 genes. The colData slot (which contains a dataframe of metadata for the cells) should have a column that contains the biological condition or grouping of interest. In this example data, that variable is the condition variable. Note that the input gene set needs to be a matrix of normalized counts. We run the function pdd

data_counts <- SummarizedExperiment::assays(sim_dat)$counts
conditions <- SummarizedExperiment::colData(sim_dat)$conditions
rownames(data_counts) <- seq_len(1000)

##here we use 2 cores to compute the distance matrix
bp <- BiocParallel::MulticoreParam(2)
D_c <- calD(data_counts,bp)

ProbDD <- pdd(data = data_counts,cd = conditions, bp = bp, D = D_c)

There are 4 input parameters needed to be specified by user, the dataset, the condition label, number of cpu cores used for computation and a distance matrix of cells. Other input parameters have default settings.

3.1 clustering of cells

We provide a default method of getting the distance matrix, archived by calD, in general pdd accept all valid distance matrix. User can also input a cluster label rather than distance matrix for the argument D, but the random distancing mechanism which relies on distance matrix will be disabled and random should be set to false.

For the number of sutypes, we provide a default function detK, which consider the smallest number of sutypes such that the ratio of difference within cluster between difference between clusters become smaller than a threshold (default setting is 1).

If user have other ways to determine \(K\), \(K\) should be specified in pdd.

## determine the number of subtypes
K <- detK(D_c)

If we set threshold to be 5% then we have estimated DD genes

EDD <- which(ProbDD > 0.95)

Notice that, pdd is actually local false discovery rate, this is a conservative estimation of DD genes. We could gain further power, let index gene by \(g = 1,2,...,G\) and let \(p_g = P(DD_g | \text{data})\), \(p_{(1)},...,p_{(G)}\) be ranked local false discovery rate from small to large. To control the false discovery rate at 5%, our positive set is those genes with the \(s^*\) smallest lFDR, where \[s^* = \text{argmax}_s\{s,\frac{\Sigma_{i = 1}^s p_{(i)}}{s} \leq 0.05\}\]

EDD <- getDD(ProbDD,0.05)

Function getDD extracts the estimated DD genes using the above transformation.

Session Information

sessionInfo()
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.16-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] scDDboost_1.0.0  ggplot2_3.3.6    BiocStyle_2.26.0
## 
## loaded via a namespace (and not attached):
##  [1] MatrixGenerics_1.10.0       Biobase_2.58.0             
##  [3] pkgload_1.3.1               sass_0.4.2                 
##  [5] jsonlite_1.8.3              gtools_3.9.3               
##  [7] brio_1.1.3                  bslib_0.4.0                
##  [9] assertthat_0.2.1            BiocManager_1.30.19        
## [11] highr_0.9                   stats4_4.2.1               
## [13] EBSeq_1.38.0                GenomeInfoDbData_1.2.9     
## [15] yaml_2.3.6                  pillar_1.8.1               
## [17] lattice_0.20-45             glue_1.6.2                 
## [19] RcppEigen_0.3.3.9.2         digest_0.6.30              
## [21] GenomicRanges_1.50.0        XVector_0.38.0             
## [23] colorspace_2.0-3            htmltools_0.5.3            
## [25] Matrix_1.5-1                pkgconfig_2.0.3            
## [27] magick_2.7.3                bookdown_0.29              
## [29] zlibbioc_1.44.0             scales_1.2.1               
## [31] BiocParallel_1.32.0         tibble_3.1.8               
## [33] generics_0.1.3              farver_2.1.1               
## [35] IRanges_2.32.0              cachem_1.0.6               
## [37] withr_2.5.0                 SummarizedExperiment_1.28.0
## [39] blockmodeling_1.1.3         BiocGenerics_0.44.0        
## [41] cli_3.4.1                   magrittr_2.0.3             
## [43] mclust_6.0.0                evaluate_0.17              
## [45] fansi_1.0.3                 Oscope_1.28.0              
## [47] gplots_3.1.3                tools_4.2.1                
## [49] lifecycle_1.0.3             matrixStats_0.62.0         
## [51] stringr_1.4.1               S4Vectors_0.36.0           
## [53] munsell_0.5.0               cluster_2.1.4              
## [55] DelayedArray_0.24.0         compiler_4.2.1             
## [57] jquerylib_0.1.4             GenomeInfoDb_1.34.0        
## [59] caTools_1.18.2              rlang_1.0.6                
## [61] grid_4.2.1                  RCurl_1.98-1.9             
## [63] SingleCellExperiment_1.20.0 bitops_1.0-7               
## [65] labeling_0.4.2              rmarkdown_2.17             
## [67] testthat_3.1.5              gtable_0.3.1               
## [69] codetools_0.2-18            DBI_1.1.3                  
## [71] R6_2.5.1                    knitr_1.40                 
## [73] dplyr_1.0.10                fastmap_1.1.0              
## [75] utf8_1.2.2                  rprojroot_2.0.3            
## [77] desc_1.4.2                  KernSmooth_2.23-20         
## [79] stringi_1.7.8               parallel_4.2.1             
## [81] Rcpp_1.0.9                  vctrs_0.5.0                
## [83] tidyselect_1.2.0            xfun_0.34