CN111564183B - 融合基因本体和神经网络的单细胞测序数据降维方法 - Google Patents

Abstract

本发明提供了一种融合基因本体和神经网络的单细胞测序数据降维方法。首先，提取基因本体术语作为深层生物信息先验知识；接着，提取细胞之间的Must‑Link约束作为细胞层面上的先验知识；然后，将上述两种先验知识与自编码器模型结合，提出simGOAE模型；最后，根据simGOAE模型对单细胞测序数据进行训练降维。本发明提出的simGOAE模型不仅能适应大样本数据集的训练，还能更好地挖掘细胞的生物信息，实现更好的单细胞测序数据降维效果。

Description

融合基因本体和神经网络的单细胞测序数据降维方法

技术领域

本发明属生物信息处理技术领域，具体涉及一种融合基因本体和神经网络的单细胞测序数据降维方法。

背景技术

随着高通量测序技术的发展，科学家们提出单细胞测序技术并广泛应用于转录组学研究。该技术测定单个细胞的序列信息，提供更高的差异分辨率，使人们更好地理解单个细胞在其微环境中的功能。单细胞测序技术的出现为人类在生物信息学领域中进一步研究细胞功能、细胞差异性表达提供了可能性。单细胞测序数据(scRNA-seq) 最重要的应用之一是对细胞聚类，聚类结果有助于鉴定新的细胞类型、分析细胞功能、细胞分化轨迹。基于单细胞测序数据的聚类问题中主要挑战是高维度(约10000– 30000个基因)，降维技术能提高细胞聚类的准确性与可靠性，为数据分析提供保证。

迄今为止，已经提出了许多方法来消除噪声的影响并减小数据维度，根据是否使用神经网络模型大致可分为两类。在未使用神经网络的方法中，最简单且常用的方法是主成分分析(PCA)，假设数据是正态分布的，则PCA使用正交变换将一组可能相关的变量的观测值转换为一组线性不相关变量的值，这些值被称为主成分。但是，对于scRNA-seq数据集，它们并不是完全线性可分离的。文献“Maaten L,Hinton G. Visualizing data using t-SNE[J].Journal of machine learning research,2008,9(Nov): 2579-2605”提出了利用T分布随机邻居嵌入(t-SNE)非线性降维技术。t-SNE基于高斯内核，将高维数据转换为低维空间。但是，它通常将多维数据映射到适合人类观察的二维或三维。因此，它总是伴随着维度数目限制。

由于神经网络在其他领域的成功应用，越来越多的使用神经网络模型对单细胞测序数据进行降维处理的方法被提出。如Lin等人在文献“Lin C,Jain S,Kim H,et al.Using neural networks for reducing the dimensions of single-cell RNA-Seq data[J].Nucleic acids research,2017,45(17):e156-e156.”中提出了一种有监督的基于神经网络生成 scRNA-seq数据的低维表示方法。Lin等人提出的NN模型将神经网络与蛋白质相互作用(PPI)网络结合起来，对大量细胞进行分类，取得了较好地结果。此外，还有基于无监督的神经网络模型的，如自编码器。使用对称的编码器与解码器，自编码器模型无需人工提取特征，可以直接实现对高维数据的降维处理。Eraslan等人在文献“Eraslan G,Simon LM,Mircea M,et al.Single-cell RNA-seq denoising using a deep count autoencoder[J].Nature communications,2019,10(1):390.”中提出了深度计数自动编码器 (DCA)，该模型以基于零扩展的负二项式分布的损失函数代替了传统的MSE损失函数，实现了对高维单细胞数据的降维处理。Lopez等人在文献“Lopez R,Regier J,Cole M B,et al.Deepgenerative modeling for single-cell transcriptomics[J].Nature methods, 2018,15(12):1053.”提出了单细胞数据变分自编码器模型(scVI)，这是一种用于估算和拟合这些因素的深层变分自编码器模型，该模型已应用于一系列单细胞测序数据降维及其他基本分析任务并获得了非凡的结果。

虽然上述方法对单细胞测序数据降维取得了一定的效果，但其缺点在于：一方面，没有充分利用现有的生物信息知识，解释性不强，需要降维及其他下游分析才能获得细胞的生物信息；另一方面，现如今的单细胞测序数据降维方法效果一般，且对不同测序平台的单细胞测序数据有一定局限性。

