Bayesian segmented Gaussian copula factor model for single-cell sequencing data

Let $\bm{x}=(x_{1},\ldots,x_{p})^{\top}$ denote $p$ random variables of interest, i.e., read counts of an individual cell in single-cell sequencing data. In order to account for an excessive number of low counts, including zeros, we propose the segmented Gaussian copula model. We begin by reviewing the standard Gaussian copula model developed by Liu et al. (2009).

Our proposed segmented Gaussian copula model accounts for dropouts (low counts) by associating them with segments of underlying latent continuous values. Let $m\in\{0,1,\ldots\}$ be a hyperparameter determining the maximum low count being inflated in the observed data. Figure 2 illustrates a schematic of the proposed segmented Gaussian copula model when $m=1$, and the formal definition for any $m$ follows below.

As illustrated in Figure 2 and formally introduced in Definiton 2, the observed variables $\bm{x}$ are generated through a two-stage process. Initially, a latent random vector $\bm{x}^{*}$ is generated from the Gaussian copula model, serving as an underlying source for the observed data. Subsequently, the observed data $\bm{x}$ are derived from the latent random vector $\bm{x}^{*}$ via a segmentation process. Specifically, the latent random variable $x_{j}^{*}$ is segmented into multiple intervals $(c_{jd},c_{j,d+1}]$, $d=0,\ldots,m$, and $(c_{j,m+1},\infty)$. If $x_{j}^{*}$ exceeds a specific threshold $c_{j,m+1}$, corresponding to the interval $(c_{j,m+1},\infty)$, the latent source $x_{j}^{*}$ is directly observed (i.e., $x_{j}=x_{j}^{*})$. Otherwise, if $x_{j}^{*}\in(c_{jd},c_{j,d+1}]$ where $d\in\{0,1,\ldots,m\}$, it results in the observation of a low count value $x_{j}=d$. Hence, $\bm{c}=(c_{j0},\ldots,c_{j,m+1})$ represents unknown thresholds that truncate the latent random variable $x_{j}^{*}$ to near-zero counts for the observed variable $x_{j}$. Additionally, given that the transformations $f_{j}$ are strictly increasing, (1) can be equivalently expressed using the latent Gaussian vector $\bm{z}=(f_{1}(x_{1}^{*}),\ldots f_{p}(x_{p}^{*}))^{\top}$,

where $\delta_{jd}=f_{j}(c_{jd})$ is the latent Gaussian level threshold corresponding to $c_{jd}$.

When $m=0$, the proposed segmented Gaussian copula model is reduced to the truncated Gaussian copula model of Yoon et al. (2020), which focuses on addressing zero inflation. When $m>0$, the segmented Gaussian copula model is capable of accommodating the inflation of low counts besides zero. For example, if $m=1$, the latent variable $x_{j}^{*}$ is segmented into three intervals $(c_{j0},c_{j1}],(c_{j1},c_{j2}]$, and $(c_{j2},\infty)$. If $x_{j}^{*}$ falls into either one of the first two intervals $(c_{j0},c_{j1}]$ and $(c_{j1},c_{j2}]$, our observation will be either $x_{j}=0$ or $1$. This allows for a flexible representation of zero and one counts in the observed data.

Next, we extend the proposed segmented Gaussian copula to scFM to capture underlying dependencies among observed genes in single-cell sequencing data while simultaneously addressing the issue of inflated near-zero counts. We adopt a low-rank decomposition of the copula correlation matrix $\bm{\Omega}$ so that the dependence among the latent variables can be explained by low-dimensional latent factors.

Note that the covariance matrix $\bm{\Omega}$ of the latent Gaussian $\bm{z}$ is constrained to be a correlation matrix for identifiability (Liu et al., 2009), and thus $\bm{\Omega}=\bm{\Psi}^{-1/2}(\bm{\Lambda}\bm{\Lambda}^{\top}+\bm{\Sigma})\bm{\Psi}^{-1/2}$ implies $\bm{\Psi}=\mbox{diag}\{\bm{\lambda}_{j\cdot}^{\top}\bm{\lambda}_{j\cdot}+\sigma_{j}^{2}\}_{j=1}^{p}$, where $\bm{\lambda}_{j\cdot}$ is the $j$-th row of the factor loading matrix $\bm{\Lambda}$. In practice, it is typically assumed that the number of latent factors $k$ is much smaller than the number of observed variables $p$, leading to a low-rank characterization of $\bm{\Omega}$. Therefore, once the factor scores $\bm{U}$ are estimated, the dimensionality of the dataset can be effectively reduced from $p$ to a lower dimension $k$ $(\ll p)$, which facilitates downstream analyses. A fundamental challenge here is that the number of latent factors $k$ is typically unknown a priori. Hence, the selection of an appropriate $k$ is of great significance. In the next subsection, we introduce a meticulously selected prior distribution for the factor loadings that enables the automatic selection of the number of latent factors as well as facilitates the identifiability of latent factors.

