FEM: mining biological meaning from cell level in single-cell RNA sequencing data

Motivation for this approach

Functional Gene Set Enrichment Analysis (GSEA) is often the last step of expression data analysis. A gene set usually represents a biological function, and the function will be used to refer to a gene set in subsequent analyses. Many R packages and some online sites, such as DAVID, provide enrichment analysis tools. The software takes a set of genes given by the user as input and returns functions (pathways) that are significantly enriched. Enrichment analysis methods for individual samples have also been proposed (Foroutan et al., 2018; Hänzelmann, Castelo & Guinney, 2013).

Zhang et al. (2020) summarized the widely used GSEA algorithms. We divided their methods into three categories: methods based on dimension reduction, methods based on single samples, and methods based on phenotype or cluster.

Methods based on dimension reduction include Pagoda2 (Fan et al., 2016) and PLAGE (Tomfohr, Lu & Kepler, 2005). These methods convert the gene expression matrix into FEM by performing dimension reduction (SVD for PLAGE and PCA for Pagoda2) on each gene set (pathway). Dimension reduction uses all samples and is not entirely based on a single cell.

Methods based on phenotype or cluster include Vision (DeTomaso et al., 2019) and z-score (Lee et al., 2008). Vision methods clusters GEM into clusters by KNN, and GSEA is based on these clusters. The Z-score method depends on phenotypes, which define each pathway activity score based on the top gene of t-test score between two conditions. These gradually join the genes in pathways to maximize the activity score and are not based on a single sample.

Methods based on the single sample include AUCell (Aibar et al., 2017), ssGSEA (Barbie et al., 2009), and GSVA (Hänzelmann, Castelo & Guinney, 2013), which are all non-parametric test statistical methods based on gene rank. Data sparsity is the main problem when analysing data from a single cell. Therefore, the level of gene expression is not the primary consideration for a pathway, but rather, whether the gene is expressed. A method based on Fisher’s exact test or the chi-square test may have better results than the no-parameter test method based on ranking.

scRNA-seq has low capture efficiency and high dropouts due to the small amount of initial material used. scRNA-seq expression data tend to be sparse (Hicks et al., 2018; Eraslan et al., 2019). A true zero represents the lack of expression of a gene in a specific cell due to the non-trivial distinction between true- and false-zero counts. Meanwhile, a false zero is a technical deviation. The RNA-seq technique only detects mRNA molecules that are present, therefore, a gene in the scRNA-seq dataset with a non-zero expression value indicates that at least one mRNA molecule is present. These sparse non-zero expression values provided inspiration for the possibility of performing functional GSEA analysis at the single-cell level.

Functional enrichment analysis may be performed at the single-cell transcriptome expression level given the characteristics of scRNA-seq data. The algorithm is as follows: first, replace the initial expression matrix with a FEM that allows for the direct exploration of the differences between cells from different biological perspectives. Second, the top variable function of a cell cluster obtained by a gene expression matrix (GEM) can be represented as a cell scatter plot. Third, the coverage of a functional gene set is more important than the expression of individual genes for GSEA since genes often overlap different functional groups, which makes single high-expression genes difficult to explain. For example, the activation of a signalling pathway is the result of interactions between all genes; the high expression of one gene does not mean that the pathway is activated.

The current factors affecting the discovery of cell groups are as follows: first, technical covariates must be regressed out before downstream analysis can occur since these factors introduce systematic error and confound the technical and biological variability. This error may lead to systematic differences in gene expression profiles between batches (Leek et al., 2010; Hicks et al., 2018; Chen, Ning & Shi, 2019). The most prominent technical covariates in single-cell data are count depth and batch. The FEM method is based on a gene set, which makes it more robust than the GEM. The measured expression levels for the genes related to a certain biological process may fluctuate between batches or single cells due to variations in sequencing processing. That is, detection for one specific high-expression gene in single cell has a certain randomness. If we can consider the expression of a set of genes related to a specific process, we can increase the detection stability for the biological process of single cells, since the number of expressed genes of a specific biological process obey a hypergeometric distribution. Second, some biological effects may affect the results of the cluster algorithms. For example, the cell cycle phase of a given cell within a cell type population may cause separation from the same cells in the clustering analysis. However, correcting for biological covariates is not always helpful in interpreting scRNA-seq data (Kolodziejczyk et al., 2015). These influencing factors often do not have uniform filtering criteria. In some cases, the cell cycle may be part of the study, or there may be a relationship between the cell cycle and other functions (Haghverdi, Buettner & Theis, 2015; Vento-Tormo et al., 2018; McDavid, Finak & Gottardo, 2016; Blasi et al., 2017). The FEM method can be used to systematically survey the biological aspect of each cell before downstream analysis. To this end, we developed an scRNA-seq FEM algorithm.