发明内容

为了克服现有技术的不足，本发明提供一种融合基因本体和神经网络的单细胞测序数据降维方法。首先，提取基因本体术语作为深层生物信息先验知识；接着，提取细胞之间的Must-Link约束作为细胞层面上的先验知识；然后，将上述两种先验知识与自编码器模型结合，提出融合基因本体术语相似性的自编码器(simGOAE)模型；最后，根据simGOAE模型对单细胞测序数据进行训练降维。本发明提出的simGOAE模型不仅能适应大样本数据集的训练，还能更好地挖掘细胞的生物信息，实现更好的单细胞测序数据降维效果。

一种融合基因本体和神经网络的单细胞测序数据降维方法，其特征在于步骤如下：

步骤一，数据预处理：删除单细胞测序数据中基因表达小于3个细胞的基因；删除单细胞测序数据中在基因本体术语注释数据库中没有表达的基因；对单细胞测序数据在基因上进行z-score标准化；

步骤二，筛选基因本体术语：删除冗余本体术语和低表达本体术语；所述的冗余本体术语是指特异性分数在[0.5,1)范围内的两个基因本体术语中具有较少注释基因的那个术语；所述的低表达本体术语是指多样性表达分数在[0,0.1]范围内的术语；

步骤三，计算基因本体术语相似性：使用Schlicker方法计算经步骤二筛选后的基因本体术语相似性；

步骤四，对基因本体术语进行聚类：设定每个簇包含的最大基因本体个数，采用层次聚类方法对基因本体术语进行聚类，得到基因本体术语簇；所述的最大基因本体个数的设定为[1,4]之间的整数；

步骤五，提取Must-Link约束：重复以下步骤1-3共100次，如果两个细胞始终在同一聚类中，则认为这两个细胞存在Must-Link约束：

步骤1：使用PCA方法降低单细胞测序数据的维度；

步骤2：使用拐点法确定细胞类别数目；

步骤3：使用kmeans方法对细胞进行聚类，得到聚类结果；

步骤六，构建simGOAE模型：共包含五层神经元模型，第一层为输入层，神经元代表经过步骤一处理后的单细胞测序数据集中的基因；第二层为本体术语层，神经元代表经过步骤二和步骤三处理过后的基因本体术语，第二层中代表基因本体术语的神经元与第一层中代表该基因本体术语的注释基因的神经元连接，第二层还包括代表基因本体术语簇的神经元，实现对同一簇内的术语神经元的取均值操作；第三层为中间层，包括100个普通神经元，与第二层全连接；第四层与第二层完全相同，与第三层全连接；第五层与第一层完全相同，与第四层部分连接，连接规则和第二层与第一层的连接规则相同；

设定simGOAE模型的损失函数为：

loss＝L_MSE+αL₂+βL_m (1)

其中，loss表示总损失，L_MSE表示输入层与第五层输出层之间的均方误差项，L₂表示正则化项，α表示正则化项权重系数，取值范围为(0,1)，L_m表示Must-Link约束项，β表示Must-Link约束项的权重系数，取值范围为(0,1)；

所述的Must-Link约束项的表达式如下：

其中，x_i和x_j分别为第i和第j个细胞的向量表示，M表示存在Must-Link约束的零一对称矩阵，M_i,j＝1表示第i个细胞和第j个细胞具有Must-Link约束，M_i,j＝0表示第i个细胞和第j个细胞不存在Must-Link约束，n表示细胞个数；h(x_i)表示输入层特征向量x_i在第三层中间层的向量表示；h(x_j)表示输入层特征向量x_j的在第三层中间层的向量表示；

步骤七，对单细胞测序数据进行降维：利用经步骤一处理后的单细胞测序数据对simGOAE模型进行训练，训练后模型第三层的输出结果即为该单细胞测序数据集最终的降维结果。

