一种可解释的基于基因选择的细胞数据分析方法与流程

1.本发明属于单细胞rna序列分析领域，具体涉及一种可解释的基于基因选择的细胞数据分析方法。


背景技术：

2.单细胞rna测序（sc-rna seq）的兴起，使得在分子水平上能够检测细胞类型，其在分析糖尿病、阿尔兹海默病和癌症的病理学中发挥了独特的作用，包括揭示患病组织中基因表达的动态，剖析高度异质细胞的功能或功能障碍，以及分析癌细胞的演变和转移。而聚类是sc-rna seq分析中的关键步骤，可以识别隐藏的细胞亚型，推断细胞的拓扑排列，描述复杂疾病中的细胞异质性。然而，由于细胞rna序列数据维度极高，具有大量噪声，同时具有复杂的数据结构和属性，这大大增加了聚类的难度。
3.现有的单细胞rna序列分析技术多有缺陷。由于数据具有复杂的空间分布，巨大的噪声和测量误差，seurat方法（一种分析单细胞转录组的方法）和sc3（single cell consensus clustering，单细胞一致性聚类）方法，难以实现高于80%的聚类精度。而基于深度学习的聚类方法，由于其具有极强的表示性而显著提高了聚类精度，但其问题依然存在。其一，深度学习是黑箱模型，其操作的可解释性较差。其二，深度学习通过梯度下降训练大量的参数用于数据的良好表示，这一过程非常耗时。其三，深度学习需要针对不同规模，不同类型的细胞设计网络的结构和超参数，这限制了深度学习框架的可扩展性。


技术实现要素：

