1 Installation and help

1.1 Install FEAST

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

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

1.2 Help for FEAST

If you have any FEAST-related questions, please post to the GitHub Issue section of FEAST at https://github.com/suke18/FEAST/issues, which will be helpful for the construction of FEAST.

2 Introduction

2.1 Background

Cell clustering is one of the most important and commonly performed tasks in single-cell RNA sequencing (scRNA-seq) data analysis. An important step in cell clustering is to select a subset of genes (referred to as “features”), whose expression patterns will then be used for downstream clustering. A good set of features should include the ones that distinguish different cell types, and the quality of such set could have significant impact on the clustering accuracy. All existing scRNA-seq clustering tools include a feature selection step relying on some simple unsupervised feature selection methods, mostly based on the statistical moments of gene-wise expression distributions. In this work, we develop a feature selection algorithm named FEAture SelecTion (FEAST), which provides more representative features.

2.2 Citation

For consensus clustering, cite Strehl and Ghosh (2002).
For using (FEAST) feature selection procedure, Su, Yu, and Wu (2021)

2.3 Quick start

We use Yan dataset(Yan et al. 2013) for the demonstration. This Yan dataset includes 6 cell types about human preimplantation embryos and embryonic stem cells. The users can run the FEAST or FEAST_fast function to obtain the top ranked genes. These featured genes are deemed to better contribute to the clustering accuracy. Here, we demonstrate that some of the top 9 selected genes are very informative (or known as marker genes e.g., CYYR1 ).

library(FEAST)
data(Yan)
k = length(unique(trueclass))
Y = process_Y(Y, thre = 2)
ixs = FEAST(Y, k=k)
# look at the features
Ynorm = Norm_Y(Y)
par(mfrow = c(3,3))
for (i in 1:9){
  tmp_ix = ixs[i]
  tmp_gene = rownames(Ynorm)[tmp_ix]
  boxplot(as.numeric(Ynorm[tmp_ix, ])~trueclass, main = tmp_gene, xlab="", ylab="", las=2)
}

3 Using FEAST for scRNA-seq clustering analysis

FEAST requires a count expression matrix with rows representing the genes and columns representing the cells. To preprocess the count matrix, FEAST filters out the genes depending on the dropout rates (sometimes known as zero-expression rates). Here, we use one template Yan (Yan et al. 2013) dataset. It is important to note that FEAST requires the number of clusters (k) as an input, which can be obtained by the prior knowledge of the cell types.

3.1 Load the data

The Yan dataset is loaded with two objects including the count expression matrix (Y) and corresponding cell type information (trueclass). Other sample datasets can be downloaded at https://drive.google.com/drive/u/0/folders/1SRT7mrX7ziJoSjuFLLkK8kjnUsJrabVM. If the users want to use one Deng dataset(Deng et al. 2014), the users can load the data and access the phenotype information by using colData(Deng) function, and the count expression matrix by using assay(Deng,"counts") function from the SummarizedExperiment Bioconductor package.

data(Yan)
dim(Y)
#> [1] 19304   124
table(trueclass)
#> trueclass
#>  2-cell  4-cell  8-cell    Late Morulae  Oocyte  Zygote    hESC 
#>       6      12      20      30      16       3       3      34

3.2 Consensus clustering

To preprocess the count expression matrix (Y), FEAST filters out the genes based on the dropout rate (sometimes known as zero-expression rate). This consensus clustering step will output the initial clustering labels based on a modified CSPA algorithm. It will find the cells that tightly clustered together and possibly filter out the cells that are less correlated to the initial clustering centers.

Y = process_Y(Y, thre = 2) # preprocess the data if needed
con_res = Consensus(Y, k=k)

3.3 Calculate the gene-level significance

After the consensus clustering step, FEAST will generate the initial clustering labels with confidence. Based on the initial clustering outcomes, FEAST further infers the gene-level significance by using F-statistics followed by the ranking of the full list of genes.

F_res = cal_F2(Y, con_res$cluster)
ixs = order(F_res$F_scores, decreasing = TRUE) # order the features

3.4 Clustering and validation

After ranking the genes by F-statistics, FEAST curates a series of top number of features (genes). Then, FEAST input these top (by default top = 500, 1000, 2000) features into some established clustering algorithm (e.g. SC3 (Kiselev et al. 2017) and TSCAN(Ji and Ji 2016)). It is worth noting that one needs to confirm the k to be the same for every iteration.

After the clustering steps, with the case of unknown cell types, FEAST evaluates the quality of the feature set by computing the average distance between each cell and the clustering centers, which is the same as the MSE. It is noting that FEAST uses all the features (genes) for distance calculation for the purpose of fair comparisons. The clustering centers are obtained by using previous clustering step.