本发明的有益效果是：由于采用基因本体术语相似性作为生物先验知识加入模型，使得模型部分连接，能够减少训练参数，并使模型具有可解释性，可以更好地挖掘细胞的生物信息；由于提取细胞Must-Link约束作为细胞层面的先验知识加入模型损失函数，可以加速模型训练，增强模型降维效果；由于使用自编码器模型，能够适应大量样本的单细胞测序数据集，实现快速有效地降维。

附图说明

图1是本发明的融合基因本体和神经网络的单细胞测序数据降维方法流程图。

图2是本发明的基因本体相似性提取过程示意图。

图3是本发明的Must-Link约束提取过程示意图。

图4是本发明的simGOAE模型结构图。

具体实施方式

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。

单细胞测序数据可以看做一个横纵坐标分别为细胞和基因的矩阵，矩阵内的数字代表某个细胞在该基因表达量，一般用实数表示。如图1所示，本发明提供了一种融合基因本体和神经网络的单细胞测序数据降维方法，其基本实现过程如下：

1、数据预处理

一般的原始单细胞测序数据中都为自然数，对其进行预处理。

(1)删除单细胞测序数据中基因表达小于3个细胞的基因(基因表达值为 0,1,2···其中细胞的某个基因表达值为0代表该基因在该细胞中未表达)；

(2)删除单细胞测序数据中在基因本体术语注释数据库中没有表达的基因；

(3)对单细胞测序数据在基因上进行z-score标准化。

2、筛选基因本体术语

基因本体(GO)是一种流行的词汇系统，用于系统地描述基因和基因产物的属性。每个GO术语都可以注释一组基因。GO被构造为有向无环图。首先，选取第三层的所有术语。其次，根据特异性分数U_ij删除冗余术语。由于术语本身的注释基因也存在相交的情况，若两个本体术语的注释基因极为相似，认为这两个基因本体术语在生物功能上可能相似，存在冗余。特异性分数U_ij定义如下：

其中，GO_i:{gene₁,gene₂,…,gene_n}代表第i个基因本体术语GO_i包含n个注释基因gene₁,gene₂,…,gene_n，i＝1,2,…，i的上限为第三层GO本体术语的个数。如果两个本体术语的特异性分数U_ij在[0.5,1)范围内，认为这两个基因本体术语不是独特的，删去两个术语中具有较少注释基因的那个术语。最后，根据多样性表达分数H_i删除低表达术语。第i个基因本体术语GO_i多样性表达分数H_i定义如下：

其中，std_j为GO_i的第j个注释基因gene_j的表达标准差。如果一个本体术语的多样性表达分数在[0,0.1]范围内，则认为该术语是低表达术语，删除该术语。

3、计算基因本体术语相似性

一般的，信息量(Information Content,IC)常用来表示某个节点所包含的信息多少。一个术语的信息量定义如下：

其中，IC_term代表术语P的信息量，G代表GO中根节点的基因总数，G_P代表 GO术语P的注释基因个数。

最低共同祖先(Lowest common ancestor,LCA)是一种非常流行的GO术语相似性度量。给定两个GO术语t_a和t_b，定义G_LCA表示这两个GO本体术语的最低共同祖先术语的注释基因个数。GO术语t_a和t_b的相似性计算公式如下：

其中，IC(LCA)表示t_a和t_b的最低共同祖先术语的信息量。

Schlicker等人提出如下术语相似性计算方法：

其中，sim_schlicker(t_a,t_b)表示t_a和t_b的schlicker相似性，IC(LCA)表示t_a和t_b的最低共同祖先术语的信息量，IC(t_a)和IC(t_b)表示t_a和t_b术语的信息量，|G_LCA|表示t_a和t_b的最低共同祖先术语的注释基因个数，|G|表示根节点本体术语的注释基因个数。公式(7) 中前一部分表示标准化术语t_a和t_b信息量和其最低公共祖先的信息量，第二部分是由其最低共同祖先术语信息量在GO中的级别决定的加权分数。

使用如上Schlicker方法计算经步骤2筛选后的基因本体术语相似性。

4、对基因本体术语进行聚类