4.为解决现有技术的不足，实现提高基因聚类精度，更准确的度量细胞间相似度，提升序列分析的可解释性、计算效率和可扩展性的目的，本发明采用如下的技术方案：
5.一种可解释的基于基因选择的细胞数据分析方法，包括如下步骤：
6.步骤一，通过基因选择，得到一组细胞数据；具体地，根据基因最值和方差初步选择基因；
7.步骤二，通过计算细胞间序距离，并基于序距离构造细胞数据的邻接矩阵；
8.步骤三，通过邻接矩阵计算分类伪标签，并基于基因与伪标签的同质性进行进一步基因选择；
9.步骤四，为进一步基因选择的细胞数据计算细胞孤立系数，并基于孤立系数去除异常细胞；
10.步骤五，将去除异常的细胞数据，通过序距离构建的无向图转化出新的邻接矩阵，新的邻接矩阵基于无向图中的连通分支计算出分类标签，同类标签再次构建无向图，并计算前后无向图的一致性，基于一致性得到细胞类数；
11.步骤六，根据细胞类数和稀疏邻接矩阵对细胞进行聚类，得到聚类结果。
12.进一步地，所述步骤二中，根据两个细胞的序距离及其他细胞分别与两个细胞的序距离比，构建邻接矩阵。
13.进一步地，所述步骤三包括如下步骤：
14.步骤3.1，计算邻接矩阵的正规化矩阵，并对正规化矩阵进行特征值分解；
15.步骤3.2，基于最大的一组特征值对应的特征向量，构建谱投影；
16.步骤3.3，对谱投影的行向量进行分类，得到对应的伪标签；
17.步骤3.4，获取细胞数据中基因对应的行向量，计算行向量关于伪标签的量子化；
18.步骤3.5，计算量子化的基因与伪标签的同质性，得到细胞数据关于伪标签的互信息；
19.步骤3.6，基于互信息对应的基因，得到进一步基因选择的细胞数据。
20.根据步骤三的方法进行进一步的基因选择，能够选择几乎所有的类指示基因。
21.进一步地，所述步骤四包括如下步骤：
22.步骤4.1，为进一步基因选择的细胞数据构建距离矩阵；
23.步骤4.2，基于距离矩阵中，与细胞距离最近的一组细胞对应的元素，计算细胞的紧密系数；
24.步骤4.3，将紧密系数最小的一组细胞作为异常细胞并移除。
25.进一步地，所述步骤五包括如下步骤：
26.步骤5.1，在去除异常的细胞数据中，将细胞作为顶点，基于序距离进行顶点连接，构建第一无向图，并将第一无向图转化为第一图邻接矩阵；基因选择和异常细胞处理后，构造的系数邻接矩阵具有更准确、更强的类内连接，和更稀疏、更弱的类间连接；
27.步骤5.2，寻找无向图的所有连通分量，满足连通分量内顶点之间存在一条连接路径，且连通分量间不存在任何边；
28.步骤5.3，计算第一图邻接矩阵的正规化矩阵，并对正规化矩阵进行特征值分解；
29.步骤5.4，基于连通分量个数设定最大的一组特征值，通过特征值对应的特征向量，构建谱投影；
30.步骤5.5，对谱投影的行向量进行分类，得到对应的标签；
31.步骤5.6，对于标签对应的同类细胞，建立顶点连接，构建第二无向图，并将第二无向图转化为第二图邻接矩阵；
32.步骤5.7，获取第一无向图和第二无向图的一致性；
33.步骤5.8，基于标签对应的所述一致性，以及标签对应的分类数，选取聚类类数，得到细胞类数。
34.进一步地，所述细胞间序距离的计算如下：
35.通过基因表达水平对细胞进行排列，基于排列后的位置，计算细胞间的相关性，根据相关性构建细胞数据的距离矩阵；
36.确定两个细胞的间距，从两个细胞周围，分别获取距离小于间距的细胞数，将少的细胞数作为两个细胞间的序距离。相比于传统的欧氏距离或斯皮尔曼spearman距离，序距离更准确表达了细胞的高维分布，为后续分析过程奠定坚实基础。
37.进一步地，基于其他细胞分别与两个细胞的序距离对邻接矩阵稀疏化。构造的稀疏邻接矩阵具有弱类间连接和强类内连接。
38.进一步地，所述步骤六中，根据聚类结果和稀疏邻接矩阵对细胞数据得到不同分辨率的可视化；
39.所述方法还包括步骤七，根据聚类结果与可视化结果，选取不同分辨率下不同类型细胞的标记基因。
40.进一步地，所述步骤六中，设定细胞数据的初始二维可视化样本，并定义样本间有理相似度，基于有理相似度定义样本与所述邻接矩阵间引力散度、斥力散度，并将可视化目标优化问题转化为同时最小化引力散度和斥力散度，采用梯度下降优化样本，选择所述新的邻接矩阵，求解可视化目标优化问题得到低分辨率可视化结果，对新的邻接矩阵稀疏化，进而求解可视化目标优化问题得到高分辨率可视化结果。
41.进一步地，所述步骤七包括如下步骤：
42.步骤7.1，对于聚类结果中某类细胞的标记基因，分别计算初步基因选择后并去除异常细胞的数据集在该类和其他类别上某个基因的均值与方差，基于均值与方差计算该基因与聚类结果单变量方差分析的p值，选取一组最大的p值对应的基因为最终选取的该类细胞的标记基因；
43.步骤7.2，针对高分辨率可视化结果，计算样本间距离矩阵；
44.步骤7.3，构造高分辨率无向图，图中每一个顶点对应可视化结果中的一个行向量，即对应一个细胞的二维投影，两个顶点间边的连接根据样本间距离阈值设定，定义无向图的连通分量满足连通分量内顶点之间存在一条连接路径，且连通分量间不存在任何边，寻找无向图的所有连通分量结果；
45.步骤7.4，将连通分量结果作为聚类结果，选取该聚类结果中某类细胞的标记基因，分别计算初步基因选择后并去除异常细胞的数据集在该类和其他类别上某个基因的均值与方差，基于均值与方差计算该基因与聚类结果单变量方差分析的p值，选取一组最大的p值对应的基因为最终选取的该类细胞的标记基因，得到高分辨率下不同类型细胞的标记基因。
46.本发明的优势和有益效果在于：
47.相比于传统方法，本发明能够筛选出有利于得到更准确聚类结果的基因并更准确的度量细胞之间的相似度，实现具有明显优势的高精度单细胞rna序列数据聚类以及更能反映数据结构的数据可视化方法；相比于近年兴起的深度学习方法，本发明提供了具有可解释性，计算快速简单，可扩展性优异的分析框架；本发明能够在多种类型的单细胞rna序列数据集上，实现简单，快速，准确的种类识别和可视化。
附图说明
48.图1是本发明实施例中方法的流程图。
49.图2a是本发明实施例中数据集yan上欧氏距离、spearman距离和序距离对应轮廓系数的示意图。
50.图2b是本发明实施例中数据集goolam上欧氏距离、spearman距离和序距离对应轮廓系数的示意图。
51.图3是本发明实施例中数据集yan中在进一步基因选择步骤前构造的稀疏邻接矩阵的示意图。
52.图4a是本发明实施例中数据集yan中两次基因选择步骤时剩余基因数目对比的示意图。
53.图4b是本发明实施例中数据集goolam中两次基因选择步骤时剩余基因数目对比的示意图。
54.图5a是本发明实施例中数据集yan中选择基因与未被选择基因的类内、类间方差对比图。
55.图5b是本发明实施例中数据集goolam中选择基因与未被选择基因的类内、类间方差对比图。
56.图6是本发明实施例中数据集yan在基因选择和异常细胞处理后稀疏邻接矩阵的示意图。
57.图7是本发明实施例中数据集yan选择聚类类数与其他类数聚类精度ari对比的示意图。
58.图8a是本发明实施例中数据集yan上，本发明和scdha的聚类精度ari稳定性对比的示意图。
59.图8b是本发明实施例中数据集goolam上，本发明和scdha的聚类精度ari稳定性对比的示意图。
60.图9a是本发明实施例中数据集usoskin上低分辨率可视化结果的示意图。
61.图9b是本发明实施例中数据集usoskin上高分辨率可视化结果的示意图。
62.图10a是本发明实施例中数据集usoskin上高分辨率np3类型细胞的标记基因gas2可视化结果的示意图。
63.图10b是本发明实施例中数据集usoskin上高分辨率nf4/5类型细胞的标记基因wdr16可视化结果的示意图。
64.图10c是本发明实施例中数据集usoskin上高分辨率pep1类型细胞的标记基因prg2可视化结果的示意图。
65.图10d是本发明实施例中数据集usoskin上高分辨率nf1类型细胞的标记基因mir704可视化结果的示意图。
66.图11是本发明实施例中装置的结构示意图。
具体实施方式
67.以下结合附图对本发明的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明，并不用于限制本发明。
68.一种可解释的基于基因选择的细胞数据分析方法，首先通过计算细胞间无量纲的序距离构造稀疏相似度矩阵，并利用谱投影得到伪标签，进而比对伪标签和基因之间相似性，筛选出少于300个基因用于聚类，并对比细胞邻域内相似度去除异常细胞，然后使用稀疏邻接矩阵的连通分量个数，以及计算伪标签和邻接矩阵的匹配度，得到聚类类数，进而得到聚类结果。最后通过求解改进交叉熵损失函数得到二维可视化投影，完成单细胞rna序列分析。
69.为了使本发明的目的、技术方案和技术效果更加清楚明白，以下结合说明书附图，以人类胚胎数据集yan、小鼠胚胎数据集goolam、小鼠大脑皮层细胞数据集usoskin为例，为本发明作进一步相似说明。如图1所示，序列分析方法具体包括如下步骤：
70.步骤一，选取细胞序列数据中的基因，根据基因最值和方差初步选择基因；对于单
细胞rna序列数据，去除基因最大值小于log23的基因，并保留基因方差大于1.5的基因，得到细胞数据c={c1,
…
,cn}，其中ci表示数据中第i个细胞，n为细胞个数，设此时基因个数为m1。
71.步骤二，计算细胞间序距离并构造稀疏邻接矩阵，包括如下步骤：
72.步骤2.1，对于细胞ci和cj，i≠j，将细胞ci和cj中基因表达水平（表达量）从高到低排列，记排列后的位置为c’i
和c’j
，进而计算细胞之间斯皮尔曼spearman相关系数s(ci,cj)为：
[0073][0074]
其中||
·
||2为向量的二范数，定义为：
[0075][0076]
构造斯皮尔曼spearman距离矩阵d
(1)
，使得d
(1)
的第i行第j列元素为d
ij(1)
=s(ci,cj)。
[0077]
步骤2.2，定义细胞cj关于细胞ci的序order(cj;ci)为所有细胞中距离细胞ci的斯皮尔曼spearman距离小于细胞cj到细胞ci的斯皮尔曼spearman距离的细胞个数；同时定义细胞ci关于细胞cj的序order(ci;cj)为所有细胞中距离细胞cj的斯皮尔曼spearman距离小于细胞ci到细胞cj的斯皮尔曼spearman距离的细胞个数；进而计算细胞ci与细胞cj间的序距离为：
[0078][0079]
定义的序距离具有更好表达细胞空间分布的特性，如图2a、图2b所示，横轴为轮廓系数，其取值范围为-1到1，越接近1说明距离越准确表述细胞空间分布，纵轴为累积细胞个数定义为达到该轮廓系数时的细胞个数。图2a、图2b中绘制了欧氏距离、斯皮尔曼spearman距离和序距离三者的累计概率随轮廓系数的变化，其中点线为斯皮尔曼spearman距离，虚线为欧氏距离，粗实线于斯皮尔曼spearman距离的序距离。图2a、图2b中序距离最慢，意味着更多的细胞在此距离下有更高的轮廓系数，即相比于传统的欧氏距离或斯皮尔曼spearman距离，序距离更准确表达了细胞的高维分布，为后续分析过程奠定坚实基础。
[0080]
对于已知距离矩阵d，使其第i行第j列元素d
ij
为ci与cj之间的距离，设真实标签对应细胞分划为{c1,
…
,ck}，对于样本ci，定义其类内距离a(i)类间距离b(i)为：
[0081][0082]
以及的轮廓系数为：
[0083][0084]
其中card{
·
}为集合中元素的个数，cm表示细胞真实分划中第m类细胞，cn表示细胞真实分划中第n类细胞；
[0085]
步骤2.3，构造邻接矩阵w，其第i行，第j列元素为w
ij
定义为：
[0086]
[0087]
其中代表其他细胞到细胞ci的序距离中从小到大排列中排行第7位的值，代表其他细胞到细胞cj的序距离中从小到大排列中排行第7位的值。进而将邻接矩阵w稀疏化为当序距离满足：
[0088][0089]
此时设置对应的w
ij
=0，将此时得到的稀疏邻接矩阵仍记为。构造的稀疏邻接矩阵具有弱类间连接和强类内连接，如图3所示，yan数据集的真实分类结果共有六类，每一个方框内代表真实的类内连接，从左上到右下分别为第一类到第六类，黑色区域代表强连接，白色区域代表无连接，每一个方框内代表真实的类内连接，可以看出yan数据集前三类之间均具有类间连接，而第四类的类内连接较弱。
[0090]
步骤三，通过邻接矩阵计算伪标签并基于伪标签进行进一步基因选择，包括如下步骤：
[0091]
步骤3.1，计算邻接矩阵的正规化矩阵，其中度矩阵，为以向量为对角线的对角矩阵，并对进行特征值分解。
[0092]
步骤3.2，选取的最大个特征值对应特征向量，构成谱投影，这里参数分别取3、4、5，得到的谱投影分别记为。
[0093]
步骤3.3，使用余弦距离下的k-means算法分别将的行向量分为3类，4类，5类，得到的分类结果分别记为伪标签。并重复运行k-means算法10次以消除算法本身的随机性。
[0094]
步骤3.4，将细胞数据的每一个基因对应的行向量作为新视角，设细胞数据，其中为第个基因。首先定义向量整数化函数为：
[0095][0096]
其中表示向量的第个分量，为下取整符号，即选取不大于的最大整数。定义基因的标准化为：
[0097][0098]
即将基因的取值范围控制在0到1之间，进而定义基因关于标签的量子化为：
[0099][0100]
计算每一个基因关于伪标签的量子化，分别记为。
[0101]
步骤3.5，设同样长度的整数向量的取值个数分别为，记分划为：
[0102][0103]
其中为向量取第个值时的元素指标，为向量取第个值时的元素指标。记：
[0104][0105]
其中card为集合元素个数。整数向量之间同质性ari为：
[0106][0107]
其中，分别定义为
[0108][0109]
分别计算整数化基因和伪标签的同质性，并将细胞数据关于伪标签互信息结果记为。对于伪标签，选取集合中最大的个元素对应的位置。这里参数取为100，分别取3，4，5。
[0110]
步骤3.6，记，保留中元素对应的基因，记进一步基因选择的细胞数据为。
[0111]
这种智能基因选择方法首先只保留了不到300个基因，如图4a、图4b所示，例如yan数据集和goolam数据集中原有20214和41428个基因，在基因过滤，即步骤一后剩余8063和8641个基因，但在基因选择后仅保留208和190个基因，仅占整体基因的1%和0.5%；另外，步骤三选择了几乎所有的类指示基因，如图5a、图5b所示，设细胞的真实分类结果为，定义第类细胞在第个基因上的类内方差为：
[0112][0113]
其中为类别中细胞数目，为类别中第个细胞在第个基因的表达水平，为类别中在第个基因的表达水平平均值。此时可以定义个基因上的类内方差为：
[0114][0115]
定义第类细胞在第个基因上的类间方差为：
[0116][0117]
这里为第中细胞个数，分别为类别中在第个基因的表达水平平均值。定义个基因上的类内方差为：
[0118][0119]
图5a、图5b中横坐标代表基因类内方差，纵坐标代表基因类间方差，当一个基因的类内方差小、类间方差大时，可以认为该基因是某一类或某几类细胞的标记基因，图中绘制了yan数据集和goolam数据集中步骤三选择基因与其他基因类内、类间的方差，其中空心圆点代表本发明选择的基因，散点代表本发明未选择的基因。可以说明，步骤三中选择的基因几乎包含了所有小类内方差，大类间方差的标记基因。
[0120]
步骤四，定义细胞孤立系数并根据孤立系数去除异常细胞，具体包含如下四个步骤：
[0121]
步骤4.1，构造的斯皮尔曼spearman距离矩阵使得矩阵的第行，第列元素为。
[0122]
步骤4.2，计算细胞的紧密系数为：
[0123][0124]
其中表示与细胞的斯皮尔曼spearman距离最近的个细胞（不包括），这里参数设置为10。
[0125]
步骤4.3，设置()最小的个细胞为异常细胞并移除，这里参数设置为，其中为上取整。设去除异常细胞后的数据仍记为，n1表示去除异常细胞后的细胞个数；
[0126]
步骤4.4，提取斯皮尔曼spearman距离矩阵中去除异常细胞的部分仍记为，并如步骤2.2中定义序距离方法，计算去除异常细胞后数据集对应序距离矩阵，使得其第行第列元素为：
[0127][0128]
步骤五，再次构造稀疏邻接矩阵并基于连通分支估计细胞类数，具体包括如下八个步骤：
[0129]
步骤5.1，构造无向图，图中第个顶点对应细胞，如果细胞间距离满足：
[0130][0131]
此时连接顶点和顶点，反之则不连接，其中代表集合中第小的值，这里参数选取为3。将无向图转化为图邻接矩阵仍记为，使得其第行第列元素为：
[0132][0133]
并构造另一个稀疏邻接矩阵，使得其第行第列元素为：
[0134][0135]
其中代表集合中第小的值，这里参数选取为7。进而将邻接矩阵稀疏化为当序距离满足：
[0136][0137]
设置，将此稀疏邻接矩阵仍记为。基因选择和异常细胞处理后，构造的系数邻接矩阵具有更准确、更强的类内连接，和更稀疏、更弱的类间连接，如图6所示，图6为yan数据集的稀疏邻接矩阵，黑色区域代表强连接，白色区域代表无连接，每一个方框内代表真实的类内连接，可以看出yan数据集前三类之间均具有类间连接，而第四类的类内连接较弱。相比于图3中，yan数据集在基因选择和异常细胞处理前的稀疏邻接矩阵，一类到四类之间的类间连接明显减弱，第五类细胞的类内连接加强，使其不再具有分块结构。
[0138]
步骤5.2，定义无向图的连通分量满足连通分量内顶点之间存在一条连接路径，且连通分量间不存在任何边，寻找无向图的所有连通分量，并记连通分量个数为。
[0139]
步骤5.3，计算邻接矩阵的正规化矩阵，其中。并对特征值分解。
[0140]
步骤5.4，选取的最大个特征值对应特征向量，构成谱投影，这里参数分别取，以及谱投影分别记为。
[0141]
步骤5.5，使用余弦距离下的并行k-means算法分别的行向量分为nc类，nc+1类，nc+2类，nc+3类，得到的分类结果分别记为标签。并重复运行k-means算法10次以消除算法本身的随机性。
[0142]
步骤5.6，基于标签，构造无向图，图中第个顶点对应细胞，如果细胞满足第个细胞和第个细胞在为同类细胞，连接图中第个顶点和第个顶点，反之则不连接。进而将无向图转化为图邻接矩阵仍记为，使得其第行，第列元素为：
[0143][0144]
这里参数分别取。基于邻接矩阵，构造无向图，图中第个顶点对应细胞，如果满足，连接图中第个顶点和第个顶点，反之则不连接。进而将无向图转化为图邻接矩阵仍记为，使得其第行，第列元素为：
[0145][0146]
步骤5.7，定义无向图在无向图下的平均同质性为：
[0147][0148]
平均异质性为：
[0149][0150]
这里为元素全为1的矩阵。为矩阵内积，定义为：
[0151][0152]
从而定义无向图在无向图下的平均一致性为：
[0153][0154]
并计算分别取下的。
[0155]
步骤5.8，设定备选聚类类数集合为，选取聚类类数为：
[0156][0157]
步骤五准确找到了数据集的真实分类结果，yan数据集中估计类数，真实类数为6类，并且通过比较yan数据集上估计类数与其他类数对应聚类结果的聚类精度如图7
所示，精度指标为ari，直线指向标记为步骤五选取聚类类数对应精度，可以得出基于稀疏邻接矩阵在步骤五估计的聚类类数条件下，得到最准确的聚类结果。
[0158]
步骤六，按照细胞类数和稀疏邻接矩阵对细胞进行聚类，并根据聚类结果和稀疏邻接矩阵，对细胞数据得到不同分辨率的可视化，具体包括以下五个步骤：
[0159]
步骤6.1，设，其中为矩阵的第个行向量，应用余弦距离下的并行k-means算法将的行向量分为类，并将结果作为去除异常细胞后的单细胞rna序列数据集的聚类结果和谱投影结果。并重复运行k-means算法10次以消除算法本身的随机性。事实上基于稀疏邻接矩阵的谱聚类稳定性高于其他方法，以最新颖的scdha方法为例，两个方法的稳定性比较如图8a、图8b所示，通过在yan数据集和goolam数据集重复运行本发明与scdha方法，计算每一次聚类结果与真实分类标签的聚类精度ari，并统计为箱型图表示，可以看出，本发明在20次重复实验中结果稳定，保持不变，而scdha方法产生了较大震荡，说明了本发明的聚类稳定性。
[0160]
步骤6.2，将二维单位圆周平均分为段，记节点坐标为，对于在标签中样本属于第类，并设第类中共个细胞，m表示第j类样本的数量。定义权重为：
[0161][0162]
定义在二维的投影为：
[0163][0164]
这里参数设置为二维标准正态分布生成的随机数。设置单细胞rna序列数据的初始二维可视化为。
[0165]
步骤6.3，对于样本集，定义样本间有理相似度为：
[0166][0167]
其中为自适应参数。构造有理相似度矩阵使得其第行第列元素为。这里参数设置为通过最小化初始可视化的样本间有理相似度与高斯核函数相似度得到，具体优化问题为：
[0168][0169]
其中参数设置为0.1，为元素全为1的矩阵，设置为矩阵内积，为样本集的欧氏距离矩阵，第行第列元素为：
[0170][0171]
为矩阵逐个元素指数操作，即有：
[0172][0173]
其中为矩阵的第行第列元素。
[0174]
步骤6.4，对于给定邻接矩阵满足其第行第列元素为，定义样本与邻接
矩阵间引力kl散度为：
[0175][0176]
其中为元素逐个元素对数操作，即有
[0177][0178]
其中为矩阵的第行第列元素。定义样本与邻接矩阵间斥力kl散度为：
[0179][0180]
进而定义可视化目标优化问题为同时最小化引力kl散度和斥力kl散度为：
[0181][0182]
求解可视化目标优化问题中，为加速求解过程，首先将可视化目标优化问题改写为固定其他样本优化求解每一个样本，即定义样本与样本间引力kl散度为：
[0183][0184]
定义样本与样本间斥力kl散度为：
[0185][0186]
最终定义样本更新样本的可视化优化函数为：
[0187][0188]
采用随机梯度下降思想，选取样本,计算其关于的引力梯度为：
[0189][0190]
以及斥力梯度为：
[0191][0192]
优化更新时，首先选取个引力样本，选取方式具体为在到范围内生成不等于的随机整数，并按照服从0到1间均匀分布的随机数生成器生成0到1间的随机数，若,则通过梯度更新，具体为：
[0193][0194]
其中为第n次迭代的下降步长。反之则重新在到范围内生成不等于的随机整数，直到选取到个引力样本，其中为截断函数，具体定义为：
[0195][0196]
然后选取个斥力样本，选取方式具体为在到范围内生成不等于的随机整数，并按照服从0到1间均匀分布的随机数生成器生成0到1间的随机数，若,则通过梯度更新，具体为：
[0197][0198]
其中为第n次迭代的下降步长。反之则重新在到范围内生成不等于的随机整数，直到选取到个斥力样本。为进一步加速梯度下降过程，考虑并行更新算法，并行更新时，均采用前一步迭代过程的输出样本计算引力梯度与斥力梯度，从而实现同时更新多
个样本，并保证一次迭代过程中，样本选择的不重复性，完成所有样本的单次更新后，称可视化过程完成一次迭代，实验中共完成500次迭代，设置第n次迭代的下降步长为：
[0199][0200]
步骤6.5，首先选择邻接矩阵为步骤5.1中构造邻接矩阵，进而求解可视化目标优化问题得到低分辨率可视化结果，求解过程中引力样本和斥力样本选取数目具体为：
[0201][0202]
进一步对邻接矩阵稀疏化为使得第行第列元素满足：
[0203][0204]
其中均为步骤5.1中定义。进而求解可视化目标优化问题得到高分辨率可视化结果，求解过程中引力样本和斥力样本选取数目具体为。
[0205]
为了说明不同分辨率可视化结果与细胞真实类型的关联，绘制了usoskin数据集在两个分辨率下的可视化结果，如图9a、图9b所示，usoskin数据集共有四类，分别为nf，np，pep，th，在低分辨率可视化结果中，本发明除将np数据集分为两个子集外，其余类型细胞空间分布均被准确还原，而在高分辨率可视化结果中，进一步将nf分为三个亚型，分别为nf1，nf2/3，nf4/5；np分为三个亚型，np1和np2和np3，pep分为两个亚型，分别为pep1和pep2。达到usoskin报告的亚型细胞识别分辨率，为其他单细胞rna序列数据集亚型细胞识别提供直观表示。
[0206]
步骤七，根据聚类结果与可视化结果选取不同分辨率下不同类型细胞的标记基因，具体包括以下四个步骤：
[0207]
步骤7.1，假设选聚类结果中第类细胞的标记基因，分别计算初步基因选择后并去除异常细胞的数据集在第类和其他类别上第个基因的均值与方差，分别记为，并计算第个基因与聚类结果单变量方差分析的p值，初次选取第类细胞的标记基因集合为：
[0208][0209]
最后选取前30个最大值对应基因为最终选取第类细胞的标记基因。
[0210]
步骤7.2，针对高分辨率可视化结果，计算欧氏距离矩阵，使得其第行第j列元素定义为：
[0211][0212]
其中为向量二范数。
[0213]
步骤7.3，构造高分辨率图，其中图中每一个顶点对应可视化结果中的一个行向量，即对应一个细胞的二维投影。两个顶点间有边连接，若满足：
[0214][0215]
其中代表集合中第小的值，这里取值为10。定义无向图的连通分量满足连通分量内顶点之间存在一条连接路径，且连通分量间不存在任何边，寻找无向图的所有连通分量，记连通分量结果为。
[0216]
步骤7.4，假设选聚类结果中第类细胞的标记基因，分别计算初步基因选择后并去除异常细胞的数据集在第类和其他类别上第个基因的均值与方差，分别记为，并计算第个基因与聚类结果单变量方差分析的p值，初次选取第类细胞的标记基因集合为：
[0217][0218]
最后选取前30个最大值对应基因为最终选取第类细胞的标记基因。至此得到了高分辨率下不同类型细胞的标记基因。为说明本发明在亚型细胞标记基因选取的有效性，绘制了本发明在usoskin数据集中识别的亚型细胞上的标记基因，如图10a至图10d所示，图中点为usoskin数据集在本发明下的高分辨率可视化结果，每一个点代表一个细胞，点的形状中正方形点代表该细胞中基因高表达，散点代表该细胞中基因低表达，图10a中画出的标记基因为np3类型细胞的标记基因gas2；图10b为nf4/5类型细胞的标记基因wdr16；图10c为pep1类型细胞的标记基因prg2；图10d为nf1类型细胞的标记基因mir704。可以看出，步骤七中高分辨率标记基因识别有效识别了亚型分辨率下的标记基因。
[0219]
以下举出多个实例，体现本发明中方法的优势，首先对8个公开的单细胞rna序列数据集进行聚类，所述数据集包括：人类植入前胚胎和人类胚胎干细胞rna序列数据集yan，小鼠植入前胚胎干细胞rna序列数据集goolam和deng，成人大脑皮层细胞rna序列数据集darmanis，小鼠腰椎背根神经节感觉神经元细胞rna序列数据集usoskin，人类胰岛细胞rna序列数据集xin和muraro，人类大脑皮层神经元细胞rna序列数据集lake。比较本发明方法与今年常用的单细胞rna序列分析方法：使用分层编码器的单细胞rna序列数据分析方法（scdha），加权最近邻单细胞rna序列数据分析方法（seurat），基于sps的单细胞rna序列数据分析方法（dropclust），基于核的相似度学习的单细胞rna序列数据分析方法（simlr），通过推断和降维的单细胞rna序列数据分析方法（cidr）。首先比较单细胞rna序列分析中的聚类步骤，评价的标准为聚类准确度ari，结果见表1，可得到本发明的聚类准确度均值表现最佳，且在大部分数据集上明显优于其他所有方法。
[0220]
表1 聚类准确度ari表
[0221][0222]
表2统计了这些方法在上述数据集中进行rna序列分析中的可视化步骤得到的轮廓系数。可以看出本发明的可视化步骤的结果具有更分离的异类细胞和更聚集的同类细胞，更还原细胞在高维中的分布。
[0223]
表2 轮廓系数表
[0224][0225]
表1和表2同时说明了，相比于以往的单细胞rna序列数据分析方法，本发明的方法在不同类型的数据集上，细胞聚类任务和可视化任务同时具有高性能的特点。
[0226]
与前述一种可解释的基于基因选择的细胞数据分析方法的实施例相对应，本发明还提供了一种可解释的基于基因选择的细胞数据分析装置的实施例。
[0227]
参见图11，本发明实施例提供的一种可解释的基于基因选择的细胞数据分析装置，包括存储器和一个或多个处理器，存储器中存储有可执行代码，所述一个或多个处理器执行所述可执行代码时，用于实现上述实施例中的一种可解释的基于基因选择的细胞数据分析方法。
[0228]
本发明一种可解释的基于基因选择的细胞数据分析装置的实施例可以应用在任意具备数据处理能力的设备上，该任意具备数据处理能力的设备可以为诸如计算机等设备或装置。装置实施例可以通过软件实现，也可以通过硬件或者软硬件结合的方式实现。以软件实现为例，作为一个逻辑意义上的装置，是通过其所在任意具备数据处理能力的设备的处理器将非易失性存储器中对应的计算机程序指令读取到内存中运行形成的。从硬件层面而言，如图11所示，为本发明一种可解释的基于基因选择的细胞数据分析装置所在任意具备数据处理能力的设备的一种硬件结构图，除了图11所示的处理器、内存、网络接口、以及非易失性存储器之外，实施例中装置所在的任意具备数据处理能力的设备通常根据该任意具备数据处理能力的设备的实际功能，还可以包括其他硬件，对此不再赘述。
[0229]
上述装置中各个单元的功能和作用的实现过程具体详见上述方法中对应步骤的实现过程，在此不再赘述。
[0230]
对于装置实施例而言，由于其基本对应于方法实施例，所以相关之处参见方法实施例的部分说明即可。以上所描述的装置实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本发明方案的目的。本领域普通技术人员在不付出创造性劳动的情况下，即可以理解并实施。
[0231]
本发明实施例还提供一种计算机可读存储介质，其上存储有程序，该程序被处理
器执行时，实现上述实施例中的一种可解释的基于基因选择的细胞数据分析方法。
[0232]
所述计算机可读存储介质可以是前述任一实施例所述的任意具备数据处理能力的设备的内部存储单元，例如硬盘或内存。所述计算机可读存储介质也可以是任意具备数据处理能力的设备的外部存储设备，例如所述设备上配备的插接式硬盘、智能存储卡（smart media card，smc）、sd卡、闪存卡（flash card）等。进一步的，所述计算机可读存储介质还可以既包括任意具备数据处理能力的设备的内部存储单元也包括外部存储设备。所述计算机可读存储介质用于存储所述计算机程序以及所述任意具备数据处理能力的设备所需的其他程序和数据，还可以用于暂时地存储已经输出或者将要输出的数据。
[0233]
以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的范围。

