# RZiMM-scRNA: A regularized zero-inflated mixture model framework for single-
cell RNA-seq data

Xinlei<EMAIL_ADDRESS>[ William Bekermanlabel=e2 [
<EMAIL_ADDRESS>Peter A. Simslabel=e5 [
<EMAIL_ADDRESS>Peter D. Canolllabel=e3 [
<EMAIL_ADDRESS>Jianhua Hulabel=e4 [
<EMAIL_ADDRESS>Department of Preventive Medicine, Northwestern
University Feinberg School of Medicine,  Department of Statistics and Data
Science, Cornell University,  Department of Systems Biology, Columbia
University,  Department of Pathology and Cell Biology, Columbia University,
Department of Biostatistics, Columbia University,

###### Abstract

Applications of single-cell RNA sequencing in various biomedical research
areas have been blooming. This new technology provides unprecedented
opportunities to study disease heterogeneity at the cellular level. However,
unique characteristics of scRNA-seq data, including large dimensionality, high
dropout rates, and possibly batch effects, bring great difficulty into the
analysis of such data. Not appropriately addressing these issues obstructs
true scientific discovery. Herein, we propose a unified Regularized Zero-
inflated Mixture Model framework designed for scRNA-seq data (RZiMM-scRNA) to
simultaneously detect cell subgroups and identify gene differential expression
based on a developed importance score, accounting for both dropouts and batch
effects. We conduct extensive simulation studies in which we evaluate the
performance of RZiMM-scRNA and compare it with several popular methods,
including Seurat, SC3, K-Means, and Hierarchical Clustering. Simulation
results show that RZiMM-scRNA demonstrates superior clustering performance and
enhanced biomarker detection accuracy compared to alternative methods,
especially when cell subgroups are less distinct, verifying the robustness of
our method.

Our empirical investigations focus on two brain tumor studies dealing with
astrocytoma of various grades, including the most malignant of all brain
tumors, glioblastoma multiforme (GBM). Our goal is to delineate cell
heterogeneity and identify driving biomarkers associated with these tumors.
Notably, RZiMM-scNRA successfully identifies a small group of oligodendrocyte
cells which has drawn much attention in biomedical literature on brain
cancers. In addition, our method discovers several new biomarkers which are
not discussed in the original studies, including PLP1, BCAN, and PTPRZ1 — all
associated with the development and malignant growth of glioma — as well as
CAMK2B, which is downregulated in glioma and GBM, and implicated in
neurodevelopment, brain function, learning and memory processes.

Single cell,

Mixture model,

Clustering,

Batch effect,

Dropout,

###### keywords:

, and

##  1 Introduction

Single-cell RNA sequencing (scRNA-seq) techniques have recently emerged and
enabled large volumes of transcriptome profiling of individual cells (Wang and
Bodovitz, 2010; Prakadan, Shalek and Weitz, 2017). Unlike conventional gene
analyses of bulk samples where an average of gene expression levels in a bulk
cell population is measured and individual cell characteristics are ignored,
scRNA-seq data provide detailed information of gene expressions at the level
of individual cells (Hedlund and Deng, 2018). Research using scRNA-seq has led
to detection of new cell types, delineation of cell heterogeneity,
identification of gene regulatory mechanisms, cell developmental dynamics, and
other important discoveries (Xue et al., 2013; Treutlein et al., 2014; Grün et
al., 2015; Petropoulos et al., 2016; Reinius et al., 2016).

Our empirical investigations aim to advance understanding of brain tumors,
including glioblastoma multiforme (GBM), which is fast-growing and the most
aggressive of all brain tumors. The first study we encounter focuses on
various grades of astrocytoma, the most common type of primary brain tumors
affecting adults and can occur in most parts of the brain and in the spinal
cord. These tumors are quite resistant to therapies and often progress rapidly
to high-grade lesions (de Sousa et al., 2017). The astrocytoma scRNA-seq data
is generated from ten IDH-mutant human astrocytoma tumors of different grades
by Venteicher et al., and available in the Gene Expression Omnibus (GEO)
database under accession number GSE89567 (Barrett et al., 2012). In the second
study, high-grade astrocytoma scRNA-seq count data is generated by Yuan et al.
and available in the GEO database under accession number GSE103224. In this
study, samples from eight patients were taken of the cellular milieu of high-
grade gliomas, namely the most aggressive types — GBM and anaplastic
astrocytoma. We intend to identify hidden cell subgroup structures and
discover important biomarkers which drive these structures.

Despite its great promise in advancing our knowledge of diseases, several
challenges are encountered in the analysis of such scRNA-seq data, prompting
the development of novel analytical strategies (Kiselev, Andrews and Hemberg,
2019; Petegrosso, Li and Kuang, 2019). Following quality control and
normalization, an essential aspect of scRNA-seq data analysis is unsupervised
clustering to identify hidden cell subgroup structures and the biomarkers
which drive these patterns. However, many problems arise when clustering
scRNA-seq data (Kiselev, Andrews and Hemberg, 2019), such as high
dimensionality, high rate of dropouts, and the presence of sample/batch
effects. High dimensionality refers to the large numbers of cells and genes in
a tissue sample. Dropout events are zero expressions of transcripts in single
cells, due to low expressions/undetectable levels, false quantification, or
failure in amplification (Stegle, Teichmann and Marioni, 2015; Kiselev,
Andrews and Hemberg, 2019; Petegrosso, Li and Kuang, 2019), which can not be
addressed via simple normalization (Vallejos et al., 2017). Sample/batch
effects may be induced by sample heterogeneity, or technical issues and
experimental designs (Büttner et al., 2017). Failure to appropriately account
for these effects can introduce non-trivial bias and significantly impact
clustering algorithms (Tung et al., 2017).

On the other hand, many methods have been developed to particularly deal with
dropouts based on imputation techniques and zero-inflated models. For example,
CIDR (Lin, Troup and Ho, 2017), CMF-Impute (Xu et al., 2020) and ZIFA (Pierson
and Yau, 2015) impute gene expression values in dropout events based on cell
and gene similarities or parametric modeling with latent variables. However,
imputation methods usually require various conditions and may lead to new bias
and the elimination of biologically meaningful variation. SC3 (Kiselev et al.,
2017) is a popular pipeline for scRNA-seq data using combined clustering
methods, while it does not perform any form of batch effect correction.
Existing approaches to correct for batch effects include model-based
adjustment, such as Combat (Johnson, Li and Rabinovic, 2007); canonical
correlation analysis (CCA) subspace alignment, such as Seurat (Butler et al.,
2018); and mutual nearest neighbors (MNN), such as fastMNN (Haghverdi et al.,
2018) and BEER (Zhang, Wu and Tian, 2019). Unfortunately, these existing
methods have various drawbacks. For example, Combat (Johnson, Li and
Rabinovic, 2007) was initially designed for bulk cell data, and thus provides
no special treatment for dropout events. Seurat (Butler et al., 2018), fastMNN
(Haghverdi et al., 2018), BEER (Zhang, Wu and Tian, 2019), and ZINB-WaVE
(Risso et al., 2018) first perform dimensionality reduction and correct for
batch effects, followed by cell clustering. However, these multi-staged
approaches often lose information across the separate steps, as the analysis
conducted at a given step is conditional on the summarized data obtained in
the previous stage, and therefore cannot recover any useful information that
might have been lost then. In summary, none of these discussed methods provide
a unified strategy to systematically analyze scRNA-seq data.

We also remark that the most common practice in scRNA-seq data analysis is
still to analyze a single biological sample at a time, followed by manually
summarizing the common findings shared by different samples. Yet this strategy
fails to leverage information shared across samples and potentially sacrifices
the detection power. To address this drawback, our method will take into
account all the samples jointly in the new modeling system.

Specifically, we propose RZiMM-scRNA, a Regularized Zero-inflated Mixture
Model framework to simultaneously cluster cells and identify differentially
expressed genes between cell subgroups in scRNA-seq data. RZiMM-scRNA also
accounts for dropout events and sample heterogeneity. In particular,
estimation of the latent group allocation and effect is directly linked to
cell clustering membership construction. An $l_{1}$-penalty on the pairwise
differences of group effects within a gene is also incorporated to impose
sparsity on the space of differentially expressed genes as a treatment for
high dimensionality. Lastly, the parametric nature of RZiMM-scRNA naturally
allows for the development of an importance score for biomarker
identification, whose statistical significance can also be assessed. Unlike
popular two-staged methods which are unable to fully utilize information in
the data, RZiMM-scRNA is a unified model framework that addresses all of the
discussed challenges. We demonstrate the promise of the proposed method via
extensive simulation studies and applications to two brain cancer scRNA-seq
datasets.

##  2 Method

###  2.1 A New Model Framework

We start with the notations. Let $X_{ij}$ be the normalized expression of the
$j^{\text{th}}$ gene from the $i^{\text{th}}$ cell, for $i=1,...,N$,
$j=1,...,J$. Note that the cells from all the samples (patients) are collapsed
together in our analysis. To account for sparse values in the single-cell gene
expression data, we assume that $X_{ij}$ follows a zero-inflated mixture
distribution,

$\displaystyle\operatorname{P}(X_{ij}=x|\mu_{ij},\sigma_{j}^{2})=(1-\pi_{j})\operatorname{I}(x=0)+\pi_{j}\operatorname{I}(x\neq
0)\cdot\mathcal{G}(x|\mu_{ij},\sigma_{j}^{2}),$ (1)

where $\pi_{j}$ is the probability that gene $j$ has a non-zero value and
$\mathcal{G}(\cdot|\mu_{ij},\sigma_{j}^{2})$ is a specified distribution for
non-zero values of gene $j$, with mean $\mu_{ij}$ in cell $i$ and variance
$\sigma_{j}^{2}$. The proposed framework can accommodate a wide range of data
types by specifying $\mathcal{G}(\cdot|\mu_{ij},\sigma_{j}^{2})$ with
different forms of distributions, including Gaussian, Poisson, and negative
binomial, which have been widely used in modeling scRNA-seq data.

Motivated by the astrocytoma scRNA-seq data previously described and also for
the purpose of demonstration, we first focus our discussion on the Gaussian
distribution
$\mathcal{G}(x|\mu_{ij},\sigma_{j}^{2})=\exp\left\\{-\frac{1}{2}\log
2\pi\sigma_{j}^{2}-\frac{(X_{ij}-\mu_{ij})^{2}}{2\sigma_{j}^{2}}\right\\}$.
This data set exhibits a high prevalence of dropout events and 75.3% of
expression values are zeros. Figure 4a displays the histogram of the
standardized non-zero expression values from five hundred randomly selected
genes, revealing that its distribution appears symmetric, bell-shaped, and
close to a standard normal distribution. Moreover, Figure 4b depicts
histograms and density plots of the non-zero expression values for four
randomly selected genes, RGS1, A2M, C3 and FKBP5. We observe multiple peaks in
the distribution of each gene, indicating the existence of cell subgroup
structures.

As shown in a t-SNE plot (Figure 6b) of the glioma data, cells from the same
sample are mostly clustered together, illustrating the difficulty of
identifying biologically meaningful cell subgroups without adjusting for the
sample effect. To perform simultaneous cell clustering and differential gene
expression detection while accounting for the sample effect (high within-
sample correlation), we impose the following structure on the mean expression
parameter $\mu_{ij}$,

$\displaystyle\mu_{ij}=e_{j,k(i)}+\sum_{p=1}^{P}G_{ip}m_{jp},$ (2)

Herein, $e_{j,k(i)}$ corresponds to the sample effect for the $j$th gene in
sample $k(i)$, where $k(\cdot)$ is the function that gives the sample ID for
each cell. For example, $k(i)=k$ indicates that cell $i$ comes from sample
$k$; correspondingly, $e_{jk}$ represents the effect of sample $k$ for gene
$j$. $G_{ip}$ is the indicator of cell subgroup membership for cell $i$, which
is of primary interest to estimate. Specifically, $G_{ip}=1$ indicates that
cell $i$ belongs to hidden subgroup $p$, and $G_{ip^{\prime}}=0$ for all
$p^{\prime}\neq p$. In addition, $m_{jp}$ represents the corresponding effect
of gene $j$ in cell subgroup $p$. For model identification, we set $e_{jK}=0$
for gene $j=1,...,J$ in sample $K$.

Let us denote $Z_{ij}=\operatorname{I}(X_{ij}\neq 0)$. The log-likelihood can
be written as

log-L $\displaystyle=$
$\displaystyle\sum_{j=1}^{J}\sum_{i=1}^{N}(1-Z_{ij})\log(1-\pi_{j})+Z_{ij}\left[\log\pi_{j}-\frac{1}{2}\left\\{\log\left(2\pi\sigma_{j}^{2}\right)+\frac{(X_{ij}-\mu_{ij})^{2}}{\sigma_{j}^{2}}\right\\}\right].$
(3)

The other important task is to identify which genes are the main contributors
to cell subgroup differentiation. To achieve this, we further regularize the
estimation of $m_{j1},\cdots,m_{jP}$ via the sparsity-induced penalty on the
difference between $m_{jp}$ and $m_{jp^{\prime}}$, denoted by
$\Omega(\bm{m})$. That is, we look to force the trivial genes that are not
differentially expressed among the cell subgroups, which are also estimated,
to satisfy $m_{j1}=,...,=m_{jP}$. Naturally, the genes which have the largest
differences among $m_{jp}$s are likely the main drivers of cell subgroup
formation. The penalized negative log-likelihood can be written as,

$\displaystyle\ell(X,\bm{\theta})\equiv-\text{log-L}+\Omega(\bm{m})=-\text{log-L}+\frac{\lambda}{P(P-1)}\sum_{j=1}^{J}\sum_{1\leq
p<p^{\prime}\leq P}|m_{jp}-m_{jp^{\prime}}|,$ (4)

where $\bm{\theta}=\\{\bm{\pi},\bm{m},\bm{e},\bm{\sigma^{2}},\bm{G}\\}$
includes all unknown parameters to be estimated, and the tuning parameter
$\lambda$ corresponds to the weight of the penalty term $\Omega(\bm{m})$.

Incorporating the above-described parameter regularization, the new model
framework is called Regularized Zero-inflated Mixture Models for scRNA data
(RZiMM-scRNA).

###  2.2 Parameter Estimation.

We estimate $\bm{\theta}$ via minimizing $\ell(X,\bm{\theta})$, such that

$\displaystyle\widehat{\bm{\theta}}=\operatorname*{argmin}_{\bm{\theta}}\ell(X,\bm{\theta}).$
(5)

First, we obtain the estimate of $\pi_{j}$, the probability of non-zero values
for gene $j$, by solving $\partial\ell(X,\bm{\theta})/\partial\pi_{j}=0$,

$\displaystyle\widehat{\pi}_{j}=\frac{\sum_{i=1}^{n}Z_{ij}}{n}.$ (6)

Note that minimization of the loss function (4) is convex when either $\bm{G}$
or $\\{\bm{m},\bm{e},\bm{\sigma^{2}}\\}$ is fixed, but not jointly convex in
$\bm{\theta}$. Therefore, we iteratively estimate $\bm{G}$ and
$\\{\bm{m},\bm{e},\bm{\sigma^{2}}\\}$ conditional on each other, by minimizing
the following terms separated from $\bm{\pi}$,

$\displaystyle\ell(Y,\bm{\theta})=l_{0}(\bm{\pi})+\sum_{j=1}^{J}\sum_{i=1}^{N}\frac{Z_{ij}}{2}\left\\{\log(2\pi\sigma_{j}^{2})+\frac{(X_{ij}-\mu_{ij})^{2}}{\sigma_{j}^{2}}\right\\}+\Omega(\bm{m})=l_{0}(\bm{\pi})+l_{1}+\Omega(\bm{m}).$
(7)

$l_{1}$ can be written as,

$\displaystyle
l_{1}=\sum_{j=1}^{J}\frac{1}{2}\left\\{T_{j}\log(2\pi\sigma_{j}^{2})+\frac{\left(X_{j}^{\prime}-G^{(j)}\bm{m}-S^{(j)}\bm{e}\right)^{T}\left(X_{j}^{\prime}-G^{(j)}\bm{m}-S^{(j)}\bm{e}\right)}{\sigma_{j}^{2}}\right\\},$
(8)

Herein, $T_{j}=\sum_{i=1}^{n}Z_{ij}$, $X_{j}^{\prime}$ is a vector of
$X_{ij}$’s when $Z_{ij}=1$, and $G$ is a $n\times P$ matrix for the group
allocation. The $S$ is a $n\times(K-1)$ matrix, with the element $S_{ik}=1$
when cell $i$ belongs to sample $k$, and 0 otherwise. The $G^{(j)}$ and
$S^{(j)}$ include only rows corresponding to $Z_{ij}=1$.

Since the penalty term $\Omega(\bm{m})$ is non-differentiable, we adopt
Majorization Maximization (MM) algorithm (Lange, 2013) in the optimization
procedure,

$\displaystyle\Omega(\bm{m})\leq\frac{\lambda}{P(P-1)}\sum_{j=1}^{J}\sum_{1\leq
p<p^{\prime}\leq
P}\frac{1}{2}\left\\{\frac{(m_{jp}-m_{jp^{\prime}})^{2}}{b_{jpp^{\prime}}}+b_{jpp^{\prime}}\right\\}=\Omega^{\prime}(\bm{m},\bm{b})$
(9)

where $b_{jpp^{\prime}}>0$, $\Omega^{\prime}(\bm{m},\bm{b})$ is differentiable
and the equivalence holds when $b_{jpp^{\prime}}=|m_{jp}-m_{jp^{\prime}}|$.

Conditional on $G$, we can obtain the estimate of $\bm{m}$ and $\bm{e}$
iteratively as follows,

  1. 1.

Set $b_{jpp^{\prime}}=|\widehat{m}_{jp}-\widehat{m}_{jp^{\prime}}|+\delta$,
where $\delta>0$ to avoid division by zero.

  2. 2.

Estimate $\widehat{\bm{m}}$ and $\widehat{\bm{e}}$ by setting
$\frac{\partial\left\\{(I)+\Omega^{\prime}(\bm{m},\bm{b})\right\\}}{\partial(\bm{m},\bm{e})}\equiv
0$.

We follow this by solving $\frac{\partial(I)}{\partial\bm{\sigma}^{2}}\equiv
0$ to estimate $\widehat{\bm{\sigma}^{2}}$. Then, conditional on the estimate
of $\bm{m}$, $\bm{e}$ and $\bm{\sigma}^{2}$, for $i=1,...,n$, we estimate
$G_{i}$ in the following,

$\displaystyle\widehat{p}(i)=\operatorname*{argmin}_{p}\sum_{j=1}^{J}\frac{\left(X_{ij}-\widehat{m}_{jp}-\widehat{e}_{j,k(i)}\right)^{2}}{2\widehat{\sigma}_{j}^{2}},$
(10)

where $\widehat{G}_{i}=(\widehat{G}_{i1},...,\widehat{G}_{iP})$,
$\widehat{G}_{ip}=1$ when $p=\widehat{p}(i)$, and 0 otherwise.

The complete procedure can be summarized as follows,

  1. 1.

Estimate $\widehat{\bm{\pi}}$ via (6). Initialize $\widehat{G}$ by randomly
assigning cell $i$ to group $p$, $1\leq p\leq P$.

  2. 2.

Iterate the following procedure until convergence or the maximum number of
iterations is reached.

    1. (a)

Fix $\widehat{G}$, and estimate $\widehat{\bm{m}}$, $\widehat{\bm{e}}$ and
$\widehat{\bm{\sigma}}^{2}$ via MM algorithm.

    2. (b)

Fix $\widehat{\bm{m}}$, and estimate $\widehat{\bm{e}}$ and
$\widehat{\bm{\sigma}}^{2}$, and $\widehat{G}$ via (10).

The tuning parameter $\lambda$ of the penalty term $\Omega(\bm{m})$ and the
optimal number of cell types can be automatically determined based on the
modified BIC proposed by Wang, Li and Leng,

$\operatorname{BIC}_{\gamma}=-2\operatorname{log-L}+\gamma d\log N$

where $d$ is the number of total parameters, $\gamma>1$ is a pre-specified
parameter. We set $\gamma=\log(\log J)$ as suggested by Wang, Li and Leng. As
shown by Wang, Li and Leng, this modification improves the performance over
the original BIC, which tends to obtain over-parsimonious results on
clustering.

###  2.3 Gene Importance Scoring.

The new model framework not only achieves cell subgroup (type) discovery in a
data-adaptive fashion, but also provides meaningful insight in detection of
important genes. To identify the differentially expressed genes which drive
formation of different cell types by RZiMM-scRNA, we define an importance
score for each gene as the summation of pairwise differences between cell
types,

$\displaystyle\operatorname{IS}_{j}=\frac{1}{P(P-1)}\sum_{1\leq
p<p^{\prime}\leq P}|m_{jp}-m_{jp^{\prime}}|$ (11)

Genes with large values of $\text{IS}_{j}$ are considered to be cluster-
distinctive biomarker candidates.

###  2.4 Cell-Type-Specific Driving Differential Gene Expression.

It is also straightforward to utilize RZiMM-scRNA to identify genes that are
non-trivially over-expressed or under-expressed in a specific cell type. We
define the measure of differential expression level for gene $j$ in cell type
$p$ as,

$\displaystyle\operatorname{DE}_{jp}=\frac{1}{P-1}\sum_{p^{\prime}\neq
p}\left(m_{jp}-m_{jp^{\prime}}\right)$ (12)

Genes with the largest magnitudes of $\text{DE}_{jp}$ have the most over-
expression or under-expression in cell type $p$, and are viewed as the main
contributors to separate cell type $p$ from the others. This measure is useful
in aligning cell clusters to known cell types and identifying new biomarkers
of a specific cell type.

###  2.5 A remark

The focus of this paper is to introduce and apply a new joint-sample modeling
framework for simultaneous discovery of cell subtypes and detection of driving
gene markers. This new system yields an intuitive measure of gene importance,
which can be used to directly rank genes according to their strength of
association with the discovered cell subgroups. Statistical significance can
also be assessed through bootstrapping, which will be deferred to our future
investigation.

##  3 Results

We conduct extensive simulation studies to assess the performance of RZiMM-
scRNA in some common settings and conclude with an application to the
previously described astrocytoma scRNA-seq datasets.

In addition to RZiMM-scRNA, we also implement Hierarchical Clustering
(H-Clust), K-means, Seurat, SC3, fastMNN, and RZiMM not adjusting for within-
sample effect (RZiMM-Naive) for comparison. H-Clust and K-Means are used to
cluster the cells and ANOVA are then conducted using the “anova” function in R
package “stats” (R Core Team, 2018) to detect significant genes based on
likelihood-ratio test. FastMNN extracts fifty top features in the PC space
based on PCA with batch effects corrected and then performs cell clustering
via K-means. The gene statistical significance level is then determined in the
same way as in H-Clust and K-Means. Seurat and SC3 have their own pipelines to
identify differentially expressed genes.

K-means and H-Clust are implemented using the functions “kmeans” and “hclust”,
respectively, in R package “stats”. For H-Clust, we use one minus Pearson
correlation as the dissimilarity measurement coupled with average
agglomeration. We then use “cutree” function in R package “stats” to divide
the ordered cells into clusters. For K-Means, Euclidean distance is used as
the dissimilarity measurement. FastMNN is implemented using function “fastMNN”
in R package “batchelor”. Seurat and SC3 are implemented using R packages
“Seurat” and “SC3”.

###  3.1 Simulation Study

####  3.1.1 Expression Data from Zero-inflated Gaussian Distribution

In this study, we generate zero-inflated expression data from models (1) and
(2), assuming $\mathcal{G}$ is a Gaussian distribution, which is motivated by
the first brain cancer study. We consider $N=1000$ cells, $J=1000$ genes,
$K=10$ samples, and $P=4$ hidden cell subgroups, where sample ID $k(i)$ and
subgroup membership $p(i)$ are randomly assigned. Each of the genes $1\leq
j\leq 100$ is a differentially expressed (DE) gene, where a randomly selected
cell subgroup is designed to have a higher mean expression than the other
three subgroups by a value drawn from the uniform distribution
$\operatorname{U}(\delta,2\delta)$, while the other three groups have a
constant effect $m_{j}$ drawn from $\operatorname{U}(4,8)$. The effect size,
$\delta$ ranging from 0.5 to 1.5, indicates the difference in DE genes between
subgroups. For genes $100<j\leq J$, we set $m_{j1}=...=m_{jP}=m_{j}$ where
$m_{j}\sim\operatorname{U}(4,8)$. That is, these genes are not differentially
expressed among the cell subgroups. In each scenario of $\delta$ value, one
hundred simulated data sets are used for investigation. The gene-specific
sample effect $e_{jk}$ is drawn from the normal distribution
$\operatorname{N}(0,0.2)$, for $k=1,...,K$, $j=1,...,J$. Therefore, the mean
expression $\mu_{ij}=m_{j,p(i)}+e_{j,k(i)}$. Further, the probability of zero-
expression $\pi_{ij}=\operatorname{expit}(\mu_{ij}/2-4)$, and the gene
expression intensity values are generated as $X_{ij}=Z_{ij}W_{ij}$, where
$Z_{ij}\sim\operatorname{Bin}(1,\pi_{ij})$ and
$W_{ij}\sim\operatorname{N}(\mu_{ij},1)$. The within-sample correlation
structure at $\delta=0.5$ is visualized in Supplement Figure S1, where the
positive correlation among cells within a sample (block-wise) is clearly
observed.

Figure 1: Performance Plots for the Simulation Study with Normalized Data. a.
ARI for the clustering result. b. AUC for important gene identification.

We assume four cell clusters (also the truth) for all methods, as our focus is
to compare cell subgroup composition and biomarker detection among the various
methods. To evaluate the clustering performance, we use the Adjusted Rand
Index (ARI) (Hubert and Arabie, 1985). In addition, we use the area under ROC
curve (AUC) as the measure of biomarker detection accuracy, where larger
values of AUC correspond to higher biomarker detection accuracy. Figure 1a-b
plots the clustering and biomarker detection accuracies for each of the
methods averaged over the one hundred simulated data sets against various
values of $\delta$. It is clear that RZiMM-scRNA consistently outperforms all
the other methods in both clustering and biomarker detection accuracy. As the
effect size $\delta$ increases, the ARI for clustering and AUC for DE gene
identification from all methods are improved. The value of $\delta$ needed for
RZiMM-scRNA to reach perfect ARI and AUC is around 1, and that for other
methods is greater than 1.5. This promising performance on zero-inflated
expression data drawn from the Gaussian distribution suggests that RZiMM-scRNA
is a powerful method to analyze the astrocytoma scRNA-seq expression data.

####  3.1.2 Count Data from Zero-inflated Poisson Distribution

In this study, the popular Poisson distribution is considered to mimic the
count nature of the data generated by single-cell RNA-sequencing. Similar to
the previous simulation setting, we consider $N=1000$ cells, $J=1000$ genes,
$K=10$ samples, and $P=4$ hidden cell subgroups. Each of the genes $1\leq
j\leq 100$ is a differentially expressed (DE) gene, where a randomly selected
cell subgroup $p^{(j)}$ is designed to have a higher mean expression than the
other three subgroups by a value drawn from the uniform distribution
$\operatorname{U}(\delta,2\delta)$, while the other three groups have a
constant effect parameter $m_{j}$ drawn from $\operatorname{U}(4,8)$. The
effect size, $\delta$ ranging from 1 to 5, indicates the degree of difference
in DE genes between subgroups. For genes $100<j\leq J$, we set
$m_{j1}=...=m_{jP}=m_{j}$ where $m_{j}\sim\operatorname{U}(4,8)$. The gene-
specific sample effect $e_{jk}$ is allowed to vary, drawn from the normal
distribution $\operatorname{N}(0,2)$, for $k=1,...,K$, $j=1,...,J$. The mean
expression $\mu_{ij}=\max(m_{j,p(i)}+e_{j,k(i)},0.5)$. Further, the
probability of zero-expression $\pi_{ij}=\operatorname{expit}(\mu_{ij}/6-2)$,
and thus $X_{ij}=Z_{ij}W_{ij}$, where
$Z_{ij}\sim\operatorname{Bin}(1,\pi_{ij})$ and
$W_{ij}\sim\operatorname{Poisson}(\mu_{ij})$, resulting in the count data
following zero-inflated Poisson distributions.

As discussed in Section 3, the proposed model framework can flexibly
accommodate non-Gaussian distributed data. We use RZiMM-Pois as the
abbreviation of out method implemented under zero-inflated Poisson
distribution. We are also interested to examine the performance of adopting
zero-inflated negative binomial distribution (RZiMM-NB) and zero-inflated
Gaussian distribution (RZiMM-scRNA). From the comparison, we hope to assess
the robustness of results to various distributions. Note that the result under
the true Poisson distribution (RZIMM-Pois) is treated as the benchmark.

Figure 2 displays the ARI and AUC results of all the considered methods. Note
that RZIMM-NB’s performance is very close to the benchmark RZIMM-Pois, which
is not surprising due to their intrinsic mathematical connection.
Interestingly, the next best performers are RZIMM-scRNA and fastMNN. The
fairly good performance of RZIMM-scRNA suggests that the Gaussian distribution
works robustly and reasonably well, regardless of the true underlying data
distributions. In addition, both popular pipeline methods Seurat and SC3
demonstrate inferior performance. Overall, the result demonstrates that the
RZiMM model framework is a powerful substitute to the most common alternatives
of scRNA-seq data analysis.

Figure 2: Performance Plots for the Simulation Study with Count Data. a. ARI
for the clustering result. b. AUC for important gene identification.

###  3.2 A Semi-simulation Scenario based on the Real Astrocytoma Data

We would like to make the simulated data as close to the reality as possible,
while being able to generate some truth structure for performance assessment.
Instead of generating the data from parametric models, we fully utilize the
first astrocytoma scRNA-seq dataset via permutation. Specifically, we first
randomly permute expression intensities of each gene across cells within a
sample, as to remove the unknown underlying cell subgroup structure while
preserving the within-sample correlation. Next, we add the true structure of
cell subgroups and the corresponding driving gene biomarkers as follows. We
randomly divide the cells into four clusters and select one hundred genes as
the pseudo-important biomarkers. For each of these important genes, we elevate
the expression intensities in a randomly assigned cell subgroup by a positive
value drawn from the uniform distribution $\operatorname{U}(v,v+1)$, where
several values of $v$ between $0$ and $4$ are examined.

Note that $v$ corresponds to the effect size of the hidden cell subgroups,
where larger values of $v$ correspond to more distinct cell subgroups and less
difficulty for a method to correctly identify the underlying cell subgroup
structure. The cell correlation map at $v=4$ is presented in Supplement Figure
S2, in which cells are ordered by sample and then by true hidden subgroup
within each sample. The within-sample correlation is quite noticeable in the
glioma data, as shown in different blocks. For example, the second sample to
the left MGH61 (shown in gray) reveals a strong between-cell correlation. The
small blocks in the map manifest the non-trivial correlation within a cell
subgroup.

We again compare the clustering and biomarker detection accuracies across
methods (see Figure 1c-d). We observe that RZiMM-scRNA is the most powerful of
all the methods, as it starts accurately detecting the underlying true cell
subgroup structure with the smallest effect size $v$. Specifically, we observe
that RZiMM-scRNA can detect true cell clusters with 100% accuracy for effect
sizes of at least $1.1$. In contrast, Seurat can only achieve perfect accuracy
when the effect size is greater than $4$. Since SC3, RZiMM-Naive, K-Means, and
H-Clust are unable to account for the sample effect, they completely fail to
detect the true cell subgroups when $v\leq 3$ and cannot identify the
important genes (low AUC in the right panel). Only as the effect size becomes
very large (specifically, $v=4$) can RZiMM-Naive fully detect the true
subgroup structure and important genes, whereas SC3, K-Means and H-Clust
cannot improve their performance.

This study also suggests that RZiMM-scRNA, compared to the other methods, will
likely provide more meaningful insight from the astrocytoma data, as the
empirical distributions of gene expression data in that real study are
essentially preserved in our simulation.

Figure 3: Performance Plots for the Semi-simulation Study. a. ARI for
clustering. b. AUC for important gene identification.

###  3.3 Analysis of Astrocytoma Data

Finally, we apply each of the methods to the first astrocytoma scRNA-seq
dataset. In the study done by Venteicher et al., ten IDH-mutant human
astrocytoma tumors spanning grades II to IV were disaggregated, sorted into
single cells, and profiled by Smart-seq2 (Picelli et al., 2014). Expression
levels were quantified as
$\operatorname{E}_{jk}=\log_{2}(\text{TPM}_{jk}/10+1)$, where
$\text{TPM}_{jk}$ refers to transcript-per-million for gene $j$ in sample $k$,
as calculated using the software RSEM (RNA-Seq by Expectation-Maximization, Li
and Dewey (2011)). Cells with fewer than 3,000 detected genes or an average
expression level of a curated list of housekeeping genes below 2.5 were
excluded. Consequently, 6,341 single cells were retained in the study. In
accordance with the analysis conducted by Venteicher et al., we include 8,231
genes among the total of 23,686 for further analysis, whose expression levels
satisfy $\log_{2}[\text{Average}(\text{TPM}_{j,1...N})+1]\geq 4$.

Figure 4: Distribution of Expression in Astrocytoma scRNA-seq Data. a.
Histogram of standardized non-zero expression values from five hundred
randomly selected genes. b. Histograms (in bar) and density (in line) plots of
non-zero expression values from four genes, RGS1, A2M, C3 and FKP5.

An important first task we must address is determining the optimal number of
cell clusters for each method. For RZiMM-scRNA, we adopt modified BIC (Wang,
Li and Leng, 2009) as the model selection criterion integrated with the model-
based method, with a grid search over $P\in\\{2,3,...,15\\}$ and
$\lambda\in\\{0,20,40,...,200\\}$ (see Supplement Figure S4). As the result,
we obtain $P=6$ as the optimal number of clusters. For K-Means, H-Clust, and
fastMNN, we utilize average silhouette value (Rousseeuw, 1987) and gap
statistics (Tibshirani, Walther and Hastie, 2001), examining the range of $P$
from $2$ to $15$. The average silhouette value method is implemented using the
“fviz_nbclust” function in R package “factoextra” (Kassambara and Mundt, 2019)
and the gap statistic is implemented using MATLAB package
“clustering.evaluation”. The average silhouette value method suggests $2$ cell
clusters for the three methods, which is too small to capture the important
pattern in our data. On the other hand, the gap statistics suggests $6$, $15$,
and $14$ clusters for K-Means, H-Clust and fastMNN, respectively (see
Supplement Figure S5). Taking into account the discordance between the average
silhouette value and gap statistic suggestions, and for fair comparison, we
use $6$ clusters for these methods.

Figure 6a is a compositional plot displaying the proportions of cells
clustered into each cell group for each sample. We notice that for several
patient samples, such as MGH56 and MGH61, SC3, K-Means and H-Clust tend to
cluster cells together by sample. This is also observed for RZiMM-Naive and
SC3. For example, the cells of sample MGH56 are mostly clustered into Cluster
5 (purple) by RZiMM-Naive and SC3, or Cluster 3 (yellow) by both K-Means and
H-Clust. Once again, we see the importance of addressing sample effect.
Without the ability to handle the sample effect, the clustering patterns of
RZiMM-Naive, SC3 K-Means, and H-Clust are dominated by sample identity, rather
than biologically important distinctions. These same clustering patterns are
observed in a correlation map (Figure 5), which reveals prominent within-
sample correlations, and in t-SNE plots (Figure 6b), where cells appear to be
clustered by sample identity.

Figure 5: Correlation Map for Astrocytoma scRNA-seq Data. The genes are
ordered by RZiMM-scRNA and H-Clust.  Figure 6: Clustering Analyses of the
Astrocytoma Data. a. Compositional plots for cluster membership of each sample
estimated by the five methods. b. t-SNE plots colored by sample or clustering
methods.  Figure 7: Pairwise ARI for Sample Index and Clustering Results

The cell correlation map (Figure 5) presents the pairwise Pearson correlation
between cells and the heatmap (Supp Figure S3) shows the expression intensity
of each gene in each cell, both of which are annotated with sample labels and
the clustering membership estimated by each method. The cells are first
ordered by the six clusters identified by RZiMM-scRNA, and then ordered within
each cluster by H-Clust. Note that Cluster 1 from RZiMM-scRNA (red) is
identified by six clustering methods, as it contains highly correlated cells,
but not by Seurat as shown by the correlation map in Figure 5. The heatmap in
Supp Figure S3 also reveals that a large group of genes are highly expressed
in Cluster 1, which is made easily identifiable by all the clustering methods.
Interestingly, Cluster 1 matches the microglial group of cells reported in
Venteicher et al.. Microglial cells mediate immune responses in the central
nervous system and are re-programmed by tumor cells to create a favorable
microenvironment for cancer progression (Abels et al., 2019). Moreover,
Cluster 2 from RZiMM-scRNA (blue) and Seurat is composed of a small group of
highly correlated cells, as shown in the correlation map. In fact, Cluster 2
matches another biologically meaningful cell type in the oligodendrocyte group
reported by Venteicher et al.. Oligodendrocytes produce myelin in the central
nervous system and its precursor cells help modulate particular
microenvironments for glioma cells (Hide et al., 2018). However, none of the
aforementioned methods other than RZiMM-scRNA or Seurat are able to
distinguish this cell type from the other cells in our data. Specifically,
these other methods tend to combine cells in Cluster 2 with some cells in
Clusters 5 and 6, as identified by RZiMM-scRNA, into the same subgroup,
indicated by blue for all the other methods.

We also visualize this result using two-dimensional t-SNE plots in Figure 6b.
In the top left panel, the plot is colored by sample; while in other panels,
the plot is colored by the clustering result for each method. Again, Cluster 2
is revealed to be clearly separated from the other cells. However, only RZiMM-
scRNA and Seurat can sharply identify this pattern. All the other methods pool
many other cells, which appear to be distant, into this cluster (blue). This
demonstrates the superior performance of RZiMM-scRNA, as well as Seurat in
identifying rare cell types.

In Figure 7, pairwise ARIs across different clustering methods and sample ID
are presented. It is noteworthy that when comparing RZiMM-scRNA, fastMNN and
Seurat results to sample ID, the ARI is close to zero, while other methods
have fairly high ARIs to sample ID. This result indicates that RZiMM-scRNA,
fastMNN and Seurat can effectively remove the batch effects, while the other
methods fail. In addition, RZiMM-scRNA shows better agreement with RZiMM-Naive
and fastMNN, indicated by fairly large ARI values. In contrast, Seurat has
overall low agreement with all the other methods.

For the particularly interesting Cluster 2 corresponding to the
oligodendrocyte cell type, we also compute the measure of differential gene
expression between Cluster 2 and all the other clusters according to equation
(13) in RZiMM-scRNA, which is displayed in the top panel of Figure 8. The
under-expressed genes in Cluster 2 are colored in blue, the top five over-
expressed genes in Cluster 2 ($\operatorname{DE}_{j2}>6$) are colored in
orange, and the rest are colored in gray. We observe a set of genes highly
over-expressed in Cluster 2 with the differential expression measure over
four. In particular, the top five over-expressed genes in Cluster 2, PLP1,
CLDN11, MOG, MBP and TF, are colored in orange. These five genes are also
among the top forty important genes identified by RZiMM-scRNA, which yield the
overall largest differentiation across the six cell clusters, presented in the
bottom panel of Figure 8 and Supplement Table S1. In Supplement Table S1, we
observe that these forty genes are either over-expressed in Cluster 1 or 2 or
under-expressed in Cluster 1 or 2, as compared to all the other clusters.

Furthermore, we investigate the biological significance of the identified top
forty genes. Promisingly, proteins encoded by most of these genes play roles
in the nervous system ranging from immune response to cell migration and
proliferation, and many are associated with forms of brain cancer, brain
disease, and other ailments. We include a brief description of the biological
significance of some of these genes in Table LABEL:tab:imp, specifically,
their relationship to cancers such as low-grade glioma (LGG), high-grade
glioma (HGG), and glioblastoma (GBM), and neurological and neurodegenerative
diseases. For a more extensive account of the function and biological
significance of these genes, please refer to Table S2 in the Supplement.

Figure 8: RZiMM-scRNA Importance. Top panel: mean difference in normalized
expression between Cluster 2 and the rest (orange: top five over-expressed
genes in Cluster 2, PLP1, CLDN11, MOG, MBP and TF; grey: over-expressed genes
except the top five; blue: the under-expressed genes in Cluster 2). Bottom
panel: importance scores for differentially expressed genes (orange: top forty
differentially expressed genes; grey: other genes)

.

Table 1: Top Important Genes Gene | Significance | Ref.
---|---|---
SPP1 |  \- Upregulated by glioma-associated microglia/macrophages (GAMs) and associated with poor patient prognosis in LGG and GBM |  Chen et al. (2019a) Szulzewsky et al. (2015)
PLP1 |  \- Key for glioma development and associated with patient prognosis |  Chen et al. (2019b)
BCAN |  \- May promote the growth and motility of brain tumor cells |  O’Leary et al. (2016)
RGS1 |  \- Gene polymorphisms associated with MS |  Caballero-Franco and Kissler (2016) Sawcer et al. (2011)
CD74 |  \- Upregulated in HGG and positively associated with patient survival |  Zeiner et al. (2015)
GPM6A |  \- Identified as a top gene candidate as a causal gene for SZ |  Ma et al. (2018)
PTPRZ1 |  \- Upregulated in LGG and GBM, and associated with malignant growth of GBM and poor patient prognosis |  Xia et al. (2019) Shi et al. (2017)
CCL3 |  \- Downregulated in HGG |  Vidyarthi et al. (2019)
SEPP1 |  \- Inactivation causes severe neurological dysfunction and brainstem neurodegeneration in mice |  Byrns et al. (2014)
GSN |  \- Significantly downregulated in many cancers including GBM |  Miyauchi et al. (2018)
FXYD6 |  \- Downregulated in brain of AD mouse model |  George et al. (2010)
MOG |  \- Anti-MOG antibodies associated with demyelinating diseases such as acute disseminated encephalomyelitis and neuromyelitis optica spectrum disorder |  Sinmaz et al. (2016)
MBP |  \- Increased expression levels detected in AD cortex and could be involved in amyloid plaque formation |  Zhan et al. (2015)
CCL4 |  \- Upregulated by GAMs |  Zeiner et al. (2019)
A2M |  \- Inhibits proliferation, migration, and invasion of astrocytoma cells |  Lindner et al. (2010)
APOE |  \- Upregulated in grade II-IV astrocytoma |  Rousseau et al. (2006)
CRYAB |  \- Upregulated in IDH1${}^{\text{R132H}}$ mutant anaplastic astrocytoma and GBM |  Avliyakulov et al. (2014) Kore and Abraham (2014)
CLDN11 |  \- Absence in CNS myelin perturbs behavior in mice and could lead to neuropsychiatric disease |  Maheras et al. (2018)
FEZ1 |  \- Abnormal, unphosphorylated co-aggregates present in brains of transgenic AD mice models |  Butkevich et al. (2016)
HLA-E |  \- Upregulated in grade II-IV astrocytoma |  Mittelbronn et al. (2007)
TF |  \- May play role in increased iron deposition in MS brain |  Khalil et al. (2014)
HLA-DRA |  \- Low expression levels associated with poor prognosis in pediatric adrenocortical tumors |  Leite et al. (2014)
C1orf61 |  \- May participate in occurrence and development of GBM |  Jia et al. (2019)
APLP1 |  \- May contribute to pathogenesis of AD-associated neurodegeneration |  Bayer et al. (1997)
CD24 |  \- Over-expression stimulates glioma cell migration, invasion, colony formation, and tumor growth in mice |  Barash et al. (2019)
LAPTM5 |  \- May be used by glioma to negatively modulate antitumoral response of immune system |  M. Hajj et al. (2020)
C1QB |  \- Upregulated in diverse grade of gliomas and plays important role in pathogenesis |  Mangogna et al. (2019)
IL1B |  \- Highly expressed in GBM and promotes glioma spheroid formation |  Wang et al. (2019)

We see in Table LABEL:tab:imp that PLP1, BCAN, and PTPRZ1 are all associated
with the development and malignant growth of glioma. This is a striking
observation since, in terms of gene importance score, three of the seven top
differentially expressed genes play notable roles in glioma progression and
malignancy. Along with this consideration, we observe that a large proportion
of the top forty important genes are associated with brain cancer and
neurological or neurodegenerative disease. This shows that RZiMM-scRNA is
quite effective in detecting biologically important genes in real scRNA-seq
data.

It is also noteworthy that many of the top over-expressed genes in each
cluster established by our new method serve similar biological functions. Upon
evaluation of cluster-specific expression difference (12), we determine the
top over- and under-expressed genes in specific cell clusters. For example, we
find that the top upregulated genes in Cluster 1 are RGS1, CD74, CCL3, CCL4,
and HLA-DRA. We see in Table S2 that these genes are involved in immune
function and regulation, revealing that particular cells dealing with immune
system mediation in glioma are grouped together by RZiMM-scRNA. Likewise, we
observe that PLP1, CLDN11, MOG, and MBP are all highly over-expressed in
Cluster 2 and note that these genes play roles in oligodendrocyte development
and myelination. The fact that the top upregulated genes in each cluster serve
similar biological functions suggests that RZiMM-scRNA is effective in
clustering individual cells into groups in accordance with their biological
functions.

In addition to new biomarkers discovered by our method, our findings are also
supported by the previous analysis done by Venteicher et al.. They reported
the presence of an inflammatory program active in GAMs consisting of cytokines
and chemokines, in which genes such as CCL3, CCL4, and IL1B are active. This
is consistent with our result that these genes, along with other genes
involved in immune regulation, are over-expressed in Cluster 1. Likewise, the
previous analysis discovered a group of nonmalignant oligodendrocytes in their
clustering that preferentially expressed markers such as MBP and CLDN11. These
genes are also upregulated in our Cluster 2 along with other oligodendrocytic
markers. Moreover, the efficacy of our method is proven by our finding that
the most under-expressed genes across these two cell groups include BCAN,
PTPRZ1, and GPM6A. Since these genes are involved in neuronal migration and
the growth of brain tumor cells, the downregulation of these genes is
consistent with these cells being nonmalignant, which is confirmed by the
properties we discuss.

###  3.4 Analysis of High-Grade Astrocytoma Data

In this study of aggressive brain tumor types, conducted by Yuan et al.,
samples from eight patients were taken of the cellular milieu of high-grade
gliomas, namely glioblastoma and anaplastic astrocytoma. Taken together,
$N=23,793$ total cells were profiled across $J=60,725$ total genes, among
which the majority have very sparse data. Posterior to the similar filtering
procedure conducted in the first study, we retain 3550 genes in our data
analysis.

Based on the empirical distribution of this dataset, a zero-inflated Poisson
distribution seems to provide better fit than the Gaussian distribution. We
thus implement our algorithm assuming the zero-inflated Poisson distribution,
in which the number of cell subgroups is fixed to six, as suggested in the
original study by Yuan et al. for fair comparisons.

Out of the top $40$ genes identified by our method, $16$ of these are
important genes found in the previous analysis—KLK4, KIF18B, NCAPH, NCAPG,
SUSD3, TROAP, ABCB1, MYBPC1, KIFC1, GTSE1, PRC1, MAD2L1, TPX2, and EFHD1. This
finding helps support the meaningfulness of our obtained result. In addition,
our method demonstrates its potential to make new discoveries as we continue
to investigate the biological significance of the top genes which are not
discussed in the previous work. Excitingly, we find that many of these genes
are associated with cancers as well as neurological and neurodegenerative
diseases. We include a brief description of the biological significance of
some of these genes in Table LABEL:tab:imp2.

Table 2: Top Genes Not Previously Discussed Gene | Significance | Ref.
---|---|---
CAMK2B |  \- Downregulated in GBM and glioma, implicated in neurodevelopment, brain function, learning and memory processes |  Xiong et al. (2019) Johansson, Göransson and Westermark (2005) Nicole and Pacary (2020)
LZTS2 |  \- Tumor suppressor gene influencing proliferation of cancer cell lines such as breast, colon, prostate, and glioma |  Thyssen et al. (2006) Kim et al. (2011)
CBR1 |  \- Decreased expression associated with tumor metastasis and poor prognosis \- Upregulated in AD brains \- Partially responsible for down syndrome cognitive impairments |  Murakami et al. (2011) Balcz et al. (2001) Arima-Yoshida et al. (2020)
PAQR8 |  \- Downregulated in GBM and malignant endometrial cancer tissue |  Moral-Morales et al. (2020) Sinreih et al. (2018)
KANK2 |  \- Potential tumor suppressor gene involved in regulation of apoptosis, cell cycle, and transcription |  Zhang et al. (2018) Wang et al. (2014)
ZNF713 |  \- Misregulation in brain may be involved in pathogenicity of autistic disorder |  Metsu et al. (2014)
USP54 |  \- Overexpressed in colorectal cancer stem cells and promotes intestinal tumorigenesis |  Fraile et al. (2016)
MSANTD3 |  \- May be associated with familial squamous cell lung carcinoma |  Byun et al. (2018)
DYNLL1P1 |  \- May be associated with schizophrenia |  Peyrot and Price (2021)

We observe in Table LABEL:tab:imp2 that many of our identified new biomarkers
are connected to tumorigenicity, neurodevelopment, and neurological diseases.
For example, CAMK2B is found to be downregulated in GBM and is implicated in
neuro-development, brain function, learning and memory processes. Likwise, we
find that PAQR8 is downregulated in GBM and malignant endometrial cancer
tissues. The application of our method to this high-grade astrocytoma data
further demonstrates the effectiveness of RZiMM-scRNA in detecting
biologically meaningful genes in real scRNA-seq data.

##  4 Discussion

Single-cell RNA-seq is a valuable technique for gaining important biological
insight into the formation and progression of heterogeneous diseases. However,
using inappropriate methods to analyze such data can result in incomplete or
misleading biological conclusions. Compared to bulk sample RNA-seq, the data
generated from scRNA-seq possess unique characteristics, including the
presence of dropouts and batch/sample effects, which brings additional
challenges into data analysis. Without proper treatment for these issues,
biologically meaningful discoveries, such as underlying cell heterogeneity
patterns, cannot be made by unsupervised clustering algorithms. In this paper,
we present a novel analytical method — RZiMM-scRNA, a unified statistical
framework to simultaneously discover cell subtypes and identify differential
gene expression, while accounting for both dropouts and sample effects. We
illustrate its superior performance in terms of both clustering and biomarker
identification over several widely used methods including Seurat, SC3,
K-Means, H-Clust, and fastMNN, via comprehensive simulation studies in some
commmon settings and an application to astrocytoma scRNA-seq data.

In the simulation studies, RZiMM-scRNA demonstrates its effectiveness in
discovering hidden cell subgroups and detecting important biomarkers, despite
the presence of dropout events and sample effects. As demonstrated by the
count data simulation study, our new framework can be straightforwardly
implemented under non-Gaussian distributions, for example, the zero-inflated
Poisson or negative binomial models that are widely used for the count data.
Interestingly, we also observe that the performance of adopting zero-inflated
Gaussian distribution in the proposed framework is quite robust to different
underlying distribution forms that are widely considered for scRNA-Seq data.

Moreover, we find during the first real astrocytoma scRNA-seq data application
that RZiMM-scRNA identifies two biologically meaningful cell types, previously
annotated as microglia and oligodendrocytes by Venteicher et al., while the
other methods fail to identify the small oligodendrocyte type. Utilizing the
model framework of RZiMM-scRNA, we are able to identify top genes that are
significantly over-expressed in the oligodendrocyte type, some of which are
also reported by Venteicher et al.. The informative paradigm of RZiMM-scRNA
also allows us to identify top genes which are differentially expressed across
all the hidden cell subgroups, most of which are confirmed to be GBM-related
by the related literature. Likewise, the second real data application to high-
grade astrocytoma tumors further demonstrates the promise of our model
framework in terms of biomarker detection.

In our future investigation, we will develop bootstrapping methods to assess
the statistical significance of identified genes that differentiate between
discovered cell subtypes.

The work was partially supported by the National Institute of Health Grants
NCI 5P30 CA013696, NCI P01 CA098101, NIAID 1R01 AI143886, and NCI 1R01
CA219896.

Supplementary Information We provide additional plots and tables for the
results section in Supplementary Information. Code and Software We have
developed an R package “RZiMM”, to implement our proposed method, and scripts
to perform simulational and real applications.

## References

  * Abels et al. (2019) [author] Abels, Erik RE. R., Maas, Sybren LNS. L., Nieland, LisaL., Wei, ZhiyunZ., Cheah, Pike SeeP. S., Tai, EricE., Kolsteeg, Christy-JoyC.-J., Dusoswa, Sophie AS. A., Ting, David TD. T., Hickman, SuzanneS. et al. (2019). Glioblastoma-associated microglia reprogramming is mediated by functional transfer of extracellular miR-21. Cell reports 28 3105–3119.
  * Arima-Yoshida et al. (2020) [author] Arima-Yoshida, FumikoF., Raveau, MatthieuM., Shimohata, AtsushiA., Amano, KenjiK., Fukushima, AkihiroA., Watanave, MasashiM., Kobayashi, ShizukaS., Hattori, SatokoS., Usui, MasayaM., Sago, HaruhikoH. et al. (2020). Impairment of spatial memory accuracy improved by Cbr1 copy number resumption and GABA B receptor-dependent enhancement of synaptic inhibition in Down syndrome model mice. Scientific reports 10 1–15.
  * Avliyakulov et al. (2014) [author] Avliyakulov, Nuraly KN. K., Rajavel, Kavitha SK. S., Le, Khanh Minh TK. M. T., Guo, LeaL., Mirsadraei, LeiliL., Yong, William HW. H., Liau, Linda ML. M., Li, SichenS., Lai, AlbertA., Nghiemphu, Phioanh LP. L. et al. (2014). C-terminally truncated form of $\alpha$B-crystallin is associated with IDH1 R132H mutation in anaplastic astrocytoma. Journal of neuro-oncology 117 53–65.
  * Balcz et al. (2001) [author] Balcz, BB., Kirchner, LL., Cairns, NN., Fountoulakis, MM. and Lubec, GG. (2001). Increased brain protein levels of carbonyl reductase and alcohol dehydrogenase in Down syndrome and Alzheimer’s disease. In Protein Expression in Down Syndrome Brain 193–201. Springer.
  * Barash et al. (2019) [author] Barash, UriU., Spyrou, ArgyrisA., Liu, PeiP., Vlodavsky, EuvgeniE., Zhu, ChenchenC., Luo, JuanjuanJ., Su, DongshengD., Ilan, NetaN., Forsberg-Nilsson, KarinK., Vlodavsky, IsraelI. et al. (2019). Heparanase promotes glioma progression via enhancing CD24 expression. International journal of cancer 145 1596–1608.
  * Barrett et al. (2012) [author] Barrett, TanyaT., Wilhite, Stephen ES. E., Ledoux, PierreP., Evangelista, CarlosC., Kim, Irene FI. F., Tomashevsky, MaximM., Marshall, Kimberly AK. A., Phillippy, Katherine HK. H., Sherman, Patti MP. M., Holko, MichelleM. et al. (2012). NCBI GEO: archive for functional genomics data sets—update. Nucleic acids research 41 D991–D995.
  * Bayer et al. (1997) [author] Bayer, Thomas AT. A., Paliga, KrzysztofK., Weggen, SaschaS., Wiestler, Otmar DO. D., Beyreuther, KonradK. and Multhaup, GerdG. (1997). Amyloid precursor-like protein 1 accumulates in neuritic plaques in Alzheimer’s disease. Acta neuropathologica 94 519–524.
  * Butkevich et al. (2016) [author] Butkevich, EugeniaE., Härtig, WolfgangW., Nikolov, MiroslavM., Erck, ChristianC., Grosche, JensJ., Urlaub, HenningH., Schmidt, Christoph FC. F., Klopfenstein, Dieter RD. R. and Chua, John Jia EnJ. J. E. (2016). Phosphorylation of FEZ1 by Microtubule Affinity Regulating Kinases regulates its function in presynaptic protein trafficking. Scientific reports 6 1–15.
  * Butler et al. (2018) [author] Butler, AndrewA., Hoffman, PaulP., Smibert, PeterP., Papalexi, EfthymiaE. and Satija, RahulR. (2018). Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nature biotechnology 36 411–420.
  * Büttner et al. (2017) [author] Büttner, MarenM., Miao, ZhichaoZ., Wolf, F AlexanderF. A., Teichmann, Sarah AS. A. and Theis, Fabian JF. J. (2017). Assessment of batch-correction methods for scRNA-seq data with a new test metric. BioRxiv 200345\.
  * Byrns et al. (2014) [author] Byrns, China NC. N., Pitts, Matthew WM. W., Gilman, Christy AC. A., Hashimoto, Ann CA. C. and Berry, Marla JM. J. (2014). Mice lacking selenoprotein P and selenocysteine lyase exhibit severe neurological dysfunction, neurodegeneration, and audiogenic seizures. Journal of Biological Chemistry 289 9662–9674.
  * Byun et al. (2018) [author] Byun, JinyoungJ., Schwartz, Ann GA. G., Lusk, ChristineC., Wenzlaff, Angela SA. S., De Andrade, MarizaM., Mandal, DiptasriD., Gaba, ColetteC., Yang, PingP., You, MingM., Kupert, Elena YE. Y. et al. (2018). Genome-wide association study of familial lung cancer. Carcinogenesis 39 1135–1140.
  * Caballero-Franco and Kissler (2016) [author] Caballero-Franco, CeliaC. and Kissler, StephanS. (2016). The autoimmunity-associated gene RGS1 affects the frequency of T follicular helper cells. Genes & Immunity 17 228–238.
  * Chen et al. (2019a) [author] Chen, JiaweiJ., Hou, ChongxianC., Zheng, ZongtaiZ., Lin, HanH., Lv, GuangzhaoG. and Zhou, DongD. (2019a). Identification of Secreted Phosphoprotein 1 (SPP1) as a Prognostic Factor in Lower-Grade Gliomas. World neurosurgery 130 e775–e785.
  * Chen et al. (2019b) [author] Chen, YSY., Guo, YY., Shen, HWH., Zhang, PP. and Chen, HH. (2019b). Screening function enrichment and related signaling pathways bioinformatics analysis of differentially expressed gene in glioma. Zhonghua yi xue za zhi 99 2311–2314.
  * de Sousa et al. (2017) [author] de Sousa, Juliana FerreiraJ. F., Torrieri, RaulR., Serafim, Rodolfo BortolozoR. B., Di Cristofaro, Luis Fernando MacedoL. F. M., Escanfella, Fábio DalbonF. D., Ribeiro, RodrigoR., Zanette, Dalila LucíolaD. L., Paçó-Larson, Maria LuisaM. L., da Silva Jr, Wilson AraujoW. A., Tirapelli, Daniela Pretti da CunhaD. P. d. C. et al. (2017). Expression signatures of DNA repair genes correlate with survival prognosis of astrocytoma patients. Tumor Biology 39 1010428317694552\.
  * Fraile et al. (2016) [author] Fraile, Julia MJ. M., Campos-Iglesias, DianaD., Rodríguez, FranciscoF., Español, YaizaY. and Freije, José MPJ. M. (2016). The deubiquitinase USP54 is overexpressed in colorectal cancer stem cells and promotes intestinal tumorigenesis. Oncotarget 7 74427\.
  * George et al. (2010) [author] George, Amee JA. J., Gordon, LaviniaL., Beissbarth, TimT., Koukoulas, IreneI., Holsinger, RM DamianR. D., Perreau, VictoriaV., Cappai, RobertoR., Tan, Seong-SengS.-S., Masters, Colin LC. L., Scott, Hamish SH. S. et al. (2010). A serial analysis of gene expression profile of the Alzheimer’s disease Tg2576 mouse model. Neurotoxicity research 17 360–379.
  * Grün et al. (2015) [author] Grün, DominicD., Lyubimova, AnnaA., Kester, LennartL., Wiebrands, KayK., Basak, OnurO., Sasaki, NobuoN., Clevers, HansH. and Van Oudenaarden, AlexanderA. (2015). Single-cell messenger RNA sequencing reveals rare intestinal cell types. Nature 525 251–255.
  * Haghverdi et al. (2018) [author] Haghverdi, LalehL., Lun, Aaron TLA. T., Morgan, Michael DM. D. and Marioni, John CJ. C. (2018). Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors. Nature biotechnology 36 421–427.
  * Hedlund and Deng (2018) [author] Hedlund, EvaE. and Deng, QiaolinQ. (2018). Single-cell RNA sequencing: technical advancements and biological applications. Molecular aspects of medicine 59 36–46.
  * Hide et al. (2018) [author] Hide, TakuichiroT., Komohara, YoshihiroY., Miyasato, YukoY., Nakamura, HideoH., Makino, KeishiK., Takeya, MotohiroM., Kuratsu, Jun-ichiJ.-i., Mukasa, AkitakeA. and Yano, ShigetoshiS. (2018). Oligodendrocyte progenitor cells and macrophages/microglia produce glioma stem cell niches at the tumor border. EBioMedicine 30 94–104.
  * Hubert and Arabie (1985) [author] Hubert, LawrenceL. and Arabie, PhippsP. (1985). Comparing partitions. Journal of classification 2 193–218.
  * Jia et al. (2019) [author] Jia, DanyunD., Lin, WeiW., Tang, HongliH., Cheng, YifanY., Xu, KaiweiK., He, YanshuY., Geng, WujunW. and Dai, QinxueQ. (2019). Integrative analysis of DNA methylation and gene expression to identify key epigenetic genes in glioblastoma. Aging (Albany NY) 11 5579\.
  * Johansson, Göransson and Westermark (2005) [author] Johansson, Fredrik KF. K., Göransson, HannaH. and Westermark, BengtB. (2005). Expression analysis of genes involved in brain tumor progression driven by retroviral insertional mutagenesis in mice. Oncogene 24 3896–3905.
  * Johnson, Li and Rabinovic (2007) [author] Johnson, W EvanW. E., Li, ChengC. and Rabinovic, ArielA. (2007). Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 8 118–127.
  * Kassambara and Mundt (2019) [author] Kassambara, AlboukadelA. and Mundt, FabianF. (2019). factoextra: Extract and Visualize the Results of Multivariate Data Analyses R package version 1.0.6.
  * Khalil et al. (2014) [author] Khalil, MM., Riedlbauer, BB., Langkammer, CC., Enzinger, CC., Ropele, SS., Stojakovic, TT., Scharnagl, HH., Culea, VV., Petzold, AA., Teunissen, CEC. et al. (2014). Cerebrospinal fluid transferrin levels are reduced in patients with early multiple sclerosis. Multiple Sclerosis Journal 20 1569–1577.
  * Kim et al. (2011) [author] Kim, Jong MyungJ. M., Song, Ji SunJ. S., Cho, Hyun HwaH. H., Shin, Keun KooK. K., Bae, Yong ChanY. C., Lee, Byung JuB. J. and Jung, Jin SupJ. S. (2011). Effect of the modulation of leucine zipper tumor suppressor 2 expression on proliferation of various cancer cells functions as a tumor suppressor. Molecular and cellular biochemistry 346 125–136.
  * Kiselev, Andrews and Hemberg (2019) [author] Kiselev, Vladimir YuV. Y., Andrews, Tallulah ST. S. and Hemberg, MartinM. (2019). Challenges in unsupervised clustering of single-cell RNA-seq data. Nature Reviews Genetics 20 273–282.
  * Kiselev et al. (2017) [author] Kiselev, Vladimir YuV. Y., Kirschner, KristinaK., Schaub, Michael TM. T., Andrews, TallulahT., Yiu, AndrewA., Chandra, TamirT., Natarajan, Kedar NK. N., Reik, WolfW., Barahona, MauricioM., Green, Anthony RA. R. et al. (2017). SC3: consensus clustering of single-cell RNA-seq data. Nature methods 14 483–486.
  * Kore and Abraham (2014) [author] Kore, Rajshekhar AR. A. and Abraham, Edathara CE. C. (2014). Inflammatory cytokines, interleukin-1 beta and tumor necrosis factor-alpha, upregulated in glioblastoma multiforme, raise the levels of CRYAB in exosomes secreted by U373 glioma cells. Biochemical and biophysical research communications 453 326–331.
  * Lange (2013) [author] Lange, KennethK. (2013). The MM algorithm. In Optimization 185–219. Springer.
  * Leite et al. (2014) [author] Leite, Fabíola AF. A., Lira, Regia CPR. C., Fedatto, Paola FP. F., Antonini, Sonir RRS. R., Martinelli Jr, Carlos EC. E., de Castro, MargaretM., Neder, LucianoL., Ramalho, Leandra NZL. N., Tucci Jr, SilvioS., Mastelaro, Maria JM. J. et al. (2014). Low expression of HLA-DRA, HLA-DPA1, and HLA-DPB1 is associated with poor prognosis in pediatric adrenocortical tumors (ACT). Pediatric blood & cancer 61 1940–1948.
  * Li and Dewey (2011) [author] Li, BoB. and Dewey, Colin NC. N. (2011). RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC bioinformatics 12 323\.
  * Lin, Troup and Ho (2017) [author] Lin, PeijieP., Troup, MichaelM. and Ho, Joshua WKJ. W. (2017). CIDR: Ultrafast and accurate clustering through imputation for single-cell RNA-seq data. Genome biology 18 59\.
  * Lindner et al. (2010) [author] Lindner, IngeI., Hemdan, Nasr YAN. Y., Buchold, MartinM., Huse, KlausK., Bigl, MarinaM., Oerlecke, IlkaI., Ricken, AlbertA., Gaunitz, FrankF., Sack, UlrichU., Naumann, AndreasA. et al. (2010). $\alpha$2-macroglobulin inhibits the malignant properties of astrocytoma cells by impeding $\beta$-catenin signaling. Cancer research 70 277–287.
  * M. Hajj et al. (2020) [author] M. Hajj, Glaucia NG. N., da Silva, Fernanda FF. F., de Bellis, BárbaraB., Lupinacci, Fernanda CSF. C., Bellato, Hermano MH. M., Cruz, Juvanier RJ. R., Segundo, Claudionor NCC. N., Faquini, Igor VI. V., Torres, Leuridan CL. C., Sanematsu, Paulo IP. I. et al. (2020). Aberrant expression of RSK1 characterizes high-grade gliomas with immune infiltration. Molecular oncology 14 159–179.
  * Ma et al. (2018) [author] Ma, ChangguoC., Gu, ChunjieC., Huo, YongxiaY., Li, XiaoyanX. and Luo, Xiong-JianX.-J. (2018). The integrated landscape of causal genes and pathways in schizophrenia. Translational psychiatry 8 1–14.
  * Maheras et al. (2018) [author] Maheras, Kathleen JK. J., Peppi, MarcelloM., Ghoddoussi, FarhadF., Galloway, Matthew PM. P., Perrine, Shane AS. A. and Gow, AlexanderA. (2018). Absence of claudin 11 in CNS myelin perturbs behavior and neurotransmitter levels in mice. Scientific reports 8 1–16.
  * Mangogna et al. (2019) [author] Mangogna, AlessandroA., Belmonte, BeatriceB., Agostinis, ChiaraC., Zacchi, PaolaP., Iacopino, Domenico GerardoD. G., Martorana, AnnaA., Rodolico, VitoV., Bonazza, DeborahD., Zanconati, FabrizioF., Kishore, UdayU. et al. (2019). Prognostic Implications of the Complement Protein C1q in Gliomas. Frontiers in immunology 10 2366\.
  * Metsu et al. (2014) [author] Metsu, SofieS., Rainger, Jacqueline KJ. K., Debacker, KimK., Bernhard, BirgittaB., Rooms, LiesbethL., Grafodatskaya, DariaD., Weksberg, RosannaR., Fombonne, EricE., Taylor, Martin SM. S., Scherer, Stephen WS. W. et al. (2014). A CGG-Repeat Expansion Mutation in ZNF 713 Causes FRA 7 A: Association with Autistic Spectrum Disorder in Two Families. Human mutation 35 1295–1300.
  * Mittelbronn et al. (2007) [author] Mittelbronn, MichelM., Simon, PeriklesP., Löffler, ChristianC., Capper, DavidD., Bunz, BenediktB., Harter, PatrickP., Schlaszus, HolgerH., Schleich, ArianeA., Tabatabai, GhazalehG., Goeppert, BenjaminB. et al. (2007). Elevated HLA-E levels in human glioblastomas but not in grade I to III astrocytomas correlate with infiltrating CD8+ cells. Journal of neuroimmunology 189 50–58.
  * Miyauchi et al. (2018) [author] Miyauchi, EisukeE., Furuta, TakuyaT., Ohtsuki, SumioS., Tachikawa, MasanoriM., Uchida, YasuoY., Sabit, HemragulH., Obuchi, WataruW., Baba, TomokoT., Watanabe, MichitoshiM., Terasaki, TetsuyaT. et al. (2018). Identification of blood biomarkers in glioblastoma by SWATH mass spectrometry and quantitative targeted absolute proteomics. PloS one 13\.
  * Moral-Morales et al. (2020) [author] Moral-Morales, DelD., González-Orozco, Juan CarlosJ. C., Capetillo-Velázquez, José MoisésJ. M., Piña-Medina, Ana GabrielaA. G., Camacho-Arroyo, IgnacioI. et al. (2020). The role of mPR$\delta$ and mPR$\varepsilon$ in human glioblastoma cells: expression, hormonal regulation, and possible clinical outcome. Hormones and Cancer 11 117–127.
  * Murakami et al. (2011) [author] Murakami, AkihiroA., Fukushima, ChikakoC., Yoshidomi, KeikoK., Sueoka, KotaroK., Nawata, ShugoS., Yokoyama, YoshihitoY., Tsuchida, ShigekiS., Ismail, EndomE., Al-Mulla, FahdF. and Sugino, NorihiroN. (2011). Suppression of carbonyl reductase expression enhances malignant behaviour in uterine cervical squamous cell carcinoma: Carbonyl reductase predicts prognosis and lymph node metastasis. Cancer letters 311 77–84.
  * Nicole and Pacary (2020) [author] Nicole, OlivierO. and Pacary, EmilieE. (2020). CaMKII$\beta$ in Neuronal Development and Plasticity: An Emerging Candidate in Brain Diseases. International Journal of Molecular Sciences 21 7272\.
  * O’Leary et al. (2016) [author] O’Leary, Nuala AN. A., Wright, Mathew WM. W., Brister, J RodneyJ. R., Ciufo, StacyS., Haddad, DianaD., McVeigh, RichR., Rajput, BhanuB., Robbertse, BarbaraB., Smith-White, BrianB., Ako-Adjei, DansoD. et al. (2016). Reference sequence (RefSeq) database at NCBI: current status, taxonomic expansion, and functional annotation. Nucleic acids research 44 D733–D745.
  * Petegrosso, Li and Kuang (2019) [author] Petegrosso, RaphaelR., Li, ZhuliuZ. and Kuang, RuiR. (2019). Machine learning and statistical methods for clustering single-cell RNA-sequencing data. Briefings in bioinformatics.
  * Petropoulos et al. (2016) [author] Petropoulos, SophieS., Edsgärd, DanielD., Reinius, BjörnB., Deng, QiaolinQ., Panula, Sarita PauliinaS. P., Codeluppi, SimoneS., Reyes, Alvaro PlazaA. P., Linnarsson, StenS., Sandberg, RickardR. and Lanner, FredrikF. (2016). Single-cell RNA-seq reveals lineage and X chromosome dynamics in human preimplantation embryos. Cell 165 1012–1026.
  * Peyrot and Price (2021) [author] Peyrot, Wouter JW. J. and Price, Alkes LA. L. (2021). Identifying loci with different allele frequencies among cases of eight psychiatric disorders using CC-GWAS. Nature genetics 53 445–454.
  * Picelli et al. (2014) [author] Picelli, SimoneS., Faridani, Omid RO. R., Björklund, Åsa KÅ. K., Winberg, GöstaG., Sagasser, SvenS. and Sandberg, RickardR. (2014). Full-length RNA-seq from single cells using Smart-seq2. Nature protocols 9 171\.
  * Pierson and Yau (2015) [author] Pierson, EmmaE. and Yau, ChristopherC. (2015). ZIFA: Dimensionality reduction for zero-inflated single-cell gene expression analysis. Genome biology 16 241\.
  * Prakadan, Shalek and Weitz (2017) [author] Prakadan, Sanjay MS. M., Shalek, Alex KA. K. and Weitz, David AD. A. (2017). Scaling by shrinking: empowering single-cell’omics’ with microfluidic devices. Nature Reviews Genetics 18 345\.
  * Reinius et al. (2016) [author] Reinius, BjörnB., Mold, Jeff EJ. E., Ramsköld, DanielD., Deng, QiaolinQ., Johnsson, PerP., Michaëlsson, JakobJ., Frisén, JonasJ. and Sandberg, RickardR. (2016). Analysis of allelic expression patterns in clonal somatic cells by single-cell RNA–seq. Nature genetics 48 1430\.
  * Risso et al. (2018) [author] Risso, DavideD., Perraudeau, FannyF., Gribkova, SvetlanaS., Dudoit, SandrineS. and Vert, Jean-PhilippeJ.-P. (2018). A general and flexible method for signal extraction from single-cell RNA-seq data. Nature communications 9 1–17.
  * Rousseau et al. (2006) [author] Rousseau, AudreyA., Nutt, Catherine LC. L., Betensky, Rebecca AR. A., Iafrate, A JohnA. J., Han, MoonjooM., Ligon, Keith LK. L., Rowitch, David HD. H. and Louis, David ND. N. (2006). Expression of oligodendroglial and astrocytic lineage markers in diffuse gliomas: use of YKL-40, ApoE, ASCL1, and NKX2-2. Journal of Neuropathology & Experimental Neurology 65 1149–1156.
  * Rousseeuw (1987) [author] Rousseeuw, Peter JP. J. (1987). Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. Journal of computational and applied mathematics 20 53–65.
  * Sawcer et al. (2011) [author] Sawcer, StephenS., Hellenthal, GarrettG., Pirinen, MattiM., Spencer, Chris CAC. C., Patsopoulos, Nikolaos AN. A., Moutsianas, LoukasL., Dilthey, AlexanderA., Su, ZhanZ., Freeman, ColinC., Hunt, Sarah ES. E. et al. (2011). Genetic risk and a primary role for cell-mediated immune mechanisms in multiple sclerosis. Nature 476 214\.
  * Shi et al. (2017) [author] Shi, YuY., Ping, Yi-FangY.-F., Zhou, WenchaoW., He, Zhi-ChengZ.-C., Chen, CongC., Zhang, LinL., Chen, LuL., Lan, XunX., Zhang, Xian-ChaoX.-C., Zhou, KaiK. et al. (2017). Tumour-associated macrophages secrete pleiotrophin to promote PTPRZ1 signalling in glioblastoma stem cells for tumour growth. Nature communications 8 1–17.
  * Sinmaz et al. (2016) [author] Sinmaz, NeseN., Nguyen, TinaT., Tea, FionaF., Dale, Russell CR. C. and Brilot, FabienneF. (2016). Mapping autoantigen epitopes: molecular insights into autoantibody-associated disorders of the nervous system. Journal of neuroinflammation 13 219\.
  * Sinreih et al. (2018) [author] Sinreih, MašaM., Knific, TamaraT., Thomas, PeterP., Grazio, Snježana FrkovićS. F. and Rižner, Tea LanišnikT. L. (2018). Membrane progesterone receptors $\beta$ and $\gamma$ have potential as prognostic biomarkers of endometrial cancer. The Journal of steroid biochemistry and molecular biology 178 303–311.
  * Stegle, Teichmann and Marioni (2015) [author] Stegle, OliverO., Teichmann, Sarah AS. A. and Marioni, John CJ. C. (2015). Computational and analytical challenges in single-cell transcriptomics. Nature Reviews Genetics 16 133–145.
  * Szulzewsky et al. (2015) [author] Szulzewsky, FrankF., Pelz, AndreasA., Feng, XiX., Synowitz, MichaelM., Markovic, DarkoD., Langmann, ThomasT., Holtman, Inge RI. R., Wang, XiX., Eggen, Bart JLB. J., Boddeke, Hendrikus WGMH. W. et al. (2015). Glioma-associated microglia/macrophages display an expression profile different from M1 and M2 polarization and highly express Gpnmb and Spp1. PloS one 10\.
  * R Core Team (2018) [author] R Core Team (2018). R: A Language and Environment for Statistical Computing R Foundation for Statistical Computing, Vienna, Austria.
  * Thyssen et al. (2006) [author] Thyssen, GregoryG., Li, Tzu-HueyT.-H., Lehmann, LynnL., Zhuo, MingM., Sharma, ManjuM. and Sun, ZijieZ. (2006). LZTS2 is a novel $\beta$-catenin-interacting protein and regulates the nuclear export of $\beta$-catenin. Molecular and cellular biology 26 8857–8867.
  * Tibshirani, Walther and Hastie (2001) [author] Tibshirani, RobertR., Walther, GuentherG. and Hastie, TrevorT. (2001). Estimating the number of clusters in a data set via the gap statistic. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 63 411–423.
  * Treutlein et al. (2014) [author] Treutlein, BarbaraB., Brownfield, Doug GD. G., Wu, Angela RA. R., Neff, Norma FN. F., Mantalas, Gary LG. L., Espinoza, F HernanF. H., Desai, Tushar JT. J., Krasnow, Mark AM. A. and Quake, Stephen RS. R. (2014). Reconstructing lineage hierarchies of the distal lung epithelium using single-cell RNA-seq. Nature 509 371–375.
  * Tung et al. (2017) [author] Tung, Po-YuanP.-Y., Blischak, John DJ. D., Hsiao, Chiaowen JoyceC. J., Knowles, David AD. A., Burnett, Jonathan EJ. E., Pritchard, Jonathan KJ. K. and Gilad, YoavY. (2017). Batch effects and the effective design of single-cell gene expression studies. Scientific reports 7 39921\.
  * Vallejos et al. (2017) [author] Vallejos, Catalina AC. A., Risso, DavideD., Scialdone, AntonioA., Dudoit, SandrineS. and Marioni, John CJ. C. (2017). Normalizing single-cell RNA sequencing data: challenges and opportunities. Nature methods 14 565\.
  * Venteicher et al. (2017) [author] Venteicher, Andrew SA. S., Tirosh, ItayI., Hebert, ChristineC., Yizhak, KerenK., Neftel, CyrilC., Filbin, Mariella GM. G., Hovestadt, VolkerV., Escalante, Leah EL. E., Shaw, McKenzie LM. L., Rodman, ChristopherC. et al. (2017). Decoupling genetics, lineages, and microenvironment in IDH-mutant gliomas by single-cell RNA-seq. Science 355 eaai8478.
  * Vidyarthi et al. (2019) [author] Vidyarthi, AurobindA., Agnihotri, TapanT., Khan, NargisN., Singh, SanpreetS., Tewari, Manoj KM. K., Radotra, Bishan DB. D., Chatterjee, DeepyanD. and Agrewala, Javed NJ. N. (2019). Predominance of M2 macrophages in gliomas leads to the suppression of local and systemic immunity. Cancer Immunology, Immunotherapy 68 1995–2004.
  * Wang and Bodovitz (2010) [author] Wang, DaojingD. and Bodovitz, StevenS. (2010). Single cell analysis: the new frontier in ‘omics’. Trends in biotechnology 28 281–290.
  * Wang, Li and Leng (2009) [author] Wang, HanshengH., Li, BoB. and Leng, ChenleiC. (2009). Shrinkage tuning parameter selection with a diverging number of parameters. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 71 671–683.
  * Wang et al. (2014) [author] Wang, JiajiaJ., Feng, XuemeiX., Bai, ZhouxianZ., Jin, Lee-WayL.-W., Duan, YongY. and Lei, HongxingH. (2014). Chromosome 19p in Alzheimer’s disease: When genome meets transcriptome. Journal of Alzheimer’s disease 38 245–250.
  * Wang et al. (2019) [author] Wang, Shao-MingS.-M., Lin, Hong-YiH.-Y., Chen, Yen-LinY.-L., Hsu, Tsung-IT.-I., Chuang, Jian-YingJ.-Y., Kao, Tzu-JenT.-J. and Ko, Chiung-YuanC.-Y. (2019). CCAAT/enhancer-binding protein delta regulates the stemness of glioma stem-like cells through activating PDGFA expression upon inflammatory stimulation. Journal of neuroinflammation 16 146\.
  * Xia et al. (2019) [author] Xia, ZhenkunZ., Ouyang, DengjieD., Li, QianyingQ., Li, MoyunM., Zou, QiongyanQ., Li, LunL., Yi, WenjunW. and Zhou, EnxiangE. (2019). The expression, functions, interactions and prognostic values of PTPRZ1: a review and bioinformatic analysis. Journal of Cancer 10 1663\.
  * Xiong et al. (2019) [author] Xiong, Dan-DanD.-D., Xu, Wen-QingW.-Q., He, Rong-QuanR.-Q., Dang, Yi-WuY.-W., Chen, GangG. and Luo, Dian-ZhongD.-Z. (2019). In silico analysis identified miRNA-based therapeutic agents against glioblastoma multiforme. Oncology reports 41 2194–2208.
  * Xu et al. (2020) [author] Xu, JunlinJ., Cai, LijunL., Liao, BoB., Zhu, WenW. and Yang, JiaLiangJ. (2020). CMF-Impute: an accurate imputation tool for single-cell RNA-seq data. Bioinformatics.
  * Xue et al. (2013) [author] Xue, ZhigangZ., Huang, KevinK., Cai, ChaochaoC., Cai, LingboL., Jiang, Chun-yanC.-y., Feng, YunY., Liu, ZhenshanZ., Zeng, QiaoQ., Cheng, LimingL., Sun, Yi EY. E. et al. (2013). Genetic programs in human and mouse early embryos revealed by single-cell RNA sequencing. Nature 500 593–597.
  * Yuan et al. (2018) [author] Yuan, JinzhouJ., Levitin, Hanna MendesH. M., Frattini, VeroniqueV., Bush, Erin CE. C., Boyett, Deborah MD. M., Samanamud, JorgeJ., Ceccarelli, MicheleM., Dovas, AthanassiosA., Zanazzi, GeorgeG., Canoll, PeterP. et al. (2018). Single-cell transcriptome analysis of lineage diversity in high-grade glioma. Genome medicine 10 1–15.
  * Zeiner et al. (2015) [author] Zeiner, Pia SP. S., Preusse, CorinnaC., Blank, Anna-EvaA.-E., Zachskorn, CorneliaC., Baumgarten, PeterP., Caspary, LixiL., Braczynski, Anne KA. K., Weissenberger, JakobJ., Bratzke, HansjürgenH., Reiß, SandyS. et al. (2015). MIF Receptor CD 74 is Restricted to Microglia/Macrophages, Associated with a M 1-Polarized Immune Milieu and Prolonged Patient Survival in Gliomas. Brain pathology 25 491–504.
  * Zeiner et al. (2019) [author] Zeiner, Pia SP. S., Preusse, CorinnaC., Golebiewska, AnnaA., Zinke, JennyJ., Iriondo, AneA., Muller, ArnaudA., Kaoma, TonyT., Filipski, KatharinaK., Müller-Eschner, MonikaM., Bernatz, SimonS. et al. (2019). Distribution and prognostic impact of microglia/macrophage subpopulations in gliomas. Brain Pathology 29 513–529.
  * Zhan et al. (2015) [author] Zhan, XinhuaX., Jickling, Glen CG. C., Ander, Bradley PB. P., Stamova, BoryanaB., Liu, DaZhiD., Kao, Patricia FP. F., Zelin, Mariko AM. A., Jin, Lee-WayL.-W., DeCarli, CharlesC. and Sharp, Frank RF. R. (2015). Myelin Basic Protein Associates with A$\beta$PP, A$\beta$ 1-42, and Amyloid Plaques in Cortex of Alzheimer’s Disease Brain. Journal of Alzheimer’s disease 44 1213–1229.
  * Zhang, Wu and Tian (2019) [author] Zhang, FengF., Wu, YuY. and Tian, WeidongW. (2019). A novel approach to remove the batch effect of single-cell data. Cell discovery 5 1–4.
  * Zhang et al. (2018) [author] Zhang, Da-LinD.-L., Qu, Li-WeiL.-W., Ma, LiangL., Zhou, Yong-ChunY.-C., Wang, Gui-ZhenG.-Z., Zhao, Xin-ChunX.-C., Zhang, ChenC., Zhang, Yan-FeiY.-F., Wang, MinM., Zhang, Mei-YingM.-Y. et al. (2018). Genome-wide identification of transcription factors that are critical to non-small cell lung cancer. Cancer letters 434 132–143.