在使用Schlicker方法对基因本体术语相似性进行计算，得到术语两两之间的相似性之后，本发明使用层次聚类方法对基因本体术语进行聚类。在这里，需要控制基因本体术语簇内的术语个数，因为如果一个簇内存在过多的术语会掩盖此类生物信息的表达，并掩盖每个术语的异质性。层次聚类顺序按照本体术语相似性从高到低的顺序进行聚类，并且，只有当两个基因本体术语的Schlicker相似性大于0.5才有可能被聚在一起，否则认为不存在相似性。本发明设定最大基因本体个数为[1,4]个术语。图2 给出了如上进行基因本体相似性计算和聚类的过程。

5、提取Must-Link约束

若两个细胞在高维空间距离很近，以至于总是被聚成同一类，那么这两个细胞在低维空间的距离也应该很近，这样的两个细胞即为成对约束，即一种细胞与细胞之间的强关系，本发明称之为Must-Link约束，即这样的细胞必须具有高度相似性或属于同一类别。

针对上述定义，对于给定的细胞原始数据，如图3所示，本发明使用下述策略提取约束：

第一，使用PCA方法降低单细胞测序数据的维度。第二，由于细胞类别数目在此刻未知，使用拐点法确定细胞类别数目。第三，使用kmeans方法对细胞进行聚类，并获得预先聚类的结果。重复前述三个步骤100次，如果两个细胞始终在同一聚类中，则认为这两个细胞存在Must-Link约束。

从预先聚类中提取约束后，可以在模型中应用Must-Link约束作为先验信息。对于同一组原始单细胞测序数据，选择的Must-Link约束是相同的。原始的自编码器模型无法对从输入层获得的约束进行建模。因此，本发明在simGOAE模型中考虑结合细胞水平上的先验约束——Must-Link约束。其数学表示如下：

其中，M是表示Must-Link约束的矩阵，该矩阵对称的零一矩阵，M_i,j＝1表示第 i个细胞x_i和第j个细胞x_j具有Must-Link约束，M_i,j＝0表示第i个细胞x_i和第j个细胞x_j不存在Must-Link约束，n表示细胞个数；h(x_i)表示网络输入层特征向量x_i的潜在层表示，h(x_j)表示网络输入层特征向量x_j的潜在层表示。其中，约束对的个数必须远小于总细胞对数：|M|＜＜|S|。假设如果在输入层中细胞x_i和x_j中存在Must-Link约束，则它们在低维空间也应距离很近或更相似。理想情况下，在编码过程后，两个Must-Link 细胞应该比输入阶段距离更相近。令d(h(x_i),h(x_j))表示在编码后(自编码器中间层) 细胞x_i和x_j的之间的距离。对于Must-Link约束，两细胞在编码前的距离d(x_i,x_j)应大于编码后的距离d(h(x_i),h(x_j))。

为了增强模型降维效果，本发明在后续构建的simGOAE模型的损失函数中增加如上Must-Link约束惩罚项。

6、构建simGOAE模型

将步骤2和步骤4得到的本体术语和本体术语簇建模为神经元节点放置在模型第二层和第四层中，其中本体术语神经元与其注释基因部分连接，本发明构建了用于对单细胞测序数据进行降维处理的simGOAE模型，如图4所示。该模型共包含五层神经元。第一层为输入层，神经元代表经过步骤1处理后的单细胞测序数据集中的基因；第二层为本体术语层，神经元代表经过步骤2和步骤3处理过后的基因本体术语，第二层神经元与输入层神经元部分连接，连接规则即为基因本体术语代表的神经元与该基因本体术语的注释基因所代表的神经元连接；第三层为中间层，该层神经元为神经网络普通的神经元，设共有100个神经元，与第二层全连接；第四层与第二层表示含义相同，为本体术语层，与第三层全连接，神经元代表与第二层相同的基因本体术语；第五层与第一层表示含义相同，为输出层，神经元代表与第一层相同的基因，与第四层部分连接，连接规则与第一层与第二层规则相同。

前三层为编码器，第三层到第五层为解码器，与编码器结构对称。

此外，在第二层中还有一类神经元，代表的是基因本体术语簇节点。根步骤3所介绍的方法，获取本体术语相似性特征，并根据术语之间的相似性，将这些术语聚成不超过4个本体术语的簇。对于这些本体术语簇神经元，它的输入即为术语神经元的输出，簇神经元的作用只是对这些相似的术语神经元输出取均值，并无其他激活函数等其他操作。同时，还可以在第二层中，增加一部分稠密神经元(在本模型中设为100 个)与输入层神经元全连接，以保证最大可能获取一些未前述本体术语神经元所捕获的生物信息，增强模型聚类效果。