We propose to use a column-wise Dirichlet-Laplace (DL, Bhattacharya et al. 2015) prior for the factor loading matrix $\bm{\Lambda}=(\lambda_{jh})$. The column-wise DL prior enables simultaneous shrinkage of all the entries of a column, thereby effectively zeroing out unnecessary latent factors. Specifically, we assume, for $j=1,\ldots,p$ and $h=1,\ldots,k_{\text{max}}$,

where $\bm{\phi}=(\phi_{1},\ldots,\phi_{k_{\text{max}}})^{\top}$, $\mbox{DE}(\tau\phi_{h})$ denotes the Laplace or double-exponential distribution with scale parameter $\tau\phi_{h}$, $\mbox{Dir}(\alpha,\ldots,\alpha)$ denotes the Dirichlet distribution with concentration parameter $\alpha$ which is considered as a hyperparameter.

The global scale $\tau$ controls the overall shrinkage towards zero, while the column-specific local scale $\phi_{h}$ allows each column of the factor loading matrix $\bm{\Lambda}$ to deviate away from zero. Here $k_{\text{max}}$ represents the maximum allowable number of latent factors. Since DL prior allows for aggressive shrinkage of all entries of a column towards zero, the effective number of factors can be much smaller than $k_{\text{max}}$, and is determined automatically; see Section 3.3.

Intuitively, the class of column-wise DL priors approximates a spike-and-slab prior (George and McCulloch, 1997) assigned to the factor loadings at the column level by a continuous density concentrated near zero with a heavy tail. A distinguishing characteristic of the DL prior, in contrast to commonly used shrinkage priors like Laplace priors (Park and Casella, 2008), is the presence of a singularity at zero in its distribution while still maintaining exponential tails. Under the proposed column-wise DL priors, the singularity at zero facilitates the concentration of most of the columns of $\bm{\Lambda}$ around zero, while the exponential tails ensure little shrinkage of a small number of significant factors. See Bhattacharya et al. (2015) for the theoretical properties of the DL prior.

In addition to the adaptive shrinkage, imposing DL priors for the factor loadings also resolves the identifiability issue of factor models. In the Bayesian framework, a typical approach for factor analysis is assigning Gaussian priors to the factor loadings (Lopes and West, 2004). Although this approach allows for a conditionally conjugate model, facilitating posterior inference through standard Gibbs sampling, it fails to address the invariance of factor models to orthogonal transformations as $\bm{\Lambda}\bm{u}=(\bm{\Lambda}\bm{P})(\bm{P}^{\top}\bm{u}):=\bm{\Lambda}^{*}\bm{u}^{*}$ for any orthogonal matrix $\bm{P}$ (Anderson and Rubin, 1955). A popular technique to deal with this invariance property is imposing conditions that remove invariance to orthogonal transformation, e.g., the lower-triangular constraint on $\bm{\Lambda}$ with an assumption of strictly positive diagonal elements (Geweke and Zhou, 1995; Aguilar and West, 2000). While this approach is effective for estimating the covariance matrix, the constraint introduces order-dependence in estimating factor loadings, making the interpretation of latent factors sensitive to the ordering of observed variables. Consequently, changes in variable ordering can influence downstream analyses, such as clustering, thereby compromising the reproducibility of their results. In contrast, our approach, which assigns the column-wise DL priors to the factor loadings, breaks this invariance in the posterior distribution without necessitating the lower-triangular constraint on the factor loadings matrix. Our factor model (3) with non-Gaussian DL prior can be viewed as a noisy independent component analysis (ICA) model, where the columns of $\bm{\Lambda}$ are the hidden sources, $\bm{u}$ is the mixing matrix, and $\bm{\epsilon}$ is the noise. Since the sources $\bm{\Lambda}$ are non-Gaussian, the theory of noisy ICA model ensures identifiability of our model parameters $\bm{\Lambda}$ and $\bm{u}$ (up to column permutations and sign changes) (Comon, 1994; Hyvärinen et al., 2004). In summary, our approach alleviates the burden of subjective decisions regarding the ordering of observed variables, which often introduces ambiguity in conventional factor models. As a result, our approach achieves greater consistency in the interpretation of latent factors, and contributes to the robustness and reliability of downstream analyses.

