scMuffin: an R package to disentangle solid tumor heterogeneity by single-cell gene expression analysis

Abstract

Introduction

Single-cell (SC) gene expression analysis is crucial to dissect the complex cellular heterogeneity of solid tumors, which is one of the main obstacles for the development of effective cancer treatments. Such tumors typically contain a mixture of cells with aberrant genomic and transcriptomic profiles affecting specific sub-populations that might have a pivotal role in cancer progression, whose identification eludes bulk RNA-sequencing approaches. We present scMuffin, an R package that enables the characterization of cell identity in solid tumors on the basis of a various and complementary analyses on SC gene expression data.

Results

scMuffin provides a series of functions to calculate qualitative and quantitative scores, such as: expression of marker sets for normal and tumor conditions, pathway activity, cell state trajectories, Copy Number Variations, transcriptional complexity and proliferation state. Thus, scMuffin facilitates the combination of various evidences that can be used to distinguish normal and tumoral cells, define cell identities, cluster cells in different ways, link genomic aberrations to phenotypes and identify subtle differences between cell subtypes or cell states. We analysed public SC expression datasets of human high-grade gliomas as a proof-of-concept to show the value of scMuffin and illustrate its user interface. Nevertheless, these analyses lead to interesting findings, which suggest that some chromosomal amplifications might underlie the invasive tumor phenotype and the presence of cells that possess tumor initiating cells characteristics.

Conclusions

The analyses offered by scMuffin and the results achieved in the case study show that our tool helps addressing the main challenges in the bioinformatics analysis of SC expression data from solid tumors.

Background

Single-cell (SC) gene expression analysis is crucial to dissect the complex cellular heterogeneity of solid tumors, which is one of the main obstacles for the development of effective cancer treatments [1]. A relevant number of software tools has been developed in recent years in the field of SC data analysis [2], a fact that stresses the key opportunities and challenges in this field. A recent study has shown that the development of tools that address common tasks (e.g., clustering of similar cells) and ordering of cells (e.g., definition of cell trajectories) is decreasing, while a greater focus is being paid on data integration and classification [2]. These observations reflect the growing availability, scale and complexity of SC datasets [2].

SC datasets of solid tumors are typical examples of complex datasets that present a series of computational challenges and whose analysis demands domain-specific and integrative approaches. In fact, solid tumors typically contain a mixture of cells with aberrant genomic and transcriptomic profiles affecting specific sub-populations that might play a pivotal role in cancer progression, whose identification eludes bulk RNA-sequencing approaches. The use of cell type-specific markers (when available) is limited, and the alterations of gene expression that mark cancer cells makes the use of markers for normal cells not completely adequate. Moreover, the molecular heterogeneity of cancer cells (due to both intra-tumor and inter-individual differences) poses intrinsic limits in the definition of such markers. In addition, solid tumor samples typically comprise cells from the surrounding tissue or infiltrating cells that need to be distinguished from tumor populations for an effective analysis. Another challenge is the identification of clinically relevant cell subtypes that may be rare in the tumor mass, such as cancer stem cells or drug resistant subclones: because of their relatively low number, these cells are typically clustered together with many others. Lastly, an intrinsic problem of many SC datasets is the sequencing depth limit at the SC level. These limitations bound the number of detectable genes to the few thousands of the highest expressed genes, which implies, for example, that some established markers may not be used for data analysis.

To address these challenges, we developed scMuffin, an R package that implements a series of complementary analyses aimed at shedding light on the complexity of solid tumor SC expression data, including: a fast and customizable gene set scoring system; gene sets from various sources, including pathways, cancer functional states and cell markers; cell cluster association analysis with quantitative (e.g., gene set scores) as well as categorical (e.g. mutation states, proliferation states) features; copy number variation (CNV) analysis; transcriptional complexity analysis; proliferation rate quantification; and gene set based multi-dataset analysis (Fig. 1). scMuffin facilitates the integrative analysis of these multiple features, thus allowing the identification of cell subtypes that elude more general clustering and classification approaches. We describe the key aspects of scMuffin implementation and then its user interface by means of a case study on a public SC expression dataset of human high-grade gliomas (HGG).

Fig. 1

Overview of scMuffin package. scMuffin offers the possibility to perform several different analyses and data integration approaches to address the main challenges of SC gene expression analysis in solid tumors

