基于混合实验组和对照组单细胞样本的差异分析方法及系统

技术领域

本发明涉及一种基于混合实验组和对照组单细胞样本的差异分析方法及系统。

背景技术

高通量的单细胞技术目前是生物信息领域的热门研究方向之一，它带领了基因组学领域的转变，即从整体组织分析转向对单个细胞的详细和全面的研究，提供了对不同生物体组织中细胞多样性的洞察力。当前许多单细胞研究的目的主要集中在鉴定细胞类型(如B细胞，T细胞，以及表达独特基因的细胞亚群)，很少聚焦于在单细胞层面分析实验组(即疾病组)和对照组的单细胞表达差异，且直接将全部单细胞按照实验组和对照组进行整体差异分析捕捉到的差异信号较弱，效果不佳。

在癌症等疾病中，由于组织中某一群细胞的病变就可能导致整个有机体的崩溃，而基于组织样本的分析仅仅获取了组织中所有细胞的平均表达信息，很难分辨细胞间的差异，也很难识别在疾病进展中起到重要作用的稀有细胞群。单细胞技术的出现有效地解决了这一问题，单细胞技术是在组织中分离单细胞并对其进行测序以获得单个细胞的基因表达的技术，它的发展使基因表达的研究达到了更高的分辨率。单细胞技术分为上游和下游两个层面，上游层面指的是细胞样本的制备和测序技术，包括单细胞分离、核酸提取和扩增、测序库制备、RNA测序。而下游层面则是针对单细胞表达矩阵的数据分析技术，包括数据质量控制、数据降维、聚类、差异分析等一系列数据分析。

目前单细胞研究的主要方向是鉴定新的细胞亚群类型和对疾病发展有重要作用的特殊细胞亚群。单细胞研究的取样策略通常有两种：一是通过实验获取疾病组织样本(如肺癌，肺纤维化的肺部组织样本)的单细胞，只针对疾病单细胞样本进行分析；二是选取疾病组织和健康组织的单细胞样本共同分析。两种方式都是通过从组织中分离单细胞进行测序，并鉴定细胞亚群。第一种取样方式更关注于鉴定疾病进程中的关键细胞亚群，即关键基因表达上调的细胞亚群；而第二种方法更关注细胞类型的鉴定，即分析实验组和对照组单细胞样本鉴定组织中不同的细胞类型(如B细胞、T细胞、巨噬细胞及其亚型等)。

上述提到的单细胞分析的主要任务是鉴定关键细胞亚群，所述关键细胞亚群可能表达某些独特基因并对疾病的发展产生重要的影响。但有一个问题一直是单细胞研究中很少关注到的，即是否能够从疾病样本和对照样本单细胞中找到有着共同表达模式(即共同表达某些基因)，且单独由疾病样本或对照样本细胞构成的细胞亚群。

寻找上述细胞亚群的意义在于可以更精细化地进行差异分析。在之前的分析中，发现许多实验组和对照组的细胞亚群共享某种基因表达模式，即在多个细胞亚群中实验组和对照组的单细胞数目相近，这也意味着这些细胞不能体现实验组和对照组的表达差异。而细胞占比差距悬殊，尤其是一个基本由疾病样本细胞构成的细胞亚群较为罕见，而往往这类罕见的细胞群体高度表达了疾病相关的重要基因。如果能够在混合两组样本的细胞中找到这样的细胞亚群，就可以为解释疾病的发病机制提供独特的见解。

在组织层面的研究中，研究人员会直接使用实验组和对照组的组织样本的基因表达数据进行差异分析，但这种分析方法在单细胞层面是行不通的。单细胞测序数据的特点是高维稀疏，仅一个组织样本就能分离出几千个单细胞，并且一个单细胞只表达一部分基因。在这样的数据背景下直接将实验组和对照组的细胞样本整体进行对比会使得到的差异信号很弱。

发明内容

有鉴于此，有必要提供一种基于混合实验组和对照组单细胞样本的差异分析方法及系统。