For the Gaussian level thresholds $\bm{\delta}=(\delta_{j0},\ldots,\delta_{j,d})^{\top}$ in (2), we assume an improper uniform prior with the following monotonicity constraint,

For the error variances $\sigma_{j}^{2}$, we assign a conditionally conjugate inverse-gamma prior $\sigma_{j}^{2}\sim\textup{IG}(a_{\sigma},b_{\sigma})$.

The latent Gaussian variables corresponding to the observed data $\bm{x}_{i}^{o}$ are defined by $z_{ij}^{o}=f_{j}(x_{ij}^{o})$. We will first compute “observed pseudodata” $\hat{z}_{ij}^{o}$ by estimating $f_{j}$ via the empirical cumulative distribution function (cdf).

Let $F_{j}$ be the (monotonically increasing) marginal cdf of the $j$-th latent random variable $x_{j}^{*}$ and $\Phi$ be the cdf of standard Gaussian. Then, we observe

which implies that $f_{j}=\Phi^{-1}\circ F_{j}$. Since $F_{j}$ is unknown, we replace it with its empirical estimate (Klaassen and Wellner, 1997),

Here, the constant term $n/(n+1)$ is necessary to restrict $\hat{F}_{j}(x)<1$ to make $\hat{z}_{ij}^{o}$ finite. Given $\hat{F}_{j}$, we compute the observed pseudodata $\hat{z}_{ij}^{o}=\hat{f}_{j}(x_{ij}^{o})=\Phi^{-1}\circ\hat{F}_{j}(x_{ij}^{o})$.

Given the observed pseudodata, we employ an MCMC algorithm to draw posterior samples of the model parameters. Let $\hat{\bm{Z}}^{o}=\{\hat{\bm{z}}_{i}^{o}\}_{i=1}^{n},\bm{Z}^{\leq 1}=\{(\bm{z}_{i}^{0},\bm{z}_{i}^{1})\}_{i=1}^{n}$, and $\bm{U}=\{\bm{u}_{i}\}_{i=1}^{n}$. Given $\hat{\bm{Z}}^{o}$, we use MCMC to draw posterior samples of $\bm{Z}^{\leq 1},\bm{\delta},\bm{\Sigma},\bm{U},\bm{\Lambda},\bm{\phi}$, and $\tau$. In particular, we utilize an efficient data-augmented Gibbs sampler, which is derived by a data-augmented representation of the proposed column-wise DL priors. Introducing auxiliary variables $\bm{\Xi}=(\xi_{jh})_{j,h}$, the DL prior (2.2) can be equivalently represented as

which facilitates Gibbs sampling by enabling the derivation of full conditionals for all parameters. Our implementation of the Gibbs sampler leverages marginalization and blocking to reduce auto-correlation in the posterior samples. It is given that, upon conditioning on the latent Gaussian vectors $\bm{Z}$, the thresholds $\bm{\delta}$ become independent of the other parameters, and similarly, conditioning on the factor loadings $\bm{\Lambda}$, the hyperparameters $(\bm{\Xi},\tau,\bm{\phi})$ are independent from the remaining parameters. Utilizing these conditional independences, our Gibbs sampler cycles through the following draws: (i) $\bm{Z}^{\leq 1}|\bm{\delta},\bm{\Sigma},\bm{U},\bm{\Lambda},\hat{\bm{Z}}^{o}$, (ii) $\bm{\delta}|\bm{Z}^{\leq 1},\hat{\bm{Z}}^{o}$, (iii) $\bm{\Sigma}|\bm{U},\bm{\Lambda},\bm{Z}^{\leq 1},\hat{\bm{Z}}^{o}$, (iv) $\bm{U}|\bm{\Sigma},\bm{\Lambda},\bm{Z}^{\leq 1},\hat{\bm{Z}}^{o}$, (v) $\bm{\Lambda}|\bm{\Sigma},\bm{U},\bm{\Xi},\tau,\bm{\phi},\bm{Z}^{\leq 1},\hat{\bm{Z}}^{o}$, and (vi) $\bm{\Xi},\tau,\bm{\phi}|\bm{\Lambda}$.