Implementation

scMuffin is implemented as an R package and provides a series of functions that allows the user to quantify various genetic and phenotypic characteristics of single cells, which can be combined to obtain insights on cell identity and function. The functions of scMuffin operate on a common data structure, the “scMuffinList”, by adding, changing or removing its elements. Computationally intensive tasks (in particular, gene set scoring and CNV inference) are parallelized. Package documentation includes R vignettes and function documentation, which are also available as “GitHub Pages” at the URL https://emosca-cnr.github.io/scMuffin. In this section, we describe package inputs and the main algorithms and definitions underlying the analyses offered by scMuffin.

Input

scMuffin is intended to be used downstream to general purpose tasks like quality control, normalization, cell clustering and dataset integration, for which there are dedicated tools, such as Seurat [3]. scMuffin requires the genes-by-cells raw counts matrix, genes-by-cells normalized expression matrix and a partition of cells in clusters. Typically, these matrices are filtered during preprocessing steps to exclude low quality cells and genes that could negatively affect the analyses. However, the characterization of cells that can be achieved with scMuffin offers insights that can be used to (further) filter the dataset and/or to decide on which cells apply particular analyses (e.g. biomarker identification). In general, according to research questions and experimental design one may want to apply strong or mild filters before using scMuffin.

Quantification of gene set expression scores at cell and cluster levels

The quantification of gene set expression scores follows the approach described in [4, 5], in which a gene set is scored on the basis of its average deviation from an empirical null. The approach is implemented enabling the user to tune various parameters (e.g., the minimum number of cells in which a gene set must be expressed, the number of bins of the null model and the possible removal of missing values) to better address the needs of a particular study in relation to the characteristics of the dataset under analysis. In particular, the code is implemented to handle missing values, which are one of the main issues in SC datasets. For example, this is important to quantify the expression of gene sets that are expressed just in a small subset of cells (and have null values elsewhere), or to account for the expression of different members of the same pathway in different cells.

The algorithm for gene set scoring is the following. Given a gene set $S$:

• 1. the genes occurring in the normalized genes-by-cells matrix are grouped into a series of bins according to their average expression across cells;

• 2. a number $k$ of random gene sets $S_{i}^{*}$: are created, of the same size of $S$, tossing genes from the same bins of $S$, in order to match the distribution of gene expression of each $S_{i}^{*}$ with that of $S$;

• 3. the averages $m_{c}$ and $m_{ic}^{*}$ are calculated, respectively, over the values of $S$ and $S_{i}^{*}$ in every cell $c$;

• 4. the expression score $Y_{c}$ is calculated as the average difference between $m_{c}$ and $m_{ic}^{*}$:

  $$Y_{c} = \frac{{\mathop \sum \nolimits_{i}^{k} \left( {m_{c} - m_{ic}^{*} } \right)}}{k}$$

• 5. the difference between the medians of $m_{c}$ and $m_{ic}^{*}$ values in any cluster is used as the cluster-level score of $S$.

CNV estimation based on adjacent gene expression

CNV inference in scMuffin is based on the “adjacent gene windows” approach [4, 6], which has been validated using both single nucleotide polymorphism arrays [6] and whole-exome sequencing [4]. The approach is implemented in parallel and, like gene set scoring, offers various parameter tuning and data filtering possibilities, which allow the investigator to optimize the analysis to specific needs. The CNV profile of each cell is calculated as a moving average of scaled gene expression levels ordered by genomic location. scMuffin offers the possibility of subtracting a “normal” reference profile to highlight relative CNVs. Since a SC dataset derived from a sample of a solid tumor can either contain or not a proper set of physiological cells to be used as reference, this reference can be derived from cells that are part of the dataset under analysis (as described in [4]), or using external dataset (as described in [6]). The main steps are:

• 1. the reference cells are added to the genes-by-cells matrix (optional);

• 2. the expression of each gene is scaled subtracting its average (optional);

• 3. the gene expression matrix is ordered by chromosome and gene location;

• 4. in each chromosome h, the estimated copy number $V_{ic}$ of cell c is calculated for all the ordered genes $i \in \left[ {\frac{w}{2} + 1,n_{h} - \frac{w}{2}} \right]$:

  $$V_{ic} = \mathop \sum \limits_{{j = i - \frac{w}{2}}}^{{i + \frac{w}{2}}} \frac{{e_{jc} }}{w + 1}$$