本发明提供一种基于混合实验组和对照组单细胞样本的差异分析方法，该方法包括如下步骤：a.获取数据，并对所述数据进行质量控制及标准化处理，所述数据包括实验数据和开源数据；b.从标准化处理后的数据中筛选出高变基因的数据，并对筛选出的数据进行降维；c.根据降维后的数据，对单细胞进行多轮不同参数的聚类分析，计算每轮聚类分析的细胞亚群实验组和对照组的构成比例，以得到最优结果；d.对细胞亚群进行差异分析。

优选地，所述的步骤a包括：

对原始表达矩阵进行质量控制，包括剔除低质量的细胞和基因；其中，低质量的细胞是指：表达基因数量较少，线粒体基因表达比例较高的细胞；低质量的基因是指：在细胞中表达比例较低的基因。

优选地，所述的步骤a还包括：

通过对唯一分子识别标签计数归一化缩小数值差异，采用LogNormalization方法，公式如下：

其中，i表示细胞，j表示基因；

在对数据归一化之后，继续对数据进行缩放，公式如下：

优选地，所述的步骤b包括：

通过对基因在细胞中的表达均值和方差进行建模，筛选高变基因，对单细胞在高变基因中的表达数据进行主成分分析，并可视化每个主成分的方差百分比，选择排名靠前的主成分用于后续聚类分析。

优选地，所述的步骤c具体包括：

采用基于图的聚类方法：首先计算每个细胞之间的相互距离以确定每个细胞的k最近邻，然后计算每个细胞与其k最近邻之间的邻域重叠构造共享最近邻图；公式如下，其中A，B代表两个细胞的“邻居”细胞数：

优选地，所述的步骤c还包括：

使用Louvain算法作为模块化优化函数确定细胞亚群的数量：通过扫描节点及其邻居节点，对每个节点计算模块度并衡量其加入模块后的模块收益，选择最大收益的邻居节点加入模块，重复迭代，最终形成节点簇；对于加权图，模块性计算公式如下，其中，m是图中边总数，i，j代表两个节点(即细胞)，A

优选地，所述的步骤d具体包括：

对每群细胞亚群计算差异基因：比较本群细胞和除本群细胞之外所有细胞的基因表达差异，使用Wilcoxon秩和检验方法进行差异表达分析。

本发明提供一种基于混合实验组和对照组单细胞样本的差异分析系统，该系统包括获取模块、筛选模块、聚类分析模块、差异分析模块，其中：所述获取模块用于获取数据，并对所述数据进行质量控制及标准化处理，所述数据包括实验数据和开源数据；所述筛选模块用于从标准化处理后的数据中筛选出高变基因的数据，并对筛选出的数据进行降维；所述聚类分析模块用于根据降维后的数据，对单细胞进行多轮不同参数的聚类分析，计算每轮聚类分析的细胞亚群实验组和对照组的构成比例，以得到最优结果；所述差异分析模块用于对细胞亚群进行差异分析。

本发明改变了现有的单细胞数据分析任务多数聚焦于鉴定表达关键基因的细胞亚群或新的细胞类型的传统做法，在鉴定细胞亚群的基础上，更加关注实验组和对照组单细胞样本的差异研究，即通过本申请能够更好地发现由特定组构成的细胞亚群，从而实现在单细胞层面的实验组和对照组的差异分析。

附图说明

图1为本发明基于混合实验组和对照组单细胞样本的差异分析方法的流程图；

图2为本发明实施例基于混合实验组和对照组单细胞样本的差异的流程示意图；

图3为本发明实施例提供的代表细胞亚群走向的聚类树(以GSE36831数据集中的巨噬细胞为例)；

图4为本发明实施例提供的细胞亚群构成图(以GSE36831数据集中的巨噬细胞为例)；

图5为本发明实施例提供的最佳聚类结果评估流程图；

图6为本发明基于混合实验组和对照组单细胞样本的差异分析系统的硬件架构图。

具体实施方式

下面结合附图及具体实施例对本发明作进一步详细的说明。

参阅图1、图2所示，是本发明基于混合实验组和对照组单细胞样本的差异分析方法较佳实施例的作业流程图。