It is noted that FEAST is not limited to SC3 clustering. It also shows superior performance on other clustering methods such as TSCAN. The code for running SC3 and TSCAN are embedded in FEAST package, which can be access by SC3_Clust and TSCAN_Clust.

## clustering step
tops = c(500, 1000, 2000)
cluster_res = NULL
for (top in tops){
    tmp_ixs = ixs[1:top]
    tmp_markers = rownames(Y)[tmp_ixs]
    tmp_res = TSCAN_Clust(Y, k = k, input_markers = tmp_markers)
    cluster_res[[toString(top)]] = tmp_res
}
## validation step
Ynorm = Norm_Y(Y)
mse_res = NULL
for (top in names(cluster_res)){
    tmp_res = cluster_res[[top]]
    tmp_cluster = tmp_res$cluster
    tmp_mse = cal_MSE(Ynorm = Ynorm, cluster = tmp_cluster)
    mse_res = c(mse_res, tmp_mse)
}
names(mse_res) = names(cluster_res)

3.5 Compare to the real cell type labels

3.5.1 Benchmark with the original SC3

After the validation step, the feature set associated with the smallest MSE value will be recommended for the further analysis. Here, we demonstrate that an improved clustering accuracy by specifying the optimal feature set. The benchmark comparison is the original setting of established clustering algorithm (e.g. SC3).

original = TSCAN_Clust(Y, k=k, input_markers = rownames(Y))
id = which.min(mse_res)
eval_Cluster(original$cluster, trueclass)
eval_Cluster(cluster_res[[id]]$cluster, trueclass)

3.5.2 Show the clustering improvement by using figures

As demonstrated in the figure below, the first panel includes the PCA illustration of the cell types. The second panel shows the clustering result by the original SC3, in which some cells (inside the gray circle) are mixed. The third panel shows the improved clustering outcomes by inputing the optimized feature set into SC3. This is an example for the Deng dataset.

4 Quick use by the wrapper function

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.

data(Yan)
Y = process_Y(Y, thre = 2) # preprocess the data if needed
con_res = Consensus(Y, k=k)
mod_res = Select_Model_short_TSCAN(Y, cluster = con_res$cluster, top = c(200, 500, 1000))
# to visualize the result, one needs to load ggpubr library. 
library(ggpubr)
PP = Visual_Rslt(model_cv_res = mod_res, trueclass = trueclass)
print(PP$ggobj) # show the visualization plot.

The result figure shows that the MSE validation process is concordant with the clustering accuracy trend. It demonstrated that using top 1000 featured genes can result in the highest clustering accuracy with improvement than the original setting.

5 Session Info