设定simGOAE模型的损失函数为：

loss＝L_MSE+αL₂+βL_m (9)

其中，loss表示总损失，L_MSE表示输入层与第五层输出层之间的均方误差项，L₂表示正则化项，α表示正则化项权重系数，取值范围为(0,1)，L_m表示Must-Link约束项，β表示Must-Link约束项的权重系数，取值范围为(0,1)。

为了优化公式(9)的目标函数，计算

偏导是至关重要。为了计算该偏导结果，L_m可以用下式表达：

其中，L_M＝D_M-M，D_M∈R^n×n是一个对称矩阵，

H是第三层中间层的简化表示。因此，

可以进行如下化简：

其中，W代表神经网络里的权重矩阵，B代表神经网络里的偏置矩阵。

7、对单细胞测序数据进行降维

使用simGOAE模型训练单细胞测序数据，将训练后第三层的结果作为降维结果输出。

为了证明本发明方法的有效性，分别采用主成分分析(PCA)、独立成分分析(ICA)、非负矩阵分解(NMF)、SIMLR、降噪自编码器(DAE)、深度计数自编码器模型(DAC)、单细胞数据变分自编码器模型(scVI)和本发明的simGOAE方法对单细胞测序数据进行降维处理，再使用kmeans++方法对降维结果进行聚类，使用两个常用的聚类评价指标ARI和NMI，对不同方法的效果定量计算和评价。ARI测量了真实标签和聚类标签之间的相似性，NMI是互信息指数的标准化值，ARI和NMI值越大意味着聚类结果与真实结果越吻合。实验所采用的单细胞测序数据集包括不同平台上的测序数据，分别为Zeisel数据集、Pollen数据集、Baron Human数据集和Zheng数据集。各数据集的信息如表1所示。

表1

数据集名称	平台	物种	样本数量	基因数量	类别数目
						Zeise	STRT-Seq	小鼠	3005	11404	7
Pollen	SMARTer	人类	300	8686	11
						Baron	inDrop	人类	1724	8550	13
Zheng	10X	人类	54621	5831	7

采用不同方法计算得到ARI结果如表2所示，NMI结果如表3所示。可以看出，本发明方法在所有数据集上的结果都要优于其他方法，证明了本发明方法的有效性。

表2

方法	Zeisel	Pollen	Baron	Zheng
					PCA	0.631	0.817	0.699	0.886
ICA	0.052	0.161	0.116	0.467
					NMF	0.041	0.028	0.353	0.167
SIMLR	0.688	0.606	0.491	0.522
					DAE	0.689	0.842	0.781	0.782
DCA	0.625	0.671	0.336	0.816
					scVI	0.355	0.861	0.421	0.874
simGOAE	0.894	0.909	0.826	0.925

表3

方法	Zeisel	Pollen	Baron	Zheng
					PCA	0.736	0. 917	0.808	0.903
ICA	0.265	0.390	0.396	0.633
					NMF	0.176	0.247	0.558	0.386
SIMLR	0.747	0.810	0.695	0.684
					DAE	0.737	0.926	0.826	0.861
DCA	0.704	0.797	0.597	0.856
					scVI	0.696	0.912	0.708	0.883
simGOAE	0.857	0.933	0.828	0.922

由于融合了基因本体术语，simGOAE模型可以通过对神经网络模型的权重分析得到不同细胞的生物信息，使得模型具有一定的可解释性。首先，通过中间层的降维结果，提取单个细胞的低维表示，然后，对于每种细胞类型，找到该细胞类型低维表示中具有最小标准偏差的10个维度对应的中间层节点，这是因为最小标准偏差的意味着这种细胞在该维度上表达最相似，而正是因为相似表达使这些细胞被聚成一类。通过寻找每个具有较小标准偏差前K个权重较大的GO本体术语簇节点(K＜＝10)，从而找到相应类细胞中的重要GO术语簇节点和GO术语节点。

