本发明涉及生物信息学中的数据挖掘领域，具体涉及一种基于流形学习和主曲线的轨迹推断方法。
背景技术：
重建组织或有机体内细胞间的谱系关系是生物学的一个长期目标，了解组织和有机体形成的谱系是生物学的基本问题之一。确定这些关系不仅可以提供有关正常组织发育和体内平衡的宝贵信息，而且可以提供有关发育障碍和疾病(如癌症)的宝贵信息。历史上，谱系追踪是通过在细胞中引入可遗传标记，然后追踪其后代来完成的。组成后代的不同细胞类型在发育上是相关的，因为所有这些标记细胞都来自同一个生成细胞。此外，在后代中发现的细胞类型的多样性体现了生成细胞的潜力。为了准确预测生成细胞的潜能，谱系追踪需要精确的细胞类型识别。理想情况下，人们会使用尽可能多的标记来实现更加精确的细胞类型分类。但是细胞类型识别通常基于有限数量的标记，因此潜在地掩盖了表达所选择的标记基因在细胞亚群内的变异性。因此，这种方法可能会对器官的复杂性产生偏见。
单细胞测序技术的迅速发展使我们能够以前所未有的分辨率探索生物系统。现在可以很容易地描述单个细胞而不是细胞群，这促进了我们对细胞内在异质性和动力学的基本理解。单细胞测序方案已经被开发用来测量不同的分子层，包括转录组学、表观基因组学和蛋白质组学。这些强大的测量手段的结合使得在多组学尺度上研究基因调控等重要的生物过程成为可能。单细胞组学数据，包括转录组学、蛋白质组学和表观基因组学数据，为研究细胞周期、细胞分化和细胞激活等细胞动力学过程提供了新的机会。这种动态过程可以使用轨迹推断(ti)方法(也称为伪时间分析)进行计算建模，这种方法根据细胞表达模式的相似性，可以对细胞沿轨迹进行排序。
尽管有这些技术突破，但由于单细胞测序数据的内在特征，包括细胞间的变异、数据的稀疏性、生物和技术噪音以及退出事件，分析和计算方面仍然存在一些挑战。有鉴于此，本发明提出了一种基于流形学习和主曲线的轨迹推断方法，其能够对单细胞动态分化轨迹过程进行建模和分析。
技术实现要素：
本发明提出了一种基于流形学习和主曲线的轨迹推断方法。用以重建组织或者有机体内细胞间的谱系关系，确定这些关系不仅可以提供有关正常组织发育和体内平衡的宝贵信息，而且可以提供有关疾病(如癌症)的宝贵信息。主要包括以下步骤：
(1)数据集收集阶段，收集已知单细胞rna-seq数据；
(2)进行特征提取，选择可变基因作为特征；
(3)进行数据降维，缓解维度诅咒，降低数据处理难度；
(4)局部定义主曲线，提出了一个初始化过程，提高了推断解的质量，并加快了收敛速度；
(5)分段子空间约束的均值移动算法建立最终主曲线模型；
(6)伪时间分析并建立直树拓扑图；
(7)差异表达基因检测。
1.数据收集阶段
首先，获取单细胞表达基因数据，我们收集了两个rna-seq数据集。一个是已经发布的scrna-seq数据集nestorowa，在这个数据集中，1656个来自小鼠造血系统的单细胞被分类和分析，数据集使用了单细胞rna测序来分析超过1656个单细胞热休克蛋白，而深度测序使每个细胞平均可以检测到6558个蛋白编码基因。索引分类，结合广泛的分类门，使我们能够回顾性地将细胞分配到12种通常分类的hspc表型，同时也捕获传统门控通常排除的中间细胞。我们还使用了最常用的scrna-seq数据集，最初由trapnell等人生成。该数据集包含沿着线性轨迹分化的人类骨骼肌成肌细胞(hsmm)细胞。其使用一种无监督算法monocle，使用在多个时间点收集的单细胞rna-seq数据增加转录组动力学的时间分辨率，能够将单片眼镜应用于原代人成肌细胞的分化，发现关键调控因子表达的开关样变化，基因调控的序列波，以及尚未发现在分化中起作用的调控因子的表达。
2.特征提取阶段
针对于单细胞rna-seq数据，模型输入的基因表达矩阵是将细胞作为矩阵的行，基因作为矩阵的列，矩阵的每一项值则为该基因在该细胞中的表达值，该基因表达值通过库大小归一化和log2转化进行调整。通过计算，绝大多数可变基因被选择为特征。简单的说，对于每个基因，计算每个基因的平均值和标准差，我们利用一种非参数回归方法(loess)去拟合均值和标准差之间的关系。最后，我们选择了曲线上方的差异显著的基因作为可变基因。
3.数据降维阶段
在处理单细胞rna-seq数据的时候，我们必须要面临的一个问题就是维度诅咒。维度诅咒所指的问题是，当维度增加时，空间的体积增加得非常之快(以指数增加)，以致于可用的数据变得稀疏。这种稀疏性会造成统计意义上的一些困难，这是因为为了获得一个统计上可靠的结果，支持结果所需的数据量往往随着维度的增加而成倍增长。即数据中的每个细胞都可以被认为是多维向量空间中的一个向量，在多维向量空间中，每个成分都是一个基因的表达水平。在特征选择之后，每个细胞仍然有数百个分量，这对于可靠的获取细胞之间的相似度或者距离来说非常的困难。为此，通常使用降维来缓解这一问题。然而，尽管目前提出了许多方法来进行降维，但没有一种方法在不同的数据中工作效果是良好的。因此，我们使用了多种方法进行降维，其中包括了mlle、umap、pca和se方法，这些方法被认为对大多数数据集都是有效的。最后，分支的数量和要学习的结构的复杂性决定了组件的数量，通常，高维数据包含更多的信息，这使得结果更加准确。但为了使过程和结果可视化，我们一般选择三个组件，它们可以捕捉大多数数据集的主要结构。所以，降维后我们得到了一个以n个细胞为行、d个分量为列的矩阵。
4.局部定义主曲线
主曲线被hastie算法和stuetzle算法定义为通过数据中间的有条理的一条光滑曲线。然而，hastie算法没有收敛性的证明，使得无法进行理论分析。还应该指出的是，这种主曲线的定义要求主曲线不能与自身相交，这是非常严格的。为了克服这些局限性，ozertem和erdogmus提出了一种新的主曲线定义。他们的算法提出，rd中的d维主曲面就是局部正交d-d维子空间中概率密度函数的局部极大点的集合。对于主曲线，即一维流形，只需代入d＝1。
5.基于分段子空间约束的均值漂移分割算法
我们提出了基于分段子空间约束的均值漂移分割算法(scms)来寻找满足ozertem和erdogmus所给出的主曲线定义的点，从而无需额外的工作就可以处理环路和自交。该算法对mean-shift(ms)均值漂移分割算法进行了改进，通过在轨迹当前点对应(n-1)个特征向量方向的正交空间中约束固定点迭代，从而导致更新收敛于主曲线而不是局部极大值。
与ms算法类似，scms算法假设一个数据的基础kde概率密度。假设xi∈rd，i＝1，...，n.为n个独立同分布的随机变量序列。用高斯核函数g(x)给出了任意点x的kde核密度估计p(x)。p(x)计算公式如下：
当我们使用各向异性变量核函数的一般情况下时，σi是xi的核协方差(对于各向同性核，可以使用标量值而不是全协方差)，其中：
最初，scms算法将细胞轨迹初始化为数据点，并设置j＝0。将各向异性高斯核(或各向同性高斯核的带宽σ)的高斯核协方差矩阵h输入到算法中。接下来，对每个轨迹求均值移位向量m(x(j))，其定义如下：
其中σ-1(x)为任意点x处的局部逆协方差矩阵，其定义如下：
σ-1(x)＝-p-1(x)h(x)+p-2g(x)gt(x)(4)
其中，
然后进行特征分解，其中σ-1(x(j))＝vγv。对于d＝1的特殊情况，即主曲线，我们让v⊥＝[v1…vn-1]为σ-1的(n-1)最大特征向量。在x点上，子空间的均值漂移更新是通过将x投影到受限空间来实现的，即
为了解决这个问题，为了解决这一问题，我们在采用scms算法得到离散主曲线后，构造一棵最小生成树(mst)，然后将不在mst上的所有数据点分配给最近的分支，从而得到一棵包含所有数据点的初始树。在这一步中，分支点被分配给与其相邻的每个分支。接下来，我们将scms算法分别应用于每个分支上的所有数据点。由于scms算法使数据点收敛到主曲线上，分离出初始的树分支，因此我们从任意端点开始深度优先搜索，并基于相同的分支点重新连接分离出的分支来重建最终的树结构。
6.伪时间分析并建立直树拓扑图
我们从任意一个称为原点的端点开始，在d维空间中根据欧氏距离计算每个单元的伪时间。简单地说，我们通过迭代将问题转化为具有固定起点的旅行商问题(tsp)来计算伪时间。
对于一个有n个数据点和m个分支的树形拓扑。对于从基点到其端点或分支点的每个分支bi＝{pi，p2，...，pl}，i＝1，2，...，m，初始化基点pb＝p1，j＝0。
其算法流程如下：
1.找到bi中最前面的k个最近的数据点pb1，pb2，...，pbk，并将它们从bi中移除。
2.将遗传算法(ga)应用于以基点为固定起点和他最近的k个数据点，求解一个以基点为固定起点的tsp问题，得到最短路径rij，然后j+1。
3.重新设置基点pb＝pbk，重复执行步骤2，直到bi为空。
4.将所有最短路径连接起来，得到bi的最短路径ri。
5.在得到每个分支b1，b2，...，bm的最短路径后，我们将它们连接起来，最后得到一个根为原点的最短路径树。然后我们设置原点的伪时间t0＝0，并通过计算每个数据点沿路径到根的距离来分配伪时间。
我们提出了一种以伪时间将细胞排列在平行分枝上的直树图。从原点开始，然后使用广度优先搜索对二维平面上的节点和边进行水平排列，x轴表示伪时间。然后细胞的分化情况在它们的伪时间和它们所属的分支中被映射到所属结构。最后，每个细胞都根据其细胞标签着色。如果要显示基因表达，每个细胞也可以根据其基因表达情况进行着色。
7.差异表达基因检测
我们假设差异表达基因是那些表达值沿着任意直线树的规律变化的基因。在分支众多的复杂树形结构中，变化可能是线性的，也可能是非线性的。常用的相关分析方法如斯皮尔曼秩相关系数、皮尔逊相关系数等只能得出变量之间的线性关系。这里我们用最大信息系数(maximuminformationcoefficient,mic)来衡量两个变量x和y之间的相关程度，不仅可以是线性关系，也可以是非线性关系。mic值越大，变量之间的直接关系越显著。当mic为1时，两者完全相关；当mic为0时，它们完全不相关。
mic的基本原理是基于互信息的，互信息是可模拟的，其定义如下：
其中p(x,y)为x和y之间的联合概率。但一般来说，联合概率是很难得到的。mic基于这样一种理念:如果两个变量之间存在一种关系，那么可以在两个变量的散点图上绘制网格，这两个变量对数据进行了划分，以封装这种关系。从而解决了互信息条件下的联合概率问题。于是，mic可通过以下方式计算:
这里a和b分别为x方向和y方向的划分数。b是一个变量，它的值被设置我们数据量的0.6次方左右。
在该方法中，我们将细胞的伪时间设为x，细胞的基因表达值设为y。计算x与y之间的mic值，以衡量基因表达值是否随伪时间发生显著变化。我们计算每个基因和伪时间之间的mic，并对每个基因序列的mic进行排序。因此，获得了对每条路径中影响最大的基因。对于每个基因，我们还计算了在所有常规程序上的平均mic和最大mic，这可以帮助我们找到在整个细胞中差异表达的基因。
具体实施方式
本发明是一种基于流形学习和主曲线的轨迹推断方法。下面描述本发明的具体实施方式。本领域技术人员应该理解的是，这些实施方式仅仅用于解释本发明的技术原理，并非旨在限制本发明的取证范围。
步骤1：从数据库中下载数据集nestorowa，在这个数据集中，1656个来自小鼠造血系统的单细胞被分类和分析，数据集使用了单细胞rna测序来分析超过1656个单细胞热休克蛋白，最终数据矩阵包含1656个细胞，4768个基因。还使用了最常用的scrna-seq数据集，最初由trapnell等人生成。该数据集包含沿着线性轨迹分化的人类骨骼肌成肌细胞(hsmm)细胞，最终数据矩阵包含271个细胞，47192个基因。
步骤2：根据上述提取出细胞信息进行特征提取。模型输入矩阵的每一项值则为该基因在该细胞中的表达值，该基因的表达值通过库规范化和log2转化进行调整。通过计算，绝大多数可变基因被选择为特征。
步骤3：将得到的特征矩阵进行降维，使用了多种方法进行降维，其中包括了mlle、umap、pca和se方法，最终得到一个包含3列的低秩矩阵，其包含了矩阵中的大多数特征信息。
步骤4：将步骤3得到的特征矩阵初始化，利用ozertem算法和erdogmus算法提出的一种新的主曲线定义，计算出一条局部主曲线，其为通过数据中间的有条理的一条光滑曲线。
步骤5：再将步骤3得到的特征矩阵进行应用分段子空间约束的均值漂移分割算法，来寻找满足ozertem和erdogmus所给出的主曲线定义的点。利用上述公式scms方法建立出相对分散的主曲线(initialtree)，然后构建一个最小生成树结构，其用一个二维数组表示，每一行代表一个树枝元素。再分别对每一个树枝中的数据应用一次上述公式scms方法建立出最终的主曲线模型(finaltree)。
步骤6：利用步骤4中得到的主曲线模型进行拟时间任务，建立直树拓扑结构。将上诉得到的最终主曲线模型利用tsp算法来进行伪时间分析，通过计算每个数据点沿路径到根的距离来分配伪时间。从任意一点开始，然后利用广度优先搜索对细胞进行排序，建立一个二维的拓扑结构，则细胞的分化情况便被映射到对应的伪时间和他们所属的分支结构中。
步骤7：利用步骤5中的拓扑结构进行差异表达基因检测。我们选择固定的一个起点s0，利用公式(7)和(8)计算每一个基因从s0到所有终点的mic值，求它们的最大值和平均值，然后对所有基因的mic排序，排序结果越大说明该基因就对细胞分化的影响越大。通过了解这些基因影响，对有关的发育障碍和疾病能够提供有用的宝贵信息。
本领域技术人员可以理解，本发明的保护范围不局限于所述的具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征进行等同的更改或替换，需要注意的是，更改或替换之后的技术方案都将落入本发明的保护范围之内。