sessionInfo()
#> R version 4.1.0 RC (2021-05-10 r80283)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows Server x64 (build 17763)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=C                          
#> [2] LC_CTYPE=English_United States.1252   
#> [3] LC_MONETARY=English_United States.1252
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.1252    
#> 
#> attached base packages:
#> [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
#> [8] methods   base     
#> 
#> other attached packages:
#>  [1] ggpubr_0.4.0                ggplot2_3.3.3              
#>  [3] FEAST_1.0.0                 SummarizedExperiment_1.22.0
#>  [5] Biobase_2.52.0              GenomicRanges_1.44.0       
#>  [7] GenomeInfoDb_1.28.0         IRanges_2.26.0             
#>  [9] S4Vectors_0.30.0            BiocGenerics_0.38.0        
#> [11] MatrixGenerics_1.4.0        matrixStats_0.58.0         
#> [13] BiocParallel_1.26.0         mclust_5.4.7               
#> [15] BiocStyle_2.20.0           
#> 
#> loaded via a namespace (and not attached):
#>   [1] readxl_1.3.1                snow_0.4-3                 
#>   [3] backports_1.2.1             plyr_1.8.6                 
#>   [5] igraph_1.2.6                splines_4.1.0              
#>   [7] fastICA_1.2-2               digest_0.6.27              
#>   [9] foreach_1.5.1               htmltools_0.5.1.1          
#>  [11] magick_2.7.2                fansi_0.4.2                
#>  [13] magrittr_2.0.1              cluster_2.1.2              
#>  [15] doParallel_1.0.16           ROCR_1.0-11                
#>  [17] openxlsx_4.2.3              colorspace_2.0-1           
#>  [19] rrcov_1.5-5                 WriteXLS_6.3.0             
#>  [21] haven_2.4.1                 xfun_0.23                  
#>  [23] dplyr_1.0.6                 crayon_1.4.1               
#>  [25] RCurl_1.98-1.3              jsonlite_1.7.2             
#>  [27] iterators_1.0.13            glue_1.4.2                 
#>  [29] gtable_0.3.0                zlibbioc_1.38.0            
#>  [31] XVector_0.32.0              DelayedArray_0.18.0        
#>  [33] car_3.0-10                  SingleCellExperiment_1.14.0
#>  [35] DEoptimR_1.0-8              abind_1.4-5                
#>  [37] scales_1.1.1                pheatmap_1.0.12            
#>  [39] mvtnorm_1.1-1               DBI_1.1.1                  
#>  [41] rngtools_1.5                rstatix_0.7.0              
#>  [43] Rcpp_1.0.6                  TrajectoryUtils_1.0.0      
#>  [45] xtable_1.8-4                foreign_0.8-81             
#>  [47] proxy_0.4-25                gplots_3.1.1               
#>  [49] RColorBrewer_1.1-2          ellipsis_0.3.2             
#>  [51] pkgconfig_2.0.3             farver_2.1.0               
#>  [53] sass_0.4.0                  utf8_1.2.1                 
#>  [55] labeling_0.4.2              tidyselect_1.1.1           
#>  [57] rlang_0.4.11                later_1.2.0                
#>  [59] munsell_0.5.0               cellranger_1.1.0           
#>  [61] tools_4.1.0                 generics_0.1.0             
#>  [63] broom_0.7.6                 evaluate_0.14              
#>  [65] stringr_1.4.0               fastmap_1.1.0              
#>  [67] yaml_2.2.1                  knitr_1.33                 
#>  [69] zip_2.1.1                   robustbase_0.93-7          
#>  [71] caTools_1.18.2              purrr_0.3.4                
#>  [73] nlme_3.1-152                doRNG_1.8.2                
#>  [75] mime_0.10                   compiler_4.1.0             
#>  [77] curl_4.3.1                  e1071_1.7-6                
#>  [79] ggsignif_0.6.1              tibble_3.1.2               
#>  [81] bslib_0.2.5.1               pcaPP_1.9-74               
#>  [83] stringi_1.6.2               highr_0.9                  
#>  [85] forcats_0.5.1               lattice_0.20-44            
#>  [87] TSCAN_1.30.0                Matrix_1.3-3               
#>  [89] ggsci_2.9                   vctrs_0.3.8                
#>  [91] pillar_1.6.1                lifecycle_1.0.0            
#>  [93] BiocManager_1.30.15         combinat_0.0-8             
#>  [95] jquerylib_0.1.4             data.table_1.14.0          
#>  [97] cowplot_1.1.1               bitops_1.0-7               
#>  [99] irlba_2.3.3                 httpuv_1.6.1               
#> [101] R6_2.5.0                    bookdown_0.22              
#> [103] promises_1.2.0.1            KernSmooth_2.23-20         
#> [105] rio_0.5.26                  codetools_0.2-18           
#> [107] gtools_3.8.2                assertthat_0.2.1           
#> [109] withr_2.4.2                 GenomeInfoDbData_1.2.6     
#> [111] mgcv_1.8-35                 hms_1.1.0                  
#> [113] grid_4.1.0                  tidyr_1.1.3                
#> [115] class_7.3-19                rmarkdown_2.8              
#> [117] SC3_1.20.0                  carData_3.0-4              
#> [119] shiny_1.6.0

Reference

Deng, Qiaolin, Daniel Ramsköld, Björn Reinius, and Rickard Sandberg. 2014. “Single-Cell Rna-Seq Reveals Dynamic, Random Monoallelic Gene Expression in Mammalian Cells.” Science 343 (6167): 193–96.

Ji, Zhicheng, and Hongkai Ji. 2016. “TSCAN: Pseudo-Time Reconstruction and Evaluation in Single-Cell Rna-Seq Analysis.” Nucleic Acids Research 44 (13): e117–e117.

Kiselev, Vladimir Yu, Kristina Kirschner, Michael T Schaub, Tallulah Andrews, Andrew Yiu, Tamir Chandra, Kedar N Natarajan, et al. 2017. “SC3: Consensus Clustering of Single-Cell Rna-Seq Data.” Nature Methods 14 (5): 483–86.

Strehl, Alexander, and Joydeep Ghosh. 2002. “Cluster Ensembles—a Knowledge Reuse Framework for Combining Multiple Partitions.” Journal of Machine Learning Research 3 (Dec): 583–617.

Su, Kenong, Tianwei Yu, and Hao Wu. 2021. “Accurate Feature Selection Improves Single-Cell Rna-Seq Cell Clustering.” Briefings in Bioinformatics.

Yan, Liying, Mingyu Yang, Hongshan Guo, Lu Yang, Jun Wu, Rong Li, Ping Liu, et al. 2013. “Single-Cell Rna-Seq Profiling of Human Preimplantation Embryos and Embryonic Stem Cells.” Nature Structural & Molecular Biology 20 (9): 1131.

The FEAST User’s Guide

19 May 2021

Abstract

Package