基于图协同过滤的单细胞RNA测序数据填补

Imputation of scRNA-seq Data Based on Graph Collaborative Filtering

摘要: 单细胞RNA测序(Single-cell RNA Sequencing, scRNA-seq)技术能以单细胞的分辨率分析转录组数据，在生物学研究中展现出广泛的应用前景。然而技术问题会导致scRNA-seq数据存在部分基因表达缺失的情况，称之为零膨胀事件。这种情况严重阻碍了下游分析，故需要对scRNA-seq数据进行填补。本文提出了一种基于图协同过滤的单细胞RNA测序数据填补算法，为scRNA-seq分析提供了一个深度学习框架。它通过结构邻居对比的图协同过滤方法提取细胞特征表示和基因特征表示，并将两者的内积应用于零膨胀负二项分布自编码器来填补scRNA-seq数据。仿真实验结果验证了该算法在仿真数据集上的填补能力，且通过下游聚类分析实验表明该算法在公共真实数据集上细胞聚类的性能。

Abstract: Single-cell RNA sequencing (scRNA-seq) technology can analyze transcriptome data at the single-cell level and is widely used in biology. However, technical issues can lead to missing gene expression in scRNA-seq data, which is called zero-inflation event. This situation seriously hinders downstream analysis, so it is necessary to impute the scRNA-seq data. This article proposes an imputation algorithm of scRNA-seq data based on graph collaborative filtering, providing a deep learning framework for scRNA-seq analysis. It extracts cell feature representations and gene feature representations through the graph collaborative filtering method of comparing structural neighbors, and applies the inner product of the two to the zero-inflated negative binomial distribution autoencoder to impute scRNA-seq data. The simulation experiment results have verified the imputation ability of the algorithm on the simulation dataset, and downstream clustering analysis experiments have shown the performance of the algorithm on cell clustering on public real datasets.

1. 引言

单细胞RNA测序(Single-cell RNA Sequencing, scRNA-seq) [1] 技术已成为生物信息学的研究热点。scRNA-seq技术在单细胞水平上分析转录组数据，且应用广泛 [2] [3] [4] 。

尽管scRNA-seq技术在单细胞分析中具有潜力，但其应用仍受到技术噪声和实验条件的限制。例如，由于RNA输入不足或细胞测序深度不足等技术或实验条件限制，可能会出现部分基因表达数据缺失的情况，即所谓的零膨胀事件或dropout事件。这些缺失值可能会导致重要生物学信息的丢失，并对scRNA-seq数据的下游分析造成阻碍。采取措施来估算或推断这些缺失值是处理scRNA-seq数据集的一个不可或缺的步骤，也就是scRNA-seq数据填补。

针对scRNA-seq数据的缺失值填补问题，目前研究人员已经开发了许多填补方法。例如，精确单细胞填补(Single-cell Impute, scImpute) [5] 通过Gamma-正态分布混合模型估计哪些值受到dropout的影响，最后通过借用其他相似细胞中相同基因的信息来估算。基于细胞Markov亲和图的填补(Markov Affinity-based Graph Imputation of Cells, MAGIC) [6] 基于热扩散的思想并对相似细胞中的信息加权来进行估算dropout。而自适应阈值低秩近似填补(Adaptively Thresholded Low-rank Approximation, ALRA) [7] 通过观察到的基因表达矩阵进行低秩矩阵补全，再基于奇异值分解求解进而填补。除了基于统计的方法，目前针对scRNA-seq数据的高维度、高稀疏性，基于深度学习的方法具备有效性和高效性。深度计数自编码器去噪(Deep Count Autoencoder Network, DCA) [8] 将基因表达分布建模为负二项(Negative Binomial, NB)分布或零膨胀负二项(Zero-inflated Negative Binomial, ZINB)分布，通过自编码器学习到的分布参数进而预测去噪后的基因表达矩阵。而单细胞变分推断(Single-cell Variational Inference, scVI) [9] 通过变分自编码器来指定ZINB分布进行填补。单细胞图神经网络填补(Single-cell Graph Neural Network, scGNN) [10] 通过特征自编码器学习并构建细胞图，利用图自编码器聚合细胞间关系，最后基于特征自编码器重构表达谱。

