This vignette introduces the usage of the Bioconductor package FEAST (FEAture SelecTion for scRNA-seq data), which is specifically designed for selecting most representative genes before performing the core of clustering. It is demonstrated an improved clustering accuracy when recruting the featured genes selected by FEAST.
FEAST 1.8.0
To install this package, start R (version > “4.0”) and enter:
if(!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("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.
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.
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 gene rankings (from the most significant to the least significant genes) of 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)
# The core function.
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)
}
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). Alternatively, the users can input the normalized matrix (usually log transformed of relative abundance expression matrix (scaled by average sequencing depth)). 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.
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
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)
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
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 such as SC3 (Kiselev et al. 2017), which is regarded as the most accurate scRNA-seq clustering algorithm (Duò, Robinson, and Soneson 2018). 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)
#tmp_res = SC3_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)
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)
id = which.min(mse_res)
eval_Cluster(original$cluster, trueclass)
eval_Cluster(cluster_res[[id]]$cluster, trueclass)
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.