步骤S1，获取数据，并对所述数据进行质量控制及标准化处理。其中，所述数据包括实验数据和开源数据。具体而言：

本实施例分析了两个数据集：第一个数据集是小鼠肺纤维化造血细胞(Cd45+)数据集，该数据集由实验获得，包含了两个不同时间点的实验组和对照组，共四组单细胞表达数据，每组检测了约5000+单细胞在30000+基因中的表达量；第二个数据集是从GEO数据库下载的人特发性肺纤维化肺组织单细胞数据集(GSE36831)，包括28个健康个体(对照组)和31个特发性肺纤维化(IPF)患者(实验组)的肺组织单细胞样本，共检测了380000+单细胞在40000+基因中的表达量。

首先，对原始表达矩阵进行质量控制，包括剔除低质量的细胞(表达基因数量较少，线粒体基因表达比例较高)和基因(基因在细胞中表达比例较低)。由于每个细胞在不同基因的表达读数差异非常大，影响后续数据分析，因此通过对UMIs计数归一化缩小数值差异。所述表达读数即是唯一分子识别标签(unique molecular identifier，UMI)。这里本实施例采用LogNormalization方法，公式如下，其中i表示细胞，j表示基因：

在对数据归一化之后，继续对数据进行缩放(即线性变换)，这是在执行降维和聚类之前的标准预处理步骤，公式如下：

步骤S2，从标准化处理后的数据中，筛选高度变化特征(Highly VariableFeatures)，即在不同细胞中表达差异较大的基因，并对数据降维。具体而言：

由于单细胞表达矩阵是稀疏矩阵，多数基因在细胞中的表达差异较小，且将所有基因用于分析会占用大量计算资源。因此，选择高度变化的特征(即高变基因)进行后续的降维和聚类分析会节省很多计算资源，并且有助于突出单细胞数据集中的生物信号。通过对基因在细胞中的表达均值和方差进行建模，筛选高变基因，对单细胞在高变基因中的表达数据进行主成分分析(PCA)，并可视化每个主成分(PC)的方差百分比(SD)，选择排名靠前的主成分用于后续聚类分析。

步骤S3，根据降维后的数据，对单细胞进行多轮不同参数的聚类分析，计算每轮聚类分析的细胞亚群实验组和对照组的构成比例，以得到最优结果：

聚类分析用于两个任务：第一个是细胞分型，即从组织分离的细胞样本中标记每个细胞的细胞类型(如巨噬细胞，B细胞，T细胞等)；第二个是对实验组和对照组的差异分析。这一步在已经确定每个细胞的细胞类型之后进行，通过提取同一细胞类型的实验组和对照组的单细胞样本，调整参数进行多次聚类，评估每次聚类结果，选择最佳聚类结果。

具体而言：

在本步骤中，本实施例采用基于图的聚类方法。该方法首先计算每个细胞之间的相互距离以确定每个细胞的k最近邻(k-nearest neighbors，k默认为20)，然后计算每个细胞与其k最近邻之间的邻域重叠(Jaccard指数)构造共享最近邻图(Shared NearestNeighbor，SNN)。Jaccard指数计算的是两个细胞共有的“邻居”和所有“邻居”的比值，比值越大，两个细胞越相似，公式如下，其中A，B代表两个细胞的“邻居”细胞数：

然后，使用Louvain算法作为模块化优化函数确定簇(细胞亚群)的数量。Louvain算法是一种从网络中提取模块的方法，是一种贪婪优化方法。Louvain算法通过扫描节点及其邻居节点，对每个节点计算模块度并衡量其加入模块后的模块收益，选择最大收益的邻居节点加入模块，重复迭代，最终形成节点簇。对于加权图，模块性计算公式如下，其中，m是图中边总数，i，j代表两个节点(即细胞)，A

本实施例利用聚类分析完成两个任务：

第一个任务是细胞分型，即对组织样本分离的单细胞鉴定细胞类型。这一步通过对所有细胞进行聚类分析，划分细胞亚群，计算每个亚群表达上调的基因。通过在线资源(文献、数据库)人工分析所述上调的基因是否包含某一细胞类型的markers，从而判断细胞亚群所属细胞类型。