Sample $\bm{Z}^{\leq 1}|\bm{\delta},\bm{\Sigma},\bm{U},\bm{\Lambda},\hat{\bm{Z}}^{o}$. We draw $(\bm{z}_{i}^{0},\bm{z}_{i}^{1})$ from a truncated multivariate Gaussian distribution,

where $\bm{\mu}_{i}$ and $\bm{C}_{i}$ are the conditional mean and conditional covariance matrix of $(\bm{z}_{i}^{0},\bm{z}_{i}^{1})$ given $\bm{z}_{i}^{o}=\hat{\bm{z}}_{i}^{o}$ and $(\bm{\Sigma},\bm{U},\bm{\Lambda})$. Since $\mbox{cov}(z_{ij},z_{ij^{\prime}}|\bm{\Sigma},\bm{U},\bm{\Lambda})=0$ for $j\neq j^{\prime}$, $\bm{\mu}_{i}$ and $\bm{C}_{i}$ are the sub-vector and sub-matrix of $\mbox{E}(\bm{z}_{i}|\bm{\Sigma},\bm{U},\bm{\Lambda})=\bm{\Psi}^{-1/2}\bm{\Lambda}\bm{u}_{i}$ and $\mbox{Var}(\bm{z}_{i}|\bm{\Sigma},\bm{U},\bm{\Lambda})=\bm{\Psi}^{-1/2}\bm{\Sigma}\bm{\Psi}^{-1/2}$ corresponding to the variables $(\bm{z}_{i}^{0},\bm{z}_{i}^{1})$, respectively. In (6), the vector $\bm{\delta}_{i,d}^{a}=\{\delta_{jd}:x_{ij}=a,j=1,\ldots,p\}$ consists of the thresholds $\delta_{jd}$ for the $i$-th observation, specifically corresponding to those variables $x_{ij}$ that equal $a$. The subscript $d$ signifies whether the threshold represents a lower bound $(d=a)$ or an upper bound $(d=a+1)$, respectively.

Sample $\bm{\delta}|\bm{Z}^{\leq 1},\hat{\bm{Z}}^{o}$. The conditional posterior of $\bm{\delta}|\bm{Z}^{\leq 1},\hat{\bm{Z}}^{o}$ is proportional to the $n$-product of indicator functions in (6), and thus Step (ii) can be done by sampling $(\delta_{j1},\delta_{j2}),j=1,\ldots,p$, independently from the uniform full conditional distribution

Sample $\bm{\Sigma}|\bm{U},\bm{\Lambda},\bm{Z}^{\leq 1},\hat{\bm{Z}}^{o}$. We draw $\sigma_{j}^{2}$ from $\sigma_{j}^{2}\sim\textup{IG}(a_{\sigma^{2}}+n/2,b_{\sigma^{2}}+\sum_{i=1}^{n}(z_{ij}-\bm{\lambda}_{j}^{\top}\bm{u}_{i}))$.

Sample $\bm{U}|\bm{\Sigma},\bm{\Lambda},\bm{Z}^{\leq 1},\hat{\bm{Z}}^{o}$. We draw $\bm{u}_{i}\sim\textup{N}(\bm{m}_{i},\bm{V}_{i})$ where $\bm{m}_{i}=\bm{V}_{i}(\bm{\Lambda}^{\top}\bm{\Sigma}^{-1}\bm{z}_{i})$ and $\bm{V}_{i}=(\bm{\Lambda}^{\top}\bm{\Sigma}^{-1}\bm{\Lambda}+\bm{I}_{k_{\text{max}}})^{-1}$.

Sample $\bm{\Lambda}|\bm{\Sigma},\bm{U},\bm{\Xi},\tau,\bm{\phi},\bm{Z}^{\leq 1},\hat{\bm{Z}}^{o}$. We draw each row of $\bm{\Lambda}$ from $\bm{\lambda}_{j}\sim\textup{N}(\bm{\eta}_{j},\bm{W}_{j})$, where $\bm{\eta}_{j}=\bm{W}_{j}(\sigma_{j}^{-2}\sum_{i=1}^{n}z_{ij}\bm{u}_{i})$ and $\bm{W}_{j}=(\sigma_{j}^{-2}\sum_{i=1}^{n}\bm{u}_{i}\bm{u}_{i}^{\top}+\bm{D}^{-1})^{-1}$ with $\bm{D}=\mbox{diag}\{\xi_{jh}\tau^{2}\phi_{jh}^{1}\}_{h=1}^{k_{\text{max}}}$.