where w is an even number that defines the window size, that is, the number of genes located before and after gene i which contribute to the estimation of $V_{ic}$, $n_{h}$ is the number of genes in h, and $e_{jc}$ is the gene expression value;

• 5. $V_{ic}$ values are scaled subtracting their average in a cell (optional);

• 6. cells are clustered by their CNV profile;

• 7. the average CNV profile of the normal reference cells is subtracted from all the CNV profiles (optional).

• 8. the CNV score of a cell is calculated as:

  $$C_{c} = \mathop \sum \limits_{i} \left( {V_{ic} } \right)^{2}$$

Cells are clustered based on their CNV profile using Seurat [3]. A CNV region of a cluster is composed of the adjacent genes $i$ that have absolute median $V_{ic}$ (within the cluster) higher than a given threshold. This threshold is by default defined as the standard deviation of $V_{ic}$ values of all cells.

Transcriptional complexity

The transcriptional complexity (TC) is calculated using three different approaches. The TC-R $T_{c}^{R}$ of a cell is inferred based on the number of expressed genes $(g_{c} )$ over the number of total mapped reads (transcripts) $(t_{c} )$:

$$T_{c}^{R} = \frac{{g_{c} }}{{t_{c} }} \cdot \frac{{\mathop {\max }\limits_{c} \left( {t_{c} } \right)}}{{\mathop {\max }\limits_{c} \left( {g_{c} } \right)}}$$

where $g_{c} = \# \left\{ {r_{ic} \ge \alpha } \right\}$, $\alpha$ is a threshold over the gene count $r_{ic}$ and defines the gene i as “expressed”, and $t_{c} = \sum\nolimits_{i} {r_{ic} }$. Values of $T_{c}^{R}$ higher (lower) than 1 indicate a higher (lower) number of expressed genes in relation to the total transcripts of the cell.

The TC-LMR (linear model residual) of a cell corresponds to its residual in the linear regression model between $g_{c}$ and $t_{c}$:

$${\text{log}}_{10} \left( {g_{c} } \right) = \beta_{0} + \beta_{1} {\text{log}}_{10} \left( {t_{c} } \right) + \varepsilon_{c}$$

$$T_{c}^{LMR} = {\text{log}}_{10} \left( {g_{c} } \right) - \left[ {\beta_{0} + \beta_{1} {\text{log}}_{10} \left( {t_{c} } \right)} \right].$$

Positive (negative) residuals indicate higher (lower) number of expressed genes in relation to total transcripts.

The TC-H (entropy) is quantified as the transcriptional entropy of the cell [7]:

$$T_{c}^{H} = - \mathop \sum \limits_{ic} \left[ {p_{ic} {\text{ log}}\left( {p_{ic} } \right)} \right]$$

where $p_{ic} = r_{ic} /t_{c}$. The higher the entropy the higher the number of expressed genes.

Proliferation rate

The proliferation rate $P_{c}$ of a cell is quantified as the maximum between the two gene set scores $Y_{c} \left( {S_{1} } \right)$ and $Y_{c} \left( {S_{2} } \right)$, calculated on the gene sets $S_{1}$ and $S_{2}$ that contain, respectively, genes involved in G1/S and G2/M cell cycle phases:

$$P_{c} = {\text{max}}\left( {Y_{c} \left( {S_{1} } \right), Y_{c} \left( {S_{2} } \right)} \right)$$

where $S_{1}$ and $S_{2}$ are defined as in Tirosh et al. [4].

Cluster enrichment analysis for quantitative and categorical features

The association between cell clusters and any cell-level feature is performed through Cell Set Enrichment Analysis (CSEA) and Over Representation Analysis (ORA) for, respectively, quantitative and categorical values. Importantly, these two types of analysis are implemented in parallel, which is particularly important for CSEA that uses permutations to build an empirical null distribution.

CSEA is the application of Gene Set Enrichment Analysis (GSEA) [8] to cells and cell sets in place of, respectively, genes and gene sets: a ranked cell list is used instead of a ranked gene list, and cell sets (cell clusters) are tested instead of gene sets. In addition, the code is implemented to handle missing values. Therefore, scMuffin tests whether the cells assigned to a cluster are located at the top (bottom) of a ranked list of cells.