第二个任务是本实施例的关键部分，即对混合实验组和对照组特定细胞类型的单细胞样本聚类。首先，选取上一步分析的特定细胞类型的两组细胞(如巨噬细胞)，调整参数进行多次聚类，构建聚类树，可视化每次调整参数后细胞分群走向(如图3)，并对聚类后每个细胞亚群计算实验组和对照组构成比例(如图4)，选择最佳聚类结果(最佳聚类结果评估流程如图5所示)。

所述最佳聚类结果需综合考虑两组样本构成比例差距悬殊的细胞亚群数目以及所述细胞亚群的细胞数量。如果两者都是最多，则该参数下的聚类结果为最佳聚类结果；否则，如果比例差距悬殊的细胞亚群数目不是最多，但所述细胞亚群的细胞数足够多，则优先考虑这样的聚类结果。细胞数量太少的亚群不能有力的证明其对疾病的影响。

步骤S4，在所述聚类分析完成后，对细胞亚群进行差异分析。具体而言：

对每群细胞亚群计算差异基因。本实施例的计算策略是比较本群细胞和除本群细胞之外所有细胞的基因表达差异，使用了Wilcoxon秩和检验(Wilcoxon Rank Sum test)方法进行差异表达分析。本实施例限制两个细胞群中基因的平均|log fold change|>0.25，并且限制被比较的两个细胞群中任意一群细胞的各基因表达占比>0.1，且差异基因的FDR<0.05。

参阅图6所示，是本发明基于混合实验组和对照组单细胞样本的差异分析系统10的硬件架构图。该系统包括：获取模块101、筛选模块102、聚类分析模块103、差异分析模块104。其中：

所述获取模块101用于获取数据，并对所述数据进行质量控制及标准化处理。其中，所述数据包括实验数据和开源数据。具体而言：

本实施例分析了两个数据集：第一个数据集是小鼠肺纤维化造血细胞(Cd45+)数据集，该数据集由实验获得，包含了两个不同时间点的实验组和对照组，共四组单细胞表达数据，每组检测了约5000+单细胞在30000+基因中的表达量；第二个数据集是从GEO数据库下载的人特发性肺纤维化肺组织单细胞数据集(GSE36831)，包括28个健康个体(对照组)和31个特发性肺纤维化(IPF)患者(实验组)的肺组织单细胞样本，共检测了380000+单细胞在40000+基因中的表达量。

首先，所述获取模块101对原始表达矩阵进行质量控制，包括剔除低质量的细胞(表达基因数量较少，线粒体基因表达比例较高)和基因(基因在细胞中表达比例较低)。由于每个细胞在不同基因的表达读数差异非常大，影响后续数据分析，因此通过对UMIs计数归一化缩小数值差异。所述表达读数即是唯一分子识别标签(unique molecularidentifier，UMI)。这里本实施例采用LogNormalization方法，公式如下，其中i表示细胞，j表示基因：

在对数据归一化之后，所述获取模块101继续对数据进行缩放(即线性变换)，这是在执行降维和聚类之前的标准预处理步骤，公式如下：

所述筛选模块102用于从标准化处理后的数据中，筛选高度变化特征(HighlyVariable Features)，即在不同细胞中表达差异较大的基因，并对数据降维。具体而言：

由于单细胞表达矩阵是稀疏矩阵，多数基因在细胞中的表达差异较小，且将所有基因用于分析会占用大量计算资源。因此，选择高度变化的特征(即高变基因)进行后续的降维和聚类分析会节省很多计算资源，并且有助于突出单细胞数据集中的生物信号。所述筛选模块102通过对基因在细胞中的表达均值和方差进行建模，筛选高变基因，对单细胞在高变基因中的表达数据进行主成分分析(PCA)，并可视化每个主成分(PC)的方差百分比(SD)，选择排名靠前的主成分用于后续聚类分析。

所述聚类分析模块103用于根据降维后的数据，对单细胞进行多轮不同参数的聚类分析，计算每轮聚类分析的细胞亚群实验组和对照组的构成比例，以得到最优结果：