Sample $\bm{\Xi},\tau,\bm{\phi}|\bm{\Lambda}$. Since

we sequentially sample from $\bm{\phi}|\bm{\Lambda}$, $\tau|\bm{\phi},\bm{\Lambda}$, and $\bm{\Xi}|\tau,\bm{\phi},\bm{\Lambda}$. For sampling $\bm{\phi}|\bm{\Lambda}$, we follow the strategy in Bhattacharya et al. (2015). Specifically, let $\textup{giG}(\kappa,\rho,\chi)$ denote a three-parameter generalized inverse Gaussian distribution with the density $p(y)\propto y^{\kappa-1}e^{-0.5(\rho y+\chi/y)}$ for $y>0$. We draw $R_{11},\ldots,R_{pk_{\text{max}}}$ independently, $R_{jh}\sim\mbox{giG}(\alpha-1,1,2|\lambda_{jh}|)$, and set $\phi_{jh}=R_{jh}/R$ with $R=\sum_{j,h}R_{jh}$. Given $\bm{\phi}$, we then sample $\tau|\bm{\phi},\bm{\Lambda}$ from $\mbox{giG}(pk_{\text{max}}(\alpha-1),1,2\sum_{j,h}|\lambda_{jh}|/\phi_{jh})$. Lastly, we sample $\bm{\Xi}|\bm{\phi},\tau,\bm{\theta}$ by drawing $\xi_{jh}$ independently from $\mbox{giG}(1/2,1,\lambda_{jh}^{2}/(\tau^{2}\phi_{jh}^{2}))$. An outline of our data-augmented Gibbs sampler is provided in Algorithm 1.

In the proposed scFM, $k_{\text{max}}$ is the maximum allowable number of latent factors, and the number of significant factors, denoted by $\hat{k}\leq k_{\text{max}}$, is determined by the proposed column-wise DL prior. The column-wise DL prior induces a shrinkage effect on factor loadings associated with irrelevant latent factors that lack support from the data, thereby aiding in the identification of $\hat{k}$. More specifically, our DL prior promotes the concentration of the majority of columns in the factor loadings matrix $\bm{\Lambda}$ around zero, except for a few columns corresponding to a limited number of significant factors. This implies the existence of potentially two clusters of columns in $\bm{\Lambda}$: one cluster closely centered around zero, indicating insignificant factor, and another cluster distant from zero, indicating significant factors.

Consequently, we suggest a simple automated approach for selecting $\hat{k}$. At each MCMC iteration, we cluster the columns of the factor loadings matrix $\bm{\Lambda}$ into two clusters using k-means, and estimate the number of significant factors by the size of the cluster that exhibits greater deviance from zero. The number of significant factors $\hat{k}$ is then estimated by taking the mode over all the MCMC iterations. Upon obtaing the estimated $\hat{k}$, we identify the significant factors by selecting the latent factors associated with the $\hat{k}$ largest columns (based on $\ell_{2}$-norm) of the posterior mean of $\bm{\Lambda}$.