The assessment of cluster enrichment in a particular value of a categorical feature is computed using the over-representation analysis (ORA) approach [9], which is based on the hypergeometric test. Here, this test assesses whether the occurrence of a particular value between the cells of a cluster in relation to all other clusters is higher than expected in a hypergeometric experiment.

Overlap between partitions of cells

Cell partitions resulting from different clustering analyses are compared by calculating the overlap coefficients among all-pairs of clusters. Given two partitions A and B, defined as sets of cell clusters $A = \left\{ {a_{i} } \right\}{ }$ and $B = \left\{ {b_{j} } \right\}$, the similarity between the cell clusters $a_{i}$ and $b_{j}$ is calculated as:

$$o_{ij} = \frac{{\left| {a_{i} \cap b_{j} } \right|}}{{\min \left( {\left| {a_{i} } \right|,\left| {b_{j} } \right|} \right)}}$$

Results and discussion

In this section, we present the user interface (Table 1), taking advantage of the results obtained on the SC dataset generated by Yuan et al. [10] from human high-grade glioma (HGG) samples (available on the Gene Expression Omnibus (GEO) repository [11] with identifier GSE103224, see Supplementary Methods in Additional file 1).

Table 1 Main tasks and corresponding functions in scMuffin

Gene set scoring

scMuffin provides functions to set up one or more gene set collections and perform SC-level estimation of gene set expression scores in relation to an empirical null model (see Implementation section). This can be applied to any gene set and can therefore be used to estimate several different cell phenotypes, like pathway activities or marker set expression.

The function prepare_gsls() allows the user to collect gene sets of cell types, pathways, cancer functional states, as well as other collections of gene sets (e.g. positional gene sets, hallmarks) from CellMarker [12], PanglaoDB [13], CancerSEA [14] and MSigDB [8] databases. Unlike many existing tools that are used to perform marker-based cell annotation [15], the availability of these gene sets within scMuffin package spares the user the effort of data collection and harmonization. The function, which also accepts any user-given gene sets, applies a series of criteria (e.g., minimum and maximum number of genes in a gene set) to filter the chosen gene sets.

The cell-level expression scores for these gene sets can be calculated using calculate_gs_scores(). For instance, the following code shows how to quantify the activity of “Cancer functional states” from CancerSEA [14] at cell and cell cluster level:

gsc <- prepare_gsls(gs_sources = "CancerSEA", scMuffinList = scML)

scML <- calculate_gs_scores(scMuffinList = scML, gs_list = gsc$CancerSEA)

scML <- calculate_gs_scores_in_clusters(scMuffinList = scML, partition_id = "global_expr")

where scML is the data structure that contains expression data and cell clusters (“global_expr”).

Cluster-level gene sets scores can be visualized as a heatmap to obtain a useful summary visualization that can provide insights for cluster annotation. This is accomplished by means of the function plot_heatmap_features_by_clusters(), which relies on the powerful ComplexHeatmap R package [16]). For example, the analysis of the CancerSEA functional states in the HGG sample PJ016 showed that the two groups of cell clusters that are spatially separated in the UMAP visualization, that is clusters {0, 6, 8, 9} and {1, 2, 3, 4, 5, 7} (Fig. 2a), reflect distinct functional states (Fig. 2b). Cell-level gene set scores can be visualized over the UMAP by means of plot_umap_colored_features() (Fig. 2c).

Fig. 2

Quantification of CancerSEA functional states in the HGG sample PJ016. a UMAP visualization where cells are coloured by expression clusters. b Cluster-level expression scores (z-score) of all the CancerSEA functional states. c UMAP visualization where cells are colored by “Invasion” gene set score

CNV estimation and association with CancerSEA functional states

CNV inference from SC transcriptomics in cancer provides a means to assess the presence of relevant genomic aberrations (duplications and deletions) based on the expression of adjacent genes. This knowledge offers crucial clues to address the difficult task of distinguishing normal from malignant cells and provides quantitative information to reconstruct the tumor clonal substructure. Moreover, CNV patterns allow the investigator to hypothesize link between genomic alterations and cell phenotypes. As a proof-of-concept, we describe CNV inference in scMuffin with and without a reference expression profile, and the association of CNV clusters with CancerSEA functional states.