技术特征：
1.一种可解释的基于基因选择的细胞数据分析方法，其特征在于包括如下步骤：步骤一，通过基因选择，得到一组细胞数据；步骤二，通过计算细胞间序距离，并基于序距离构造细胞数据的邻接矩阵；步骤三，通过邻接矩阵计算分类伪标签，并基于基因与伪标签的同质性进行进一步基因选择；步骤四，为进一步基因选择的细胞数据计算细胞孤立系数，并基于孤立系数去除异常细胞；步骤五，将去除异常的细胞数据，通过序距离构建的无向图转化出新的邻接矩阵，新的邻接矩阵基于无向图中的连通分支计算出分类标签，同类标签再次构建无向图，并计算前后无向图的一致性，基于一致性得到细胞类数；步骤六，根据细胞类数和邻接矩阵对细胞进行聚类，得到聚类结果。2.根据权利要求1所述的一种可解释的基于基因选择的细胞数据分析方法，其特征在于：所述步骤二中，根据两个细胞的序距离及其他细胞分别与两个细胞的序距离比，构建邻接矩阵。3.根据权利要求1所述的一种可解释的基于基因选择的细胞数据分析方法，其特征在于：所述步骤三包括如下步骤：步骤3.1，计算邻接矩阵的正规化矩阵，并对正规化矩阵进行特征值分解；步骤3.2，基于最大的一组特征值对应的特征向量，构建谱投影；步骤3.3，对谱投影的行向量进行分类，得到对应的伪标签；步骤3.4，获取细胞数据中基因对应的行向量，计算行向量关于伪标签的量子化；步骤3.5，计算量子化的基因与伪标签的同质性，得到细胞数据关于伪标签的互信息；步骤3.6，基于互信息对应的基因，得到进一步基因选择的细胞数据。4.根据权利要求1所述的一种可解释的基于基因选择的细胞数据分析方法，其特征在于：所述步骤四包括如下步骤：步骤4.1，为进一步基因选择的细胞数据构建距离矩阵；步骤4.2，基于距离矩阵中，与细胞距离最近的一组细胞对应的元素，计算细胞的紧密系数；步骤4.3，将紧密系数最小的一组细胞作为异常细胞并移除。5.根据权利要求1所述的一种可解释的基于基因选择的细胞数据分析方法，其特征在于：所述步骤五包括如下步骤：步骤5.1，在去除异常的细胞数据中，将细胞作为顶点，基于序距离进行顶点连接，构建第一无向图，并将第一无向图转化为第一图邻接矩阵；步骤5.2，寻找无向图的所有连通分量，满足连通分量内顶点之间存在一条连接路径，且连通分量间不存在任何边；步骤5.3，计算第一图邻接矩阵的正规化矩阵，并对正规化矩阵进行特征值分解；步骤5.4，基于连通分量个数设定最大的一组特征值，通过特征值对应的特征向量，构建谱投影；步骤5.5，对谱投影的行向量进行分类，得到对应的标签；步骤5.6，对于标签对应的同类细胞，建立顶点连接，构建第二无向图，并将第二无向图
转化为第二图邻接矩阵；步骤5.7，获取第一无向图和第二无向图的一致性；步骤5.8，基于标签对应的所述一致性，以及标签对应的分类数，选取聚类类数，得到细胞类数。6.根据权利要求1所述的一种可解释的基于基因选择的细胞数据分析方法，其特征在于：所述细胞间序距离的计算如下：通过基因表达水平对细胞进行排列，基于排列后的位置，计算细胞间的相关性，根据相关性构建细胞数据的距离矩阵；确定两个细胞的间距，从两个细胞周围，分别获取距离小于间距的细胞数，将少的细胞数作为两个细胞间的序距离。7.根据权利要求1所述的一种可解释的基于基因选择的细胞数据分析方法，其特征在于：基于其他细胞分别与两个细胞的序距离对邻接矩阵稀疏化。8.根据权利要求1所述的一种可解释的基于基因选择的细胞数据分析方法，其特征在于：所述步骤六中，根据聚类结果和邻接矩阵对细胞数据得到不同分辨率的可视化；所述方法还包括步骤七，根据聚类结果与可视化结果，选取不同分辨率下不同类型细胞的标记基因。9.根据权利要求8所述的一种可解释的基于基因选择的细胞数据分析方法，其特征在于：所述步骤六中，设定细胞数据的初始二维可视化样本，并定义样本间有理相似度，基于有理相似度定义样本与所述邻接矩阵间引力散度、斥力散度，并将可视化目标优化问题转化为同时最小化引力散度和斥力散度，采用梯度下降优化样本，选择所述新的邻接矩阵，求解可视化目标优化问题得到低分辨率可视化结果，对新的邻接矩阵稀疏化，进而求解可视化目标优化问题得到高分辨率可视化结果。10.根据权利要求8所述的一种可解释的基于基因选择的细胞数据分析方法，其特征在于：所述步骤七包括如下步骤：步骤7.1，对于聚类结果中某类细胞的标记基因，分别计算初步基因选择后并去除异常细胞的数据集在该类和其他类别上某个基因的均值与方差，基于均值与方差计算该基因与聚类结果单变量方差分析的p值，选取一组最大的p值对应的基因为最终选取的该类细胞的标记基因；步骤7.2，针对高分辨率可视化结果，计算样本间距离矩阵；步骤7.3，构造高分辨率无向图，图中每一个顶点对应可视化结果中的一个行向量，即对应一个细胞的二维投影，两个顶点间边的连接根据样本间距离阈值设定，定义无向图的连通分量满足连通分量内顶点之间存在一条连接路径，且连通分量间不存在任何边，寻找无向图的所有连通分量结果；步骤7.4，将连通分量结果作为聚类结果，选取该聚类结果中某类细胞的标记基因，分别计算初步基因选择后并去除异常细胞的数据集在该类和其他类别上某个基因的均值与方差，基于均值与方差计算该基因与聚类结果单变量方差分析的p值，选取一组最大的p值对应的基因为最终选取的该类细胞的标记基因，得到高分辨率下不同类型细胞的标记基因。

技术总结
本发明公开了一种可解释的基于基因选择的细胞数据分析方法，选择一组细胞数据；计算细胞间序距离，并基于序距离构造细胞数据的邻接矩阵；通过邻接矩阵计算分类伪标签，并基于基因与伪标签的同质性进行基因选择；为细胞数据计算细胞孤立系数，并基于孤立系数去除异常细胞；细胞数据通过序距离构建的无向图转化出新的邻接矩阵，新的邻接矩阵基于无向图中的连通分支计算出分类标签，同类标签再次构建无向图，并计算前后无向图的一致性，基于一致性得到细胞类数；根据细胞类数和邻接矩阵对细胞进行聚类，根据聚类结果和邻接矩阵，对细胞数据得到不同分辨率的可视化；根据聚类结果与可视化结果，选取不同分辨率下不同类型细胞的标记基因。基因。基因。


技术研发人员：倪天昊 张鑫愉 李冰杰 金开秀 李洪佳
受保护的技术使用者：杭州木攸目医疗数据有限公司
技术研发日：2023.08.31
技术公布日：2023/10/7