POWSC: power and sample size snalysis for single-cell RNA-seq

Installation and quick start

Installation

To install this package, start R (version > “4.1”) and enter:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("POWSC")

Quick start on using POWSC

We use one scRNA-seq template data about ES-MEF cell types (GSE29087). In this case, we adopt one cell type (fibro) for the two-group comparison corresponding to the 1st scenario.

suppressMessages(library(POWSC))

## Warning: multiple methods tables found for 'setequal'

## Warning: replacing previous import 'BiocGenerics::setequal' by
## 'S4Vectors::setequal' when loading 'IRanges'

## Warning: replacing previous import 'BiocGenerics::setequal' by
## 'S4Vectors::setequal' when loading 'GenomeInfoDb'

## Warning: replacing previous import 'BiocGenerics::setequal' by
## 'S4Vectors::setequal' when loading 'GenomicRanges'

## Warning: replacing previous import 'BiocGenerics::setequal' by
## 'S4Vectors::setequal' when loading 'XVector'

## Warning: replacing previous import 'BiocGenerics::setequal' by
## 'S4Vectors::setequal' when loading 'SummarizedExperiment'

## Warning: replacing previous import 'BiocGenerics::setequal' by
## 'S4Vectors::setequal' when loading 'S4Arrays'

## Warning: replacing previous import 'BiocGenerics::setequal' by
## 'S4Vectors::setequal' when loading 'DelayedArray'

## Warning: replacing previous import 'BiocGenerics::setequal' by
## 'S4Vectors::setequal' when loading 'SparseArray'

## Warning: replacing previous import 'S4Arrays::read_block' by
## 'DelayedArray::read_block' when loading 'SummarizedExperiment'

## Warning: replacing previous import 'BiocGenerics::setequal' by
## 'S4Vectors::setequal' when loading 'MAST'

data("es_mef_sce")
sce = es_mef_sce[, colData(es_mef_sce)$cellTypes == "fibro"]
set.seed(12)
rix = sample(1:nrow(sce), 500)
sce = sce[rix, ]
est_Paras = Est2Phase(sce)
sim_size = c(100, 200) # A numeric vector
pow_rslt = runPOWSC(sim_size = sim_size, est_Paras = est_Paras,per_DE=0.05, DE_Method = "MAST", Cell_Type = "PW")

Background

Determining the sample size for adequate power to detect statistical significance is a crucial step at the design stage for high-throughput experiments. Even though a number of methods and tools are available for sample size calculation for microarray and RNA-seq under the context of differential expression, this topic in the field of single-cell RNA sequencing is understudied. Moreover, the unique data characteristics present in scRNA-seq including sparsity and heterogeneity gain the challenge.

Introduction

POWSC is an R package designed for power assessment and sample size estimation in scRNA-seq experiment. It contains three main functionalities: (1). Parameter Estimation: adopted and modified from the core of SC2P. (2). Data Simulation: consider two cases: paired-wise comparison and multiple comparisons. (3). Power Evaluation: provide both stratified (detailed) and marginal powers.

For the model assumption, we adopt the mixture of zero-inflated poisson (ZIP) and log-normal poisson (LNP).

The power assessor in POWSC provides a variety of power evaluations including stratified and marginal power analyses for DEs characterized by two forms (phase transition or magnitude tuning), under different comparison scenarios. In addition, POWSC offers information for optimizing the tradeoffs between sample size and sequencing depth with the same total reads.

Use POWSC

In the context of differential expression (DE) analysis, scientists are usually interested in two different scenarios: (1) within cell type: comparing the same cell types across biological conditions such as case vs. control, which reveals the expression change of a particular cell type under different contexts. (2) between cell types: comparing different cell types under the same condition, which identifies biomarkers to distinguish cell types. In either case, the experiment starts from a number of cells randomly picked from a tissue sample consisting of a mixture of different cell types. The only factor one can control is the total number of cells.

For two-group comparison

In the first scenario, the numbers of cells for a particular cell type in different biological conditions are often similar, barring significant changes in cell composition. It uses one cell type as the benchmark data and perturbs the genes with mixture proportion (π) and log fold changes (lfc) as the DE genes according two forms.

# Users can customize how many cells they need as to change the number of n. 
simData = Simulate2SCE(n=100, estParas1 = est_Paras, estParas2 = est_Paras)
de = runDE(simData$sce, DE_Method = "MAST")
estPower1 = Power_Disc(de, simData = simData)
estPower2 = Power_Cont(de, simData = simData)

For multiple-group comparisons

In the second scenario, the numbers of cells for distinct cell types can be very different, so the power for DE highly depends on the mixing proportions. Here, we use 1000*57 expression matrix sampled from the patient ID AB_S7 from the GSE67835 data. The original data GSE67835 can be download here: https://www.dropbox.com/s/55zdktqfqiwfs3l/GSE67835.RData?dl=0. In this case, we utilize the most abundant three cell types including oligodendrocytes, hybrid, and neurons. The underlying cell type proportion is 20%, 30%, 50% respectively. We simulate a series of datasets corresponding to different sample sizes. For each dataset, we then perform pairwise comparison and report the power evaluation for each comparison