CNV inference is performed by means of the function CNV_analysis():

scML <- CNV_analysis(scMuffinList = scML, reference = GTEx_mean, center_genes = TRUE)

where GTEx_mean contains a reference expression profile, which can be obtained by means of the function process_GTEx_gene_reads() starting from data available at The Genotype-Tissue Expression (GTEx) portal [17]. The resulting scMuffinList scML contains the matrix with CNV values for each genomic regions in every cell, cell clusters by CNV, detected CNV regions and mapping information between gene locations and CNV regions.

Heatmaps that show CNV profiles of every cell, clustered by similarity, can be generated by means of plot_heatmap_CNV(). In the considered dataset (PJ030), we observed that the reference profile (normal brain samples from GTEx) is included into cluster 3, while clusters 0, 1 and 2 show large aneuploidies, some of which are typical of HGG, like the amplification of chromosome 7 as reported by dataset authors [10] (Table 2, Fig. 3a, b). The median CNV profile of every cell cluster can be visualized by means of plot_profile_CNV() (Fig. 3c) while boxplots of the distribution of CNV score by clusters can be visualized through boxplot_points() (Fig. 3d).

Table 2 CNV regions detected in chromosomes 7 and 19 of sample PJ030

Fig. 3

CNV analysis (PJ030). a CNV heatmap where cells (columns) are grouped into “CNV clusters”; the cluster of cells that include the reference profile is highlighted through a red label. b UMAP visualizations where cells are colored by CNV clusters. c z-score of cluster 0 median CNV profile. d Distribution of cell CNV score (C) in clusters

The analysis of another dataset (PJ016, without reference profile) showed two groups of clusters ({0, 1, 4} and {2, 3}) that map to separated regions of the UMAP visualization (Fig. 4a, b). Interestingly, clusters 2 and 3 are marked by peculiar amplifications in chromosome arms 1p and 19p. These “CNV clusters” correspond to expression clusters {0, 6, 8, 9} (Fig. 4c). The similarity of multiple cell partitions can be quantified using the function overlap_matrix().

Fig. 4

CNV analysis and CancerSEA functional states (PJ016). a CNV heatmap where cells (columns) are grouped into “CNV clusters”. b UMAP visualizations where cells are colored by CNV clusters. c Overlap indices between expression clusters (prefix “EXPR_”) and CNV clusters (prefix “CNV_”). d Heatmap of CancerSEA functional state expression score (z-score) in CNV clusters

An example of integrative analysis enabled by scMuffin is the functional assessment of CNV clusters. We quantified the expression scores of CancerSEA functional states throughout the CNV clusters of sample PJ016. As expected, the two aforementioned groups of CNV clusters {0, 1, 4} and {2, 3} are characterized by different functional states (Fig. 4d). This finding suggests that the peculiar amplifications of chromosomes 1 and 19 that characterize these clusters might underlie two phenotypes. A suggestive evidence of such hypothesis, with regards to the CancerSEA functional state “Invasion”, is the location of the two Invasion markers YBX1 (Y-Box Binding Protein 1) and HNRNPM (Heterogeneous Nuclear Ribonucleoprotein M) within the amplified regions of chromosomes 1p and 19p (Fig. 4a). YBX1 is a DNA/RNA-binding protein and transcription factor which plays a central role in coordinating tumor invasion in glioblastoma [18]. HNRNPM belongs to a family of spliceosome auxiliary factors and is involved in the regulation of splicing; the upregulation of these factors results in tumor-associated aberrant splicing, which promotes glioma progression and malignancy [19, 20]. In particular, HNRNPM was identified as an interactor of the DNA/RNA binding protein SON, which drives oncogenic RNA splicing in glioblastoma [21]. While it is beyond the scope of this article to validate this hypothesis, these findings clearly highlight the usefulness of the integrative analysis of CNVs and CancerSEA functional states provided by scMuffin.

Cluster enrichment analysis: association of clusters and features

scMuffin contains functions for assessing the association between cell clusters and quantitative as well as categorical features, by means of CSEA and ORA, respectively. The function assess_cluster_enrichment() performs CSEA and ORA accordingly to the feature type under consideration. For example, the following instruction runs CSEA over the gene set scores, considering the partition “global_expr”:

scML <- assess_cluster_enrichment(scMuffinList = scML, feature_name = "gene_set_scoring", partition_id = "global_expr")

The results of CSEA and ORA include several quantities, such as the normalized enrichment score (NES), the enrichment ratio (ER) and the false discovery rate (FDR) (Tables 3, 4). These results can be visualized as boxplots (quantitative features) and barplots (categorical features). For instance, the cluster 0 of sample PJ016 is particularly enriched—as expected—in “Invasion” markers (Table 3, Fig. 5a) and cells assigned to G2/M or S phases (Table 4, Fig. 5b). The function extract_cluster_enrichment_table() provides a means to obtain tables of enrichment results across clusters, which can be visualized as heatmaps (Fig. 5c), while extract_cluster_enrichment_tags() defines the most significant features as cluster labels that are therefore available for annotating UMAPs (Fig. 5d).

Table 3 The output of CSEA analysis for a quantitative feature

Table 4 The output of an ORA analysis for a categorical feature

Fig. 5

Cluster enrichment analysis (PJ016). a Output of boxplot_cluster() function for a quantitative feature—CancerSEA functional states—in cluster 0: distribution of gene set scores values (red) and corresponfing null models (gray, suffix “_ref”) (top); Normalized Enrichment Score (z-score) and FDR q-value (bottom). b Output of barplot_cluster() function for a categorical feature—cell cycle phase—in cluster 0: expected (“exp”) and observed (“obs”) occurrences of every possible category (top); enrichment ratio (ER) and hypergeometric p-value (bottom). c Heatmap of NES values for CancerSEA functional states in expression clusters; asterisks indicate FDR q-value < 0.01. d UMAP visualization where clusters are labelled with the top three most significant CancerSEA functional states

Transcriptional complexity, proliferation rate and cell state trajectories

The number of expressed genes and proliferation rate are two relevant pieces of information for the characterization of cell type and cell state in cancer.

Undifferentiated stem or progenitor cells show a high number of expressed genes (open chromatin state) compared to differentiated cell types. In cancer context, a high number of expressed genes might suggest tumor initiating cells (TICs) or might indicate de-differentiation processes of tumor progression [22, 23]. For instance, in a recent study on glioblastoma, chromatin accessibility was associated with a specific subpopulation of putative tumor-initiating Cancer Stem Cells (CSCs) with invasive phenotype and low survival rate prediction [24]. The transcriptional complexity of a cell can be calculated as follows:

scML <- transcr_compl(scMuffinList = scML, min_counts=3, min_cells=10, min_genes=500)

where min_counts is the threshold above which a gene is considered expressed, min_cells is the minimum number of cells in which a gene must be expressed and min_genes is the minimum number of genes that a cell must express. In SC datasets, transcriptional complexity quantification is complicated by the presence of a linear relationship between number of expressed genes and total transcripts detected, where biological variability is entangled with sequencing deepness. scMuffin provides three measures of transcriptional complexity, which show three different – but complementary—relations with the number of total transcripts (t): the ratio (TC-R) prioritizes cells with a relatively low t (Fig. 6a); the linear model residuals (TC-LMR) are independent from t (Fig. 6b); lastly, the transcriptional entropy (TC-H) [7] highlights genes with relatively high number of t (Fig. 6c). All three measures may provide useful insights for the identification of particular cell types or states. For example, TC-R could be used when one is not interested in cell with a relatively abundant RNA production, like can be the case of quiescent cells with stem-like properties.

Fig. 6

Transcriptional complexity, proliferation rate and cell state trajectories (PJ016). a Transcriptional Complexity Ratio (TC-R) of cells. b Transcriptional Complexity Linear Model Residuals (TC-LMR) of cells. c Transcriptional Complexity Entropy (TC-H) of cells. d Distribution of TC-R values split by cell proliferation score e–f Distribution of cells according to the first two “diffusion components” (DC)

The proliferation rate is a relevant indicator for distinguishing cell types in solid tumors and helps to identify cells with potential clinical relevance and interest as candidate therapeutic targets [25, 26]. Cell proliferation rate is quantified based on the expression of G1/S and G2/M genes:

scML <- proliferation_analysis(scMuffinList = scML).