聚类分析用于两个任务：第一个是细胞分型，即从组织分离的细胞样本中标记每个细胞的细胞类型(如巨噬细胞，B细胞，T细胞等)；第二个是对实验组和对照组的差异分析。这一步在已经确定每个细胞的细胞类型之后进行，通过提取同一细胞类型的实验组和对照组的单细胞样本，调整参数进行多次聚类，评估每次聚类结果，选择最佳聚类结果。

具体而言：

所述聚类分析模块103采用基于图的聚类方法。首先计算每个细胞之间的相互距离以确定每个细胞的k最近邻(k-nearest neighbors，k默认为20)，然后计算每个细胞与其k最近邻之间的邻域重叠(Jaccard指数)构造共享最近邻图(Shared Nearest Neighbor，SNN)。Jaccard指数计算的是两个细胞共有的“邻居”和所有“邻居”的比值，比值越大，两个细胞越相似，公式如下，其中A，B代表两个细胞的“邻居”细胞数：

然后，使用Louvain算法作为模块化优化函数确定簇(细胞亚群)的数量。Louvain算法是一种从网络中提取模块的方法，是一种贪婪优化方法。Louvain算法通过扫描节点及其邻居节点，对每个节点计算模块度并衡量其加入模块后的模块收益，选择最大收益的邻居节点加入模块，重复迭代，最终形成节点簇。对于加权图，模块性计算公式如下，其中，m是图中边总数，i，j代表两个节点(即细胞)，A

本实施例利用聚类分析完成两个任务：

第一个任务是细胞分型，即对组织样本分离的单细胞鉴定细胞类型。这一步通过对所有细胞进行聚类分析，划分细胞亚群，计算每个亚群表达上调的基因。通过在线资源(文献、数据库)人工分析所述上调的基因是否包含某一细胞类型的markers，从而判断细胞亚群所属细胞类型。

第二个任务是本实施例的关键部分，即对混合实验组和对照组特定细胞类型的单细胞样本聚类。首先，选取上一步分析的特定细胞类型的两组细胞(如巨噬细胞)，调整参数进行多次聚类，构建聚类树，可视化每次调整参数后细胞分群走向(如图3)，并对聚类后每个细胞亚群计算实验组和对照组构成比例(如图4)，选择最佳聚类结果(最佳聚类结果评估流程如图5所示)。

所述最佳聚类结果需综合考虑两组样本构成比例差距悬殊的细胞亚群数目以及所述细胞亚群的细胞数量。如果两者都是最多，则该参数下的聚类结果为最佳聚类结果；否则，如果比例差距悬殊的细胞亚群数目不是最多，但所述细胞亚群的细胞数足够多，则优先考虑这样的聚类结果。细胞数量太少的亚群不能有力的证明其对疾病的影响。

所述差异分析模块104用于在所述聚类分析完成后，对细胞亚群进行差异分析。具体而言：

所述差异分析模块104对每群细胞亚群计算差异基因。本实施例的计算策略是比较本群细胞和除本群细胞之外所有细胞的基因表达差异，使用了Wilcoxon秩和检验(Wilcoxon Rank Sum test)方法进行差异表达分析。本实施例限制两个细胞群中基因的平均|log fold change|>0.25，并且限制被比较的两个细胞群中任意一群细胞的各基因表达占比>0.1，且差异基因的FDR<0.05。

本申请通过混合特定细胞类型的实验组和对照组单细胞样本进行多轮聚类分析，评估每轮聚类结果，选择最佳聚类结果并进行差异研究。与基于组织测序数据的差异分析方法不同，本申请首先通过聚类分析找到实验组和对照组细胞占比悬殊的细胞亚群，然后分析上述细胞亚群的表达模式，比直接分析实验组和对照组的单细胞具备更强的差异基因信号。

虽然本发明参照当前的较佳实施方式进行了描述，但本领域的技术人员应能理解，上述较佳实施方式仅用来说明本发明，并非用来限定本发明的保护范围，任何在本发明的精神和原则范围之内，所做的任何修饰、等效替换、改进等，均应包含在本发明的权利保护范围之内。