We apply the proposed scFM to our motivating scRNA-seq dataset, which consists of 7,247 cells from the lymphoblastoid cell line (LCL) GM12878, for which 18,170 genes were sequenced. The data can be obtained from the 10x Genomics Datasets website (https://www.10xgenomics.com). The goal of this analysis is to characterize underlying relationships among the sequenced genes and to obtain low-dimensional factor scores summarizing the variability in the observed scRNA-seq data to aid detection of previously uncharacterized cell subtypes in LCL. We filter out low-quality cells such as empty droplets, cell doublets and multiplets using the R package Seurat (Hao et al., 2023), resulting in $n=5,135$ cells. The detailed data preprocessing procedure is provided in Appendix A. We focus on $p=100$ genes that show the largest variability across cells, which contain $\sim 58\%$ zeros and $\sim 15\%$ ones across those genes.

For the proposed scFM, we use $m=1$ to explicitly account for the notable inflation of zero and one counts observed in the motivating scRNA-seq dataset (refer to Figure 1). We set the maximum number of latent factors at $k_{\text{max}}=8$ to achieve a balance between model interpretability and avoiding model misspecification, as fewer latent factors tend to provide clearer insights for model interpretation while too few risks excluding significant factors. Our analysis supports this decision, demonstrating that the number of significant latent factors is indeed fewer than $k_{\text{max}}=8$. We use the same hyperparameter values as in Section 4 and run the proposed MCMC in Section 3 for 10,000 iterations with 5,000 burn-in. To assess the goodness-of-fit of the proposed scFM on this dataset, we perform posterior predictive checks, which compare the posterior predictive distribution with the observed data. Our proposed model shows no significant lack of fit; see Appendix A for details.

We utilize the approach outlined in Section 3.3 to determine the number of significant factors $\hat{k}$ and to discern their identity. For this dataset, scFM determines that there are two significant factors (i.e. $\hat{k}=2$), and factors corresponding to the two largest columns (in $\ell_{2}$-norm) of the posterior mean of the factor loadings matrix $\bm{\Lambda}$ are identified as the significant factors (termed factors 1 and 2).

Figure 3 illustrates the absolute loadings of all genes on these two significant factors, providing insight into a potential biological interpretation for each estimated factor. For factor 1, genes with the 10 largest absolute loadings are MKI67, CENPF, HMGB2, NUSAP1, TOP2A, CDC20, CENPA, CKS1B, UBE2C, and STMN1. Notably, several of these genes are closely linked to the cell cycle, a fundamental process where cells go through a series of events leading to cell division. Specifically, MKI67, CENPF, HMGB2, and TOP2A have been known as cell cycle markers in single-cell studies of human lymphoblastoid cell lines (Sawada et al., 2020; SoRelle et al., 2022). Therefore, we call factor 1 the cell cycle factor.

For factor 2, genes with the ten most substantial loadings are MIR155HG, FSCN1, TRAF1, EBI3, IL4I1, LINC00158, BCL2A1, CD83, MARCKSL1, and IER3, some of which are known to be related to immune functions. FSCN1 and CD83 are markers for dendritic cells, specialized immune cells crucial for presenting antigens to other immune cells and initiating immune responses (Zimmer et al., 2012; Lechmann et al., 2002). Moreover, MIR155HG and IL4I1 are frequently over-expressed in immune cells including B-cells, and play important roles in immune infiltration (Peng et al., 2019; Bod et al., 2018). Therefore, we call factor 2 the immune cell factor.

To further interpret the two major latent factors identified by the proposed scFM, we use them to identify distinct cellular subtypes or states in LCL. The scores of the two selected factors serve as a low-dimensional representation of the high-dimensional gene expression measurements in scRNA-seq data. Therefore, by clustering the data in the low-dimensional factor space, we can identify groups of cells exhibiting similar expression behavior. For simplicity, we use the k-means algorithm to cluster the latent factor scores, where the number of clusters is determined by the elbow method.

Our analysis reveals three clusters, which are depicted in Figure 4. For cluster 1, we observe a notable elevation in scores for factor 1, while scores for factor 2 are close to zero. Cluster 2, in contrast, exhibits elevated scores for factor 2 but with diminished scores for factor 1. Cluster 3 displays factor scores around zero for both factors. Combining this observation with the biological interpretation of the estimated factors in Section 5.2 suggests that cluster 1 likely comprises cycling cells, cluster 2 likely consists of immune cells, and the remaining cells are grouped into cluster 3. We further investigate gene expression measurements for relevant genes discussed in Section 5.2 to support this interpretation. Table 2 reports average gene expression measurements of MKI67, CENPF, HMGB2, TOP2A, FSCN1, CD83, MIR155HG, and IL4I1 for the entire dataset and for each cluster. Recall that the first four of these genes are associated with factor 1, while the latter four are associated with factor 2. For cluster 1, MKI67, CENPF, HMGB2, and TOP2A are over-expressed. Given that MKI67, CENPF, HMGB2, and TOP2A serve as markers for cycling cells, it can be inferred that cells in cluster 1 are currently undergoing cell cycling. On the other hand, cluster 2 exhibits high expression of the latter four genes FSCN1, CD83, MIR155HG, and IL4I1, which is known to be strongly correlated with immune cells and immune cell function. Consequently, it is plausible to conclude that cluster 2 represents immune cells.