data("GSE67835_AB_S7_sce")
sim_size = 200 # use a large sample is preferrable
cell_per = c(0.2, 0.3, 0.5)
col = colData(sce)
exprs = assays(sce)$counts
(tb = table(colData(sce)$Patients, colData(sce)$cellTypes))
# use AB_S7 patient as example and take three cell types: astrocytes hybrid and neurons
estParas_set = NULL
celltypes = c("oligodendrocytes", "hybrid", "neurons")
for (cp in celltypes){
    print(cp)
    ix = intersect(grep(cp, col$cellTypes), grep("AB_S7", col$Patients))
    tmp_mat = exprs[, ix]
    tmp_paras = Est2Phase(tmp_mat)
    estParas_set[[cp]] = tmp_paras
}
######### 
#########  Simulation part
######### 
sim = SimulateMultiSCEs(n = sim_size, estParas_set = estParas_set, multiProb = cell_per)

######### 
#########  DE analysis part
######### 
DE_rslt = NULL
for (comp in names(sim)){
    tmp = runDE(sim[[comp]]$sce, DE_Method = "MAST")
    DE_rslt[[comp]] = tmp
}

######### 
######### Summarize the power result
######### 
pow_rslt = pow1 = pow2 = pow1_marg = pow2_marg = NULL
TD = CD = NULL
for (comp in names(sim)){
    tmp1 = Power_Disc(DE_rslt[[comp]], sim[[comp]])
    tmp2 = Power_Cont(DE_rslt[[comp]], sim[[comp]])
    TD = c(TD, tmp2$TD); CD = c(CD, tmp2$CD)
    pow1_marg = c(pow1_marg, tmp1$power.marginal)
    pow2_marg = c(pow2_marg, tmp2$power.marginal)
    pow_rslt[[comp]] = list(pow1 = tmp1, pow2 = tmp2)
    pow1 = rbind(pow1, tmp1$power)
    pow2 = rbind(pow2, tmp2$power)
}

######### 
######### Demonstrate the result by heatmap
######### 
library(RColorBrewer); library(pheatmap)
breaksList = seq(0, 1, by = 0.01)
colors = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(length(breaksList))
dimnames(pow1) = list(names(sim), names(tmp1$CD))
dimnames(pow2) = list(names(sim), names(tmp2$CD))
pheatmap(pow2, display_numbers = TRUE, color=colors, show_rownames = TRUE,
         cellwidth = 30, cellheight = 40, legend = TRUE,
         border_color = "grey96", na_col = "grey",
         cluster_row = FALSE, cluster_cols = FALSE,
         breaks = seq(0, 1, 0.01),
         main = "")

Session Info

sessionInfo()

## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        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       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] pheatmap_1.0.12             RColorBrewer_1.1-3         
##  [3] POWSC_1.15.0                MAST_1.33.0                
##  [5] SingleCellExperiment_1.29.0 SummarizedExperiment_1.37.0
##  [7] GenomicRanges_1.59.0        GenomeInfoDb_1.43.0        
##  [9] IRanges_2.41.0              S4Vectors_0.45.0           
## [11] MatrixGenerics_1.19.0       matrixStats_1.4.1          
## [13] Biobase_2.67.0              BiocGenerics_0.53.1        
## [15] generics_0.1.3              BiocStyle_2.35.0           
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6            xfun_0.49               bslib_0.8.0            
##  [4] ggplot2_3.5.1           lattice_0.22-6          vctrs_0.6.5            
##  [7] tools_4.4.2             tibble_3.2.1            fansi_1.0.6            
## [10] highr_0.11              pkgconfig_2.0.3         Matrix_1.7-1           
## [13] data.table_1.16.2       lifecycle_1.0.4         GenomeInfoDbData_1.2.13
## [16] compiler_4.4.2          stringr_1.5.1           progress_1.2.3         
## [19] statmod_1.5.0           munsell_0.5.1           htmltools_0.5.8.1      
## [22] sys_3.4.3               buildtools_1.0.0        sass_0.4.9             
## [25] yaml_2.3.10             pillar_1.9.0            crayon_1.5.3           
## [28] jquerylib_0.1.4         limma_3.63.1            DelayedArray_0.33.1    
## [31] cachem_1.1.0            abind_1.4-8             digest_0.6.37          
## [34] stringi_1.8.4           reshape2_1.4.4          maketools_1.3.1        
## [37] fastmap_1.2.0           grid_4.4.2              colorspace_2.1-1       
## [40] cli_3.6.3               SparseArray_1.7.1       magrittr_2.0.3         
## [43] S4Arrays_1.7.1          utf8_1.2.4              prettyunits_1.2.0      
## [46] UCSC.utils_1.3.0        scales_1.3.0            rmarkdown_2.29         
## [49] XVector_0.47.0          httr_1.4.7              hms_1.1.3              
## [52] evaluate_1.0.1          knitr_1.48              rlang_1.1.4            
## [55] Rcpp_1.0.13-1           glue_1.8.0              BiocManager_1.30.25    
## [58] jsonlite_1.8.9          R6_2.5.1                plyr_1.8.9             
## [61] zlibbioc_1.52.0