Workflow of the proposed methods

The FEM was divided into four steps (Fig. 1). Step 1: The standardized scRNA-seq GEM was transformed into a FEM by multi-module gene enrichment analysis (gene ontology, pathway). Step 2: The p-value represents the significance of enrichment, therefore, the p-value obtained by enrichment analysis was converted into information content (FEM). Step 3: FEM/GEM was used for data standardization, dimensionality reduction clustering, and UMAP visualization (McInnes et al., 2018). Step 4: Differentially expressed genes (DEG) and differentially expressed functions (DEF) analysis was performed based on GEM clustering (see the “GEM and FEM fusion analysis” section for details).

Dataset

Peripheral blood mononuclear cells (PBMCs) are populations of immune cells that remain in the less dense upper interface of the Ficoll layer. PBMCs include lymphocytes (T cells, B cells, and NK cells), monocytes, and dendritic cells (DCs). The frequencies of these populations vary across individuals in humans. Lymphocytes are typically in the range of 70–90%, monocytes range from 10–30%, while DCs are only present at 1–2% (Norman, 1995). The PBMCs dataset used in this study was downloaded from the 10X Genomics official website (https://www.10xgenomics.com/resources/datasets), which included B cells, NK cells, CD8 T cells, memory CD4 T cells, naïve CD4 T cells, DC, CD14+ monocytes, FCGR3A+ monocytes, and a small number of platelets. This dataset contained 2,700 cells in total.

The human pancreas dataset contained 2,126 cells and 10 cell types. This included alpha, ductal, endothelial, delta, acinar, beta, gamma, mesenchymal, and epsilon cells, as well as a small number of unknown cell types (Muraro et al., 2016).

The human liver dataset consisted of 777 cells, including seven types of cells: definitive endoderm cells, immediate hepatoblast cells, induced pluripotent stem cells (IPSCs), material hepatocytic cells, hepatic endoderm cells, endothelial cells, and mesenchymal stem cells (Camp et al., 2017).

FEM algorithm

Algorithm optimization

Fisher’s exact test is a time-consuming process. For single-cell data, a statistical test would be required for each function of each cell. Therefore, when the number of cells is large, the computation time would be untenable. Therefore, the algorithm was optimized with the addition of multi-core computing.

Optimization of the algorithm can be illustrated using the symbols in section “Gene-functional group conversion”. First, because N and K were invariable for all cells, these values were stored in memory to avoid recalculation with each replicate. Second, the gene expression matrix was transformed into 0,1 matrix A, where 1 represents the expression of gene i in cell j and 0 represents the non-expression of gene i in cell j. The gene set was also transformed into 0,1 matrix B, where 1 represents the presence of gene i in set s and 0 represents the absence of gene i in set s. The element (k) in the product A × B^T of the two matrices represents the number of genes expressed by cell j in set s (Fig. 2).

Cluster and differentially expressed gene detection based on FEM and integration of GEM and FEM by Seurat

Evaluation of FEM clustering results

We defined an evaluation score, SC,  to compare the overlap between clusters based on FEM clustering and real cell types. The number of clusters was related to specific parameters, therefore, each cluster should contain only one cell type as much as possible. Additionally, the cluster’s number should be greater than or equal to the number of cell types. Therefore, the clustering parameters were set so that the number of clusters was greater than or equal to the number of actual cell types. For each cluster i, SC_i was measured by the consistency of the cell type within the cluster (the proportion of the cell type with the largest number of cells to the total number of cells in the cluster). We used the following formula to evaluate clustering results: (3)${\displaystyle S{C}_i=\frac{c_{i,max}}{c_{i,T}}}$

c_i,max represents the cell type with the largest number of cells in the i-th cluster. c_i,T represents the total number of all cells in the i-th cluster.

Data availability

Publicly available scRNA-seq datasets were used in this study. The PBMCs dataset was downloaded from the 10X Genomics dataset page (https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz). The liver and pancreas scRNA-seq data can be accessed in the NCBI Gene Expression Omnibus (GEO) under accession numbers: GSE81252, GSE85241.

Code availability

The FEM software and related scripts have been deposited in the GitHub repository (https://github.com/qingyunpkdd/single_cell_fem).

FEM can separate different cell types

The Reactome pathway gene sets have only 5,741 unique genes and a large amount of gene information was lost when using the FEM algorithm. In comparison GO gene sets had 15,578 unique genes. Therefore, GO gene set-based FEM was used for cluster analysis. The Reactome pathway gene sets were used in the analysis of GC-DEF and their cluster results are shown in Fig. S1. GO-based FEM results are shown in Figs. 3–5.

Three data sets (PBMCs, liver, pancreas) were used to replace GEM for dimensionality reduction and clustering to verify whether the GO-based FEM algorithm could separate different cell types. The number of clusters was artificially adjusted and the number of clusters was set to be greater than or equal to the number of actual cell types. This was done so that a cluster contained only one main cell type, which allowed us to better distinguish the different cell types.

We needed to discover different cell groups and determine the relationship between different cell types for single-cell scRNA-seq data. These cells may have similar functions or a common recent differentiation origin. The FEM method may help determine functionally similar cells. Cells were noted to be different subtypes of the same cell if the two groups of cells were far apart in GEM, but close or partially overlapped in FEM. This result may also indicate that the two groups of cells may perform similar functions.

We performed GEM clustering and GO-based FEM clustering on the scRNA-seq data, and compared the two clusters. The FEM clustering results were also compared with real cell types to determine the similarity of cell functions. A low SC score in a particular cluster i indicated that the cluster may contain multiple cell types. A moderate number of a certain cell type A in the cluster i was defined as the main cell type MC_i. For a cluster, all MCs may have the same or similar functions.

The results show that the FEM method can distinguish different types of cells that have different functions on three data sets (Figs. 3–5). Figures 3–5 show that FEM can find cells with similar origins or different subgroups of a cell. Figure 3 shows three clusters (G1, G2, G3) that have a low SC score and all of the T cells overlapped in cluster G3. The MC_i of NK cells and T cells in cluster G2 literature show that they have a common differentiation origin T/NK cell progenitor. G1’s MC_i contained only monocytes (Lim et al., 2013). The FEM in Fig. 4 helped determine the relationship of definitive endoderm cells and hepatic endoderm cells (Camp et al., 2017). Figure 5 illustrates that alpha, beta, gamma, epsilon, and delta cells form one cluster in the FEM clusters, and that the cell differentiation tree verifies that they have a common endocrine progenitor (Jacobson & Tzanakakis, 2017).

FC-DEF can directly detect the functional differences between clusters

GC-DEF can reveal the expression of specific functions from the UMAP(t-SNE) diagram of GEM and the results of GC-DEF may aid GEM in determining the functions of a specific group. Official data pre-processing in the PBMC dataset included the removal of cells with excessive mitochondrial genes (>5%, quality control) and those with too many (>2,500) or too few (<200) features (Bittersohl & Steimer, 2016). Some cells may have been filtered out in the data pre-processing stage whereas our method did not filter out any cells. Indeed, the filtered cells were found to be located above the NC cells. Using the GC-DEF method, this cell group was found to have a highly expressed cell proliferation-related pathway. Thus, these cells were identified as proliferative (Fig. 6).

The most significantly expressed pathway, “hemostasis”, was located in the platelet clusters. The “innate immune system” pathway was significantly expressed in the monocyte, DC cell, NK cell, and platelet cell populations, which was consistent with literature results (NORMAN, 1995). The top five highly expressed functions were consistent with the cell type characteristics. The other FC-DEF and GC-DEF results are provided in Table S1.

Table 2 shows that most of the corresponding cells were closely related to their corresponding top functions, such as the high expression of “Reactome regulation of beta cell development” in beta cells and the high expression of “Reactome gluconeogenesis” in mature stem cells. All other results are shown in Table S1.

GO-based GC-DEF results was consistent with the pathway-based methods and literature results (Norman, 1995). All other results are presented in Table 3 and Tables S1–S4.

Difference between GO enrichment analysis results based on DEG and FEM

The general GSE analysis was based on DEGs between groups (cell types or clusters). The FEM-based gene set (function) enrichment can be defined as the DEF between groups. The results (Table S5, p-value threshold ≤ 0.05) show that: (1) the enrichment analysis results obtained by the two methods are quite different. (2) There are more GO items based on the FEM method than the method based on DEGs.

The possible explanation for the first point is that the two methods are based on different data, and the FEM method is an enrichment result of a single cell. The DEG method is based on results from cell clusters. FEM uses all expressed genes in a cell to calculate GO enrichment, while the DEG method only uses those differentially expressed genes, which may explain why there are more GO items based on the FEM method than the DEGs method. The FEM method focuses on the “differentially expressed GO items”, and the DEG method focuses on “the enriched GO items of differentially expressed genes”.

Validation with an immune dataset

The Immunologic Signatures Collection (ImmuneSigDB) (Liberzon et al., 2011) was used as a validation dataset to test whether the proposed method could detect sets of genes that had been identified as up-regulated or down-regulated by traditional methods. The ImmuneSigDB is composed of gene sets that represent cell types, states, and perturbations within the immune system (Godec et al., 2016). Figure 7 shows that within the bulk RNA-seq dataset, the cell types from the up-regulated expression marker gene set were highly expressed using the method proposed in this study. This demonstrated the efficacy of the proposed method for detecting cell-type-specific gene sets.