As a proof-of-concept we jointly analysed cell TC and cell proliferation of sample PJ016 and visualized the results on a diffusion map [27]—which can be run by means of the wrapper function diff_map(). This joint analysis shows that cells with the highest values of TC-R tend to have lower proliferation scores (Fig. 6d). These cells are located more closely to the branching point of the diffusion map (Fig. 6e), while the cells with the highest proliferation rates are located at a corner of one of the branches (Fig. 6f). This evidence suggests that the cells with high TC and low proliferation are interesting candidates for further analyses aimed at finding TICs in this dataset of HGG. Interestingly, this particular pattern is captured by TC-R, but not by TC-LMR and TC-H (Additional file 1: Fig. S1).

Comparison of datasets

A SC dataset carries an extensive amount of information. The integration of multiple SC datasets is a challenging task and multiple approaches have been proposed to address it [28]. Typically, the integrated datasets are computationally demanding due to their huge size. An alternative possibility lies in cross-checking the expression of cluster markers between datasets: the expression of the cluster markers of every dataset is assessed in all considered datasets. For example, Nguyen et al. [5] used this approach to study the occurrence of the characteristic cell types of normal mammary gland across samples collected from different subjects.

scMuffin provides a function to quantify the expression of a series of gene sets across multiple datasets:

res <- inter_dataset_comparison(seu_obj_list = seu_obj_list, gsl = cluster_markers_list)

where seu_obj_list is a list of Seurat objects (each containing a SC dataset). If the gene sets are defined as cluster markers (“cluster_markers_list” in the code above), then the analysis sheds light over the possible presence of cells with similar expression profiles across datasets.

We performed this analysis on three datasets (PJ016, PJ017, PJ030) and found cell clusters that are more similar to clusters of other datasets than clusters of the same dataset, suggesting the presence of similar cells (Fig. 7a), like cluster 7 of PJ030 (“PJ030_C7”), which is close to clusters {0, 6, 8, 9} of PJ016. Importantly, the same analysis can be done using gene sets with a functional meaning, to obtain a functional characterization of similar cell types across datasets (Fig. 7b).

Fig. 7

Inter-dataset comparison of gene set expression. a Heatmap of cluster marker (columns, “_Mx” suffix, where “x” indicates cluster number) expression scores across clusters (rows, “_Cx” suffix, where “x” indicates cluster number) of datasets PJ016, PJ017 and PJ030. b Heatmap of CancerSEA functional states scores across clusters of datasets as in panel (a)

scMuffin and scCancer

scCancer [29] is an R package for automated processing of SC expression data in cancer. The two packages (scMuffin and scCancer) address the same challenge and are complementary. Indeed, each of them offers some analyses that the other does not (Additional file 1: Table S1). Examples of this are cluster association analysis (CSEA and ORA), comparison of multiple cell partitions and inter-dataset gene set expression assessment, available in scMuffin, and cell interaction, expression programs and cancer micro-environment cell classification, provided by scCancer. Other types of analyses are common to both packages but differ by the underlying algorithm and/or implementation. This is the case, for instance, of gene set scoring and CNV analysis.

Conclusions

Here, we presented scMuffin, an R package that provides a series of functions to perform and integrate various types of analyses on SC expression data. As a proof-of-concept, to illustrate the user interface and its potentialities, we analysed publicly available SC expression datasets of human HGG. We described two examples of integrative analyses which returned particularly interesting findings that would deserve further investigations: the functional characterization of CNVs highlighted a possible link between amplifications of chromosomes 1p and 19p and invasive tumor phenotype; the joint analysis of transcriptional complexity, proliferation rate and cell state trajectories spotted cells that have some characteristics of TICs. Despite the inference of CNV from SC expression data is a challenging task, subject to various issues (e.g., missing data, number of aberrancies), it provides useful insights to clarify cell identity, like the distinction between a more aberrant cluster and the others or the distinction of tumoral and normal cells. The functions of scMuffin can be combined conveniently to address various tasks. For example, different partitions of the same dataset obtained using the same clustering method with different parameters can be characterized by means of their overlap and gene set expression, to gain insights on the appropriate number of clusters. The analyses offered by scMuffin and the results achieved in the proof-of-concept show that scMuffin helps addressing the main challenges in the bioinformatics analysis of SC datasets from solid tumors.