针对scRNA-seq数据填补问题，本文提出一种基于图协同过滤的单细胞RNA测序数据填补(Imputation of scRNA-seq Data based on Graph Collaborative Filtering, scGCF)算法。通过图协同过滤分别获得细胞特征表示和基因特征表示，再将细胞特征表示和基因特征表示交互，进而输入基于ZINB分布的自编码器重建scRNA-seq数据的表达谱。

2. 模型介绍

scGCF模型主要由图协同过滤框架、ZINB分布自编码器这两个部分组成，其模型示意图如图1所示。给定基因表达矩阵 $X={\left(x_{ij}\right)}_{m\times n}\in {ℝ}^{m\times n}$ ，细胞的集合 $\mathcal{C}=\left\{c\right\}$ 且 $\left|\mathcal{C}\right|=m$ ，基因的集合 $\mathcal{G}=\left\{g\right\}$ 且 $\left|\mathcal{G}\right|=n$ ，细胞基因交互矩阵是 $A={\left(a_{ij}\right)}_{m\times n}\in {ℝ}^{m\times n}$ ，且 $a_{ij}=\{\begin{array}{l}1,\text{if}{x}_{ij}>0\\ 0,\text{if}{x}_{ij}=0\end{array}$ ，构图如下：

$\mathcal{T}=\left\{\mathcal{V},\mathcal{E}\right\}$ (1)

其中 $\mathcal{V}=\left\{\mathcal{C}\cup \mathcal{G}\right\}$ 表示节点集， $\mathcal{E}=\left\{\left(c,g\right)|c\in \mathcal{C},g\in \mathcal{G},A_{cg}=1\right\}$ 表示边集。

Figure 1. Schematic diagram of scGCF’s model structure

图1. scGCF模型结构示意图

2.1. 数据预处理

针对scRNA-seq数据的原始计数矩阵，本文采取了预处理操作以避免低质量数据对后续分析的影响。本文移除了在极少数(少于3个)细胞中表达的基因以及在少数(少于200个)基因中表达的细胞，得到计数矩阵 $X={\left(x_{ij}\right)}_{m\times n}\in {ℝ}^{m\times n}$ ，其中 $m,n$ 分别为过滤后的细胞和基因的数量。

对细胞i，大小因子为：

$s_i={\displaystyle {\sum }_{j=1}^n{x}_{ij}}/media{n}_{i=1,\cdots ,m}\left({\displaystyle {\sum }_{j=1}^n{x}_{ij}}\right),$ (2)

再作大小因子归一化来减轻不同测序深度可能带来的影响，后经log转换、伪计数加1后，得到：

$X^{\prime }=\log \left(G_c^{-1}X+I\right),$ (3)

其中 $I$ 为全一元素的矩阵， $G_c=\mathrm{diag}\left(s_1,s_2,\cdots ,s_m\right)\in {ℝ}^{m\times m}$ 且 $\mathrm{diag}\left(s_1,s_2,\cdots ,s_m\right)$ 是一个以 $s_1,s_2,\cdots ,s_m$ 作为对角线的对角矩阵。

对基因j，定义基因因子：

$g_j=\underset{i}{\max }{x^{\prime }}_{ij},$ (4)

同理全部基因因子构成的矩阵为 $G_g=\mathrm{diag}\left(g_1,g_2,\cdots ,g_n\right)\in {ℝ}^{n\times n}$ 。矩阵 $X$ 中每个细胞作为一个样本，其中的细胞序号i和所有基因的序号作为输入，样本中基因表达作为标签，输入模型。矩阵 $G_c$ 和 $G_g$ 则用于最后的scRNA-seq数据填补。

2.2. 图协同过滤框架

对观察到的细胞和基因之间的交互进行建模，通过在图 $\mathcal{T}$ 上应用传播和预测函数来生成细胞和基因特征表示：

$e_c^{\left(l+1\right)}=\text{ReLU}\left(\underset{i\in N_c}{{{\displaystyle \sum }}^{\text{}}}\frac{1}{\sqrt{\left|N_c\right|\left|N_i\right|}}{e}_i^{\left(l\right)}\right),e_g^{\left(l+1\right)}=\text{ReLU}\left(\underset{j\in N_g}{{{\displaystyle \sum }}^{\text{}}}\frac{1}{\sqrt{\left|N_g\right|\left|N_j\right|}}{e}_j^{\left(l\right)}\right),$ (5)

其中K是图神经网络(Graph Neural Network, GNN)的层数， $e_c^{\left(0\right)}\in {ℝ}^{k\times 1}$ 和 $e_g^{\left(0\right)}\in {ℝ}^{k\times 1}$ 分别是Xavier初始化后的细胞c特征表示以及基因g特征表示，参数k是细胞(基因)特征表示中的特征维度。 $e_c^{\left(l+1\right)}$ 和 $e_g^{\left(l+1\right)}$ 分别是经过 $l+1$ 层GNN后的细胞c特征表示及其基因g特征表示，且 $l\left(0\le l\le K-1\right)$ 。 $N_c$ 和 $\left|N_c\right|$ 分别表示细胞c在 $\mathcal{T}$ 上的邻居及其邻居个数。经过所有的GNN层后，采用加权和函数来得到细胞特征表示和基因特征表示：

$e_c=\frac{1}{K+1}\underset{l=0}{\overset{K}{{{\displaystyle \sum }}^{\text{}}}}{e}_c^{\left(l\right)},e_g=\frac{1}{K+1}\underset{l=0}{\overset{K}{{{\displaystyle \sum }}^{\text{}}}}{e}_g^{\left(l\right)},$ (6)

并采用 ${\xi }_{cg}=e_c^{\mathrm{T}}\cdot e_g$ 来预测基因g在细胞c中的表达量。

鉴于交互图是二分图，图协同过滤框架在图上进行偶数次GNN的传播，会自然聚合同质结构邻居的信息。本文对每个细胞及其结构邻居作对比，将细胞本身的表示和偶数层GNN相应输出的表示视为一对，基于InfoNCE来最小化每一对之间的距离：

$L^{\mathcal{C}}=\underset{c\in \mathcal{C}}{{{\displaystyle \sum }}^{\text{}}}-\log \frac{\exp \left[\left(e_c^{\left(e\right)}\cdot e_c^{\left(0\right)}\right)/\tau \right]}{\underset{\tilde{c}\in \mathcal{C}}{{{\displaystyle \sum }}^{\text{}}}\exp \left[\left(e_c^{\left(e\right)}\cdot e_{\tilde{c}}^{\left(0\right)}\right)/\tau \right]},$ (7)

其中 $\tau$ 代表softmax的温度超参，e是一个偶数， $e_c^{\left(e\right)}$ 和 $e_c^{\left(0\right)}$ 分别为第e层GNN和初始的细胞表示。同样的，可以得到针对基因的结构对比损失：

$L^{\mathcal{G}}=\underset{g\in \mathcal{G}}{{{\displaystyle \sum }}^{\text{}}}-\log \frac{\exp \left[\left(e_g^{\left(e\right)}\cdot e_g^{\left(0\right)}\right)/\tau \right]}{\underset{\tilde{g}\in \mathcal{G}}{{{\displaystyle \sum }}^{\text{}}}\exp \left[\left(e_g^{\left(e\right)}\cdot e_{\tilde{g}}^{\left(0\right)}\right)/\tau \right]},$ (8)

其中 $e_g^{\left(e\right)}$ 和 $e_g^{\left(0\right)}$ 分别为第e层GNN和初始的基因表示。总结构对比损失定义如下：

$L_S=L^{\mathcal{C}}+\alpha L^{\mathcal{G}}$ (9)

其中 $\alpha$ 是一个超参数，用于平衡权重。

2.3. ZINB分布自编码器

考虑到scRNA-seq的计数数据表现为高度稀疏和过度分散，假设其分布为零膨胀负二项分布：

$\mathrm{ZINB}\left(x|\pi ,\mu ,\theta \right)=\pi \cdot I_0\left(x\right)+\left(1-\pi \right)\cdot \mathrm{NB}\left(x|\mu ,\theta \right),$ (10)

其中负二项分布 $\mathrm{NB}\left(x|\mu ,\theta \right)=\frac{\Gamma \left(x+\theta \right)}{x!\Gamma \left(\theta \right)}{\left(\frac{\theta }{\theta +\mu }\right)}^{\theta }{\left(\frac{\mu }{\theta +\mu }\right)}^x$ 表征scRNA-seq数据的计数分布， $\mu ,\theta$ 表示均值和散度， $\pi$ 表示真实的基因表达值被观测为0的概率， $I_0\left(x\right)$ 为示性函数。

基于scRNA-seq数据的ZINB分布特征来设置自解码器。基于ZINB分布的负对数似然损失函数为：

$L_{ZINB}=-\mathrm{ZINB}\left(x|\pi ,\mu ,\theta \right),$ (11)

以上述函数作为损失函数进行解码来模拟scRNA-seq数据的分布，进而得到三个输出层 $\Pi ,{\rm M},\Theta$ ，它们分别代表dropout事件概率、NB分布的均值和散度，则有下式：

$\left\{\hat{\Pi },\widehat{{\rm M}},\hat{\Theta }\right\}=\underset{\Pi ,{\rm M},\Theta }{\arg \max }\mathrm{ZINB}\left(X|\Pi ,{\rm M},\Theta \right)=\underset{\Pi ,{\rm M},\Theta }{\arg \max }{\displaystyle \underset{c=1}{\overset{m}{\prod }}{\displaystyle \underset{g=1}{\overset{n}{\prod }}ZINB}}\left(x_{cg}|{\pi }_{cg},{\mu }_{cg},{\theta }_{cg}\right),$ (12)

其中m、n分别表示细胞数和基因数。

将提出的总结构邻居对比学习损失与ZINB分布重构损失结合，并添加GNN层中表示的2-范数惩罚损失，将算法的训练损失定义为：

$L=L_{ZINB}+{\lambda }_1{L}_S+{\lambda }_2{\Vert \bar{e}\Vert }_2,$ (13)

其中 ${\lambda }_1$ 和 ${\lambda }_2$ 是控制所提出的总结构邻居对比学习损失和正则化项的权重的超参数， $\bar{e}$ 表示GNN层中细胞和基因表示的参数集合。

连接ZINB分布自编码器后的网络架构表示如下：

$\{\begin{array}{l}H=ReLU\left(\varphi W^{\text{T}}+b\right),\\ B=\text{BatchNorm}\left(H\right),\\ M=G_c\left[\exp \left(\mathrm{sigmoid}\left(B{W}_{{\rm M}}^T+b_{{\rm M}}\right)\times G_g\right)-1\right],\\ \Pi =\mathrm{sigmoid}\left(B{W}_{\Pi }^T+b_{\Pi }\right),\\ \Theta =\mathrm{softplus}\left(\left(B{W}_{\Theta }^T+b_{\Theta }\right)\times G_g\right),\end{array}$ (14)

其中 $\varphi ={\left({\xi }_{cg}\right)}_{m\times n}\in {ℝ}^{m\times n}$ 为经过K (默认为2)个GNN层后初步估计的基因表达矩阵，r为隐藏层的神经元数。 $W\in {ℝ}^{r\times n}$ 和 $b\in {ℝ}^{m\times r}$ 分别是从交互层到隐藏层的权重和偏置， $H\in {ℝ}^{m\times r}$ 是隐藏层表示。 $\text{BatchNorm}(\cdot )$ 表示批次归一化(Batch Normalization, BN)层的函数， ${\rm B}\in {ℝ}^{m\times r}$ 是BN层输出的表示。 ${\rm M}\in {ℝ}^{m\times n}$ 、 $\Pi \in {ℝ}^{m\times n}$ 、 $\Theta \in {ℝ}^{m\times n}$ 分别指ZINB分布的三个输出层输出， $W_{{\rm M}}\in {ℝ}^{n\times r}$ 、 $b_{{\rm M}}\in {ℝ}^{m\times n}$ 分别是从神经元为r的BN层映射回神经元为n的参数网络层 $M$ 时的权重和偏置。 $W_{\Pi }\in {ℝ}^{n\times r}$ 、 $b_{\Pi }\in {ℝ}^{m\times n}$ 和 $W_{\Theta }\in {ℝ}^{n\times r}$ 、 $b_{\Theta }\in {ℝ}^{m\times n}$ 同理。指数函数和softplus函数分别应用于 $M$ 和 $\Theta$ 以确保它们的非负性，sigmoid激活函数则限制 $\Pi$ 的取值范围为。 $M$ 作为最终预测的基因表达矩阵。

3. 数值实验

3.1. scRNA-seq数据集

为了评估scRNA-seq数据的填补效果，使用R包Splatter [11] 生成scRNA-seq仿真数据。固定仿真数据中真正零表达的比例为35%，分别设计含有dropout时零值比例为85%、90%、95%的仿真观测矩阵及其对应不含dropout的真实矩阵。利用这三个稀疏度不同的仿真数据集进行仿真实验。

为了测试scGCF在聚类分析中的能力，采用真实scRNA-seq数据集：PBMC 数据集 [12] 、Zeisel数据集 [4] 、Human kidney数据集 [13] 和Adam数据集 [14] 。PBMC数据集是由10x基因组学平台测序得到的人类外周血单核细胞数据集，Human kidney 数据集是由10x基因组学平台测序得到的人类肾脏细胞数据集，两者均可从10x基因组学网站下载。Zeisel数据集是由Illumina测序平台得到的3005个来自小鼠体感皮层和海马体区域的细胞组成，Adam数据集则由3660个小鼠肾脏细胞数据组成且由Drop-seq平台测序，两者均可在基因表达数据库下载。

3.2. 对比算法及评价指标

本文以观察到的表达矩阵作基准，记为observed，并选择了六种对比算法scImpute [5] 、MAGIC [6] 、DCA [8] 、ALRA [7] 、scVI [9] 、scGNN [10] 。scImpute和MAGIC作为非深度学习对比算法。而由于算法涉及到图协同网络和ZINB自编码器，故选取ALRA、DCA、scVI和scGNN进行对比。所有对比算法均采用其默认参数来执行填补操作。

为了评价填补的效果，本文从仿真数据基因表达恢复实验和聚类分析实验来对scGCF进行评价。在填补实验中，通过Splatter [11] 包模拟得到一个设置细胞数量m、基因量数n、零表达值比例等参数的真实表达矩阵 $Y={\left(y_{ij}\right)}_{m\times n}$ 和对应含有dropout的观测基因表达矩阵 $Y^{\prime }={\left({y^{\prime }}_{ij}\right)}_{m\times n}$ 。再对 $Y^{\prime }$ 进行填补得到矩阵 $\hat{Y}={\left({\hat{y}}_{ij}\right)}_{m\times n}$ ，将计算 $\hat{Y}$ 与真实基因表达 $Y$ 之间的平均绝对误差(Mean Absolute Error, MAE)、均方根误差(Root Mean Squared Error, RMSE)、Pearson相关系数(Pearson Correlation Coefficient, PCC)这三个评价指标以衡量填补后的数据在数值和结构上与真实数据的接近程度。MAE和RMSE都表示填补误差，数值越小表征填补越精确，其计算公式如下：

$\text{MAE}=\frac{1}{mn}{\displaystyle \underset{i}{\sum }{\displaystyle \underset{j}{\sum }\left|y_{ij}-{\hat{y}}_{ij}\right|}},\text{RMSE}={\left(\frac{1}{mn}{\displaystyle \underset{i}{\sum }{\displaystyle \underset{j}{\sum }{\left(y_{ij}-{\hat{y}}_{ij}\right)}^2}}\right)}^{\frac{\text{1}}{\text{2}}},$ (15)

PCC表示填补后数据与原始数据之间的相关性，越接近1则相关性越强，其计算公式如下：

$\mathrm{PCC}=\frac{{\displaystyle \underset{i}{\sum }{\displaystyle \underset{j}{\sum }\left(y_{ij}-y_{mean}\right)\left({\hat{y}}_{ij}-{\hat{y}}_{mean}\right)}}}{\sqrt{{\displaystyle \underset{i}{\sum }{\displaystyle \underset{j}{\sum }{\left(y_{ij}-y_{mean}\right)}^2}}}\sqrt{{\displaystyle \underset{i}{\sum }{\displaystyle \underset{j}{\sum }{\left({\hat{y}}_{ij}-{\hat{y}}_{mean}\right)}^2}}}},$ (16)

其中 $y_{mean}$ 为真实基因表达矩阵 $Y$ 的均值， ${\hat{y}}_{mean}$ 表示填补后矩阵 $\hat{Y}$ 的均值。

聚类分析实验通过将细胞聚类后得到簇标签，将其和真实的细胞类型对比，以衡量聚类准确性。聚类分析的指标为可调节Rand系数(Adjusted Rand Index, ARI)和标准化互信息(Normalized Mutual Information, NMI)。ARI的定义如下：

$ARI\left(c,\hat{c}\right)\text{=}\left[{\displaystyle \underset{ij}{\sum }\left(\begin{array}{l}{n}_{ij}\\ 2\end{array}\right)}-\left({\displaystyle \underset{i}{\sum }\left(\begin{array}{l}{a}_i\\ 2\end{array}\right)}-{\displaystyle \underset{j}{\sum }\left(\begin{array}{l}{b}_j\\ 2\end{array}\right)}\right)/\left(\begin{array}{l}r\\ 2\end{array}\right)\right]/\left[\frac{1}{2}\left({\displaystyle \underset{i}{\sum }\left(\begin{array}{l}{a}_i\\ 2\end{array}\right)}+{\displaystyle \underset{j}{\sum }\left(\begin{array}{l}{b}_j\\ 2\end{array}\right)}\right)-\left({\displaystyle \underset{i}{\sum }\left(\begin{array}{l}{a}_i\\ 2\end{array}\right)}-{\displaystyle \underset{j}{\sum }\left(\begin{array}{l}{b}_j\\ 2\end{array}\right)}\right)/\left(\begin{array}{l}r\\ 2\end{array}\right)\right],$ (17)

其中 $N={\left(n_{ij}\right)}_{r\times k}$ 为混淆矩阵，r为真实细胞类型数，k为聚类簇数，且 $a_i={\displaystyle \underset{j}{\sum }{n}_{ij}}$ ， $b_j={\displaystyle \underset{i}{\sum }{n}_{ij}}$ ， $c,\hat{c}$ 。

分别为真实细胞标签和聚类簇标签。而NMI的定义如下：

$NMI=\frac{2I\left(U,V\right)}{H\left(U\right)+H\left(V\right)},$ (18)

其中 $I\left(U,V\right)$ 为互信息

$I\left(U,V\right)={\displaystyle \underset{i=1}{\overset{r}{\sum }}{\displaystyle \underset{j=1}{\overset{k}{\sum }}\frac{\left|u_i\cap v_j\right|}{m}\log \frac{m\left|u_i\cap v_j\right|}{\left|u_i\right|\times \left|v_j\right|}}},$ (19)

其中 $u_i$ 、 $v_j$ 分别为真实细胞的标签中属于第i类的集合和预测的簇标签中属于第j类的集合，而交叉熵 $H\left(U\right)$ 、 $H\left(V\right)$ 的公式为：

$H\left(U\right)=-{\displaystyle \underset{i=1}{\overset{r}{\sum }}\frac{u_i}{m}\log \frac{u_i}{m}},H\left(V\right)=-{\displaystyle \underset{j=1}{\overset{k}{\sum }}\frac{v_j}{m}\log \frac{v_j}{m}},$ (20)

ARI和NMI的数值越接近1意味着聚类效果越理想。

3.3. 仿真实验

为了验证填补效果，进行了仿真数据填补实验。本文设定将观测矩阵与其真实表达矩阵对比并计算MAE、RMSE以及PCC指标，记为observed方法。接着运用了scImpute、MAGIC、DCA、ALRA、scVI、scGNN和scGCF这七种不同的算法，对仿真的观测矩阵进行了填补操作，将这些填补后的矩阵与对应不同稀疏性的真实表达矩阵进行了对比，同样计算了三个指标。各算法计算获得的三个指标结果见图2。当零表达比例达到85%、90%和95%时，scGCF算法展现出了显著的优势，不仅取得了更高的PCC指标，还实现了更低的RMSE。特别是在零表达比例高达90%和95%的情况下，scGCF算法在MAE上表现超越了observed和其他6种填补算法。在零值比例为85%时scGCF算法的MAE值也仅次于DCA。这些实验结果充分表明scGCF能比其他对比算法更好地恢复基因表达矩阵。

(a) (b) (c)

Figure 2. Comparison of metrics (a) MAE; (b) RMSE; (c) PCC among methods on simulated data

图2. 仿真数据上各算法的指标比较(a) MAE值；(b) RMSE值；(c) PCC值

3.4. 聚类分析实验

为了研究scGCF用于聚类分析的性能，考虑PBMC数据集、Zeisel数据集、Human kidney数据集和Adam数据集。聚类时采用Louvain [15] 聚类，且聚类的分辨率设置为1，并对ARI和NMI指标进行计算。

对于observed、scImpute、MAGIC和ALRA算法，依次将观测矩阵或填补矩阵进行PCA将矩阵上的基因维度缩减至50，最后进行聚类。对于DCA、scVI、scGNN和scGCF算法，将对应算法得到的细胞隐层表示用于聚类。对于四个完整的数据集(PBMC、Zeisel、Human Kidney和Adam)，通过Louvain聚类获得的指标结果如图3所示。在这四个数据集上scGCF算法取得了更高的ARI值。同时scGCF在Zeisel数据集上取得了第二高的NMI值，在其余的三个数据集上取得了更高的NMI值。整体上，与其他对比算法相比scGCF在完整数据集上取得了较好的聚类效果。

为了评估这些算法对细胞聚类的稳健性，对这些完整的数据集使用了采样策略。随机对数据集中的每类细胞采样其各自的95%，并组合构成一个子样本集。该采样过程在一个数据集上独立执行了十次，得到了十个子样本，再将算法应用于这十个子样本聚类而得到聚类指标。每个数据集上得到的ARI值以及NMI值分别绘制为图4和图5中的箱线图。

由图4可知在Zeisel、Human kidney和Adam这三个数据集上scGCF的ARI值更高且更稳定，在PBMC数据集的ARI值上scGCF和observed均高于其他对比算法。由图5可知，scGCF在PBMC和Adam这两个数据集上的NMI值表现更好，而在Zeisel数据集上的NMI值仅低于scGNN，在Human kidney数据

Figure 3. (a) Barplot of ARI values obtained by clustering using methods on complete datasets; (b) Barplot of NMI values obtained by clustering using methods on complete datasets

图3. (a) 完整数据集上各算法聚类得到的ARI条形图；(b) 完整数据集上各算法聚类得到的NMI条形图

Figure 4. Boxplots of ARI values obtained by clustering using methods on 10 subsamples of the dataset. (a) PBMC dataset; (b) Zeisel dataset; (c) Human kidney dataset; (d) Adam dataset

图4. 在数据集10个子样本上应用各算法聚类而获得的ARI箱线图。(a) PBMC数据集；(b) Zeisel数据集；(c) Human kidney数据集；(d) Adam数据集

Figure 5. Boxplots of NMI values obtained by clustering using methods on 10 subsamples of the dataset. (a) PBMC dataset; (b) Zeisel dataset; (c) Human kidney dataset; (d) Adam dataset

图5. 在数据集10个子样本上应用各算法聚类而获得的NMI箱线图。(a) PBMC数据集；(b) Zeisel数据集；(c) Human kidney数据集；(d) Adam数据集

集的NMI值仅低于DCA。结合图3~图5的实验结果，与其他对比算法相比scGCF具有较好的聚类精度和鲁棒性。

4. 结论

综上，本文提出了一种新的scRNA-seq数据填补算法scGCF，它利用图协同过滤框架聚合得到细胞特征表示和基因特征表示，并使用基于ZINB分布的自编码器来填补scRNA-seq数据。在仿真数据集上验证了scGCF算法的有效性，并评估了其填补能力。同时在四个真实数据集上进行了聚类分析，表明scGCF良好的细胞聚类能力。scGCF不仅在scRNA-seq下游聚类分析上增强了可靠性，还拓展了scRNA-seq数据的填补算法领域。