通过对本发明的simGOAE模型的权重分析，还能得到不同细胞的生物信息，即通过降维结果提取单个细胞的低维表示，对于每种细胞类型，找到具有最小标准偏差的K个第三层中间层节点，即GO簇节点，从而找到相应类细胞中的重要GO术语。表4给出了利用simGOAE模型在Zeisel数据集上找出的重要GO术语。例如，GO： 0006810是星形胶质细胞-表皮(astrocytes-ependymal)细胞类型的重要术语，它的主要GO功能是借助某种物质(例如转运蛋白，孔蛋白或运动蛋白)形成的多细胞生物。另一个术语GO：0030030是在细胞水平上执行的过程，该过程导致装配，组件对齐或因细胞伸长或诸如鞭毛或轴突的过程分解而伸长。因此，由于本发明的simGOAE模型融合了基因本体信息，具有可解释性，可以用于进行细胞相关生物信息挖掘。

表4

Claims (1)

1.一种融合基因本体和神经网络的单细胞测序数据降维方法，其特征在于步骤如下：

步骤一，数据预处理：删除单细胞测序数据中基因表达小于3个细胞的基因；删除单细胞测序数据中在基因本体术语注释数据库中没有表达的基因；对单细胞测序数据在基因上进行z-score标准化；

步骤二，筛选基因本体术语：删除冗余本体术语和低表达本体术语；所述的冗余本体术语是指特异性分数在[0.5,1)范围内的两个基因本体术语中具有较少注释基因的那个术语；所述的低表达本体术语是指多样性表达分数在[0,0.1]范围内的术语；

步骤三，计算基因本体术语相似性：使用Schlicker方法计算经步骤二筛选后的基因本体术语相似性；

步骤四，对基因本体术语进行聚类：设定每个簇包含的最大基因本体个数，采用层次聚类方法对基因本体术语进行聚类，得到基因本体术语簇；所述的最大基因本体个数的设定为[1,4]之间的整数；

步骤五，提取Must-Link约束：重复以下步骤1-3共100次，如果两个细胞始终在同一聚类中，则认为这两个细胞存在Must-Link约束：

步骤1：使用PCA方法降低单细胞测序数据的维度；

步骤2：使用拐点法确定细胞类别数目；

步骤3：使用kmeans方法对细胞进行聚类，得到聚类结果；

步骤六，构建simGOAE模型：共包含五层神经元模型，第一层为输入层，神经元代表经过步骤一处理后的单细胞测序数据集中的基因；第二层为本体术语层，神经元代表经过步骤二和步骤三处理过后的基因本体术语，第二层中代表基因本体术语的神经元与第一层中代表该基因本体术语的注释基因的神经元连接，第二层还包括代表基因本体术语簇的神经元，实现对同一簇内的术语神经元的取均值操作；第三层为中间层，包括100个普通神经元，与第二层全连接；第四层与第二层完全相同，与第三层全连接；第五层与第一层完全相同，与第四层部分连接，连接规则和第二层与第一层的连接规则相同；

设定simGOAE模型的损失函数为：

loss＝L_MSE+αL₂+βL_m (1)

其中，loss表示总损失，L_MSE表示输入层与第五层输出层之间的均方误差项，L₂表示正则化项，α表示正则化项权重系数，取值范围为(0,1)，L_m表示Must-Link约束项，β表示Must-Link约束项的权重系数，取值范围为(0,1)；

所述的Must-Link约束项的表达式如下：

其中，x_i和x_j分别为第i和第j个细胞的向量表示，M表示存在Must-Link约束的零一对称矩阵，M_i,j＝1表示第i个细胞和第j个细胞具有Must-Link约束，M_i,j＝0表示第i个细胞和第j个细胞不存在Must-Link约束，n表示细胞个数；h(x_i)表示输入层特征向量x_i在第三层中间层的向量表示；h(x_j)表示输入层特征向量x_j的在第三层中间层的向量表示；

步骤七，对单细胞测序数据进行降维：利用经步骤一处理后的单细胞测序数据对simGOAE模型进行训练，训练后模型第三层的输出结果即为该单细胞测序数据集最终的降维结果。

Images