基于单细胞测序从头预测调控突变的系统、设备和介质

1.本发明涉及单细胞测序技术领域，具体地，涉及一种基于单细胞测序从头预测调控突变的系统、设备和介质。


背景技术：

2.随着全基因组关联研究(genome wide association study，gwas)的快速发展，与人类复杂疾病和性状具有显著关联的遗传位点的发现日益增加。然而，人们对这些遗传关联背后的生物学机制的理解远远落后于当前已有的数据库资源。在人类复杂疾病和性状的遗传关联位点中，只有小部分位于基因的编码区，而绝大部分则位于基因的非编码区。这意味着决定生物学性状的突变除了影响蛋白质的结构和功能，还可能通过影响基因的转录表达来发挥作用。位于编码区的突变功能相对容易解释，例如，外显子突变效应可以通过遗传密码规则进行评估。但是对于位于非编码区的突变功能的理解则更加困难，迄今为止，还没有一种统一的方法实现对其的精准解释。
3.表达数量性状位点(expression quantitative trait locus，eqtl)是一类能够影响基因表达水平的遗传位点，它们可以通过在全基因组尺度上对人群中的突变和基因表达水平之间的关联进行测试来鉴定。通过进一步整合eqtl与gwas的结果，许多研究将基因表达水平作为突变与生物学表型之间的关键分子表型，以揭示gwas发现的遗传关联位点的生物学功能。迄今为止，绝大多数的eqtl分析都是基于传统大样本转录组测序(bulk rna sequencing，bulk rna-seq)进行的。这种方法从百万个裂解细胞中收集rna并进行测序，所获得的基因表达数据代表的是来自样本中多个不同细胞类型的平均值。过去的许多研究表明，基因调控在不同的细胞类型之间存在非常大异质性。因此，在传统eqtl关联分析中，遗传变异对基因表达的细胞类型特异性影响可能会被掩盖。
4.细胞类型互作表达数量性状位点(cell type interaction eqtl，cell type-ieqtl)是一种依赖于特定细胞类型的表达数量性状位点。这种新兴的ieqtl算法是对传统eqtl关联分析方法的扩展。该算法首先利用去卷积分析估计组织样本中各种细胞类型的比例，然后通过在线性模型中测试基因型和估计的细胞比例，以推测细胞类型特异性的基因表达调控的遗传变异效应。这种方法在未来的科研和临床应用中表现出极大的潜力，因为它具有极高的通用性和实用性。该方法不需要研究人员额外收集大量特定细胞类型的基因组与转录组数据，只需要使用现有的组织水平遗传数据库(例如gtex数据集)以及可下载的单细胞转录图谱资源，就能通过计算ieqtl检测细胞类型水平的基因表达遗传调控，而几乎不需要额外的计算资源成本。
5.然而，由于在组织样本中估计的不同细胞类型的比例通常相互关联(或负相关)，ieqtl除了捕获测试细胞类型的特异性调控效应，还捕获了与之相关的细胞类型的特异性调控效应。这导致了ieqtl计算中一个非常大的局限：难以确定检测到的遗传调控效应的真正细胞类型来源(即因果细胞类型)。此外，由于基因组上的多个位点通常处于连锁不平衡(linkage disequilibrium，ld)，一个遗传位点的ieqtl统计值除了捕获该位点的等位基因
调控效应，还捕获了与之连锁的其他遗传位点的等位基因调控效应，因此难以确定真正具有基因表达调控功能的突变(即因果突变)。这些问题限制了ieqtl分析的精准度，阻碍了它在功能基因组学以及临床靶点转化方面的应用。


技术实现要素：

6.为了解决上述技术问题中的至少一个，本发明采用的技术方案如下：
7.本发明第一方面提供一种基于单细胞测序从头预测调控突变的系统，包括以下模块：
8.染色质特征提取模块，用于基于卷积神经网络提取基因转录起始位点
±
20kb区域的染色质特征；
9.染色质特征存储模块，与所述染色质特征提取模块连接，用于接收并存储所述染色质特征提取模块提取得到的全基因组范围内全部基因的染色质特征；
10.细胞亚群特异性基因表达谱获取模块，用于接收所述多细胞亚群的单细胞测序数据，生成细胞亚群特异性基因表达谱及相应的细胞亚群信息；
11.预测模型建立模块，分别与所述染色质特征存储模块和所述细胞亚群特异性基因表达谱获取模块连接，用于针对不同细胞亚群，分别利用回归方法，基于所述染色质特征和所述细胞亚群特异性基因表达谱训练出基于染色质特征预测基因表达水平的预测模型；
12.预测模型应用模块，与所述预测模型建立模块连接，用于接收所述预测模型建立模块训练出的针对不同亚群的预测模型；
13.所述染色质特征提取模块还直接与所述预测模型应用模块连接，用于：
14.①
接收特定细胞亚群的一个或多个基因序列，并提取基因的染色质特征输入到所述预测模型应用模块，相应地，所述预测模型应用模块还用于利用训练得到的针对相应细胞亚群的预测模型基于所述染色质特征预测基因表达水平；和/或
15.②
接收一对包含一个突变位点的等位基因，提取所述等位基因的染色质特征并计算所述等位基因的染色质特征的差值；将染色质特征的差值输入到所述预测模型应用模块，相应地，所述预测模型应用模块还用于利用训练得到的针对不同细胞亚群的预测模型基于染色质特征的差值预测所述等位基因在相应细胞亚群中的突变预测效应，若所述等位基因在某个细胞亚群中的突变预测效应的绝对值大于其它细胞亚群，则该等位基因上的突变位点为该细胞亚群的特异性调控突变位点。
16.在本发明中，“细胞亚群”是比“细胞类型”更下位的一个概念。一种“细胞类型”包括一个或多个“细胞亚群”。
17.在本发明的一些实施方案中，在预测模型建立模块中，在训练时，选择一部分基因作为训练集，另一部分基因作为测试集，进行模型训练和评估，具体地，在特定细胞亚群中利用训练集中的基因的染色质特征和基因表达水平利用线性回归模型来拟合细胞亚群特异性的基因表达水平。优选地，划分到训练集中的基因数目多于划分到测试集中的基因数目。在本发明的一些具体实施方案中，将位于8号染色体上的990个基因划分为测试集，将其余基因划分为训练集。当然本领域技术人员可以选择其他染色体的基因作为测试集，甚至随机划分基因作为训练集。若仍有部分基因没有划分为测试集或训练集，则可以利用上述
①
的方法预测其表达水平。
18.在本发明的一些实施方案中，所述等位基因是一对包括一个snp的基因，分别为参考基因序列和备选等位基因序列。相对于参考基因序列，备选等位基因序列中含有一个snp位点。
19.在本发明中，
②
中由于输入到预测模型中的是染色质特征的差值，反应的是基因表达水平的变化程度，也就是突变预测效应。因此所述突变预测效应反应突变对基因表达水平的影响程度，通过计算等位基因的突变预测效应，即可得到snp位点对基因表达的干扰效应。在本发明的一些实施方案中，当突变预测效应绝对值大于0.5时，才认为是细胞亚群的特异性调控突变位点。
20.在本发明的一些实施方案中，所述预测模型应用模块还与细胞亚群特异性基因表达谱获取模块连接，所述预测模所述染色质特征提取模块还用于：
21.③
接收一对包含一个突变位点的等位基因，提取所述等位基因的染色质特征并输入到所述预测模型应用模块；相应地，所述预测模型应用模块还用于利用训练得到的针对不同细胞亚群的预测模型基于相应的染色质特征预测所述等位基因在相应细胞亚群中的表达水平，并用于计算预测得到所述等位基因在相应细胞亚群中的表达水平的差值，若所述等位基因在某个细胞亚群中得到的表达水平的差值的绝对值大于其它细胞亚群，则该等位基因上的突变位点为该细胞亚群的特异性调控突变位点。尽管该途径的计算方法较为复杂，但在实现突变预测的效果上是等同的。
22.在本发明的一些实施方案中，所述染色质特征基于beluga提取所述染色质特征，具体地，利用2000bp的序列窗口大小作为输入，在每个200bp的基因组区域中2002个表观特征。在本发明的另一些实施方案中，所述染色质特征基于deepsea提取所述染色质特征，具体地，利用1000的序列窗口大小作为输入，在每个200bp的基因组区域中919个表观特征。在本发明的又一些实施方案中，所述染色质特征基于sei提取所述染色质特征，具体地，利用4096的序列窗口大小作为输入，在每个200bp的基因组区域中4381400个表观特征。
23.在本发明的一些实施方案中，并用于将单细胞测序数据中基因表达数进行归一化，缺失值用0进行填充，进一步对归一化后的矩阵进行填补，最后将同一细胞亚群中的所有细胞的基因表达量取平均值，得到所述细胞亚群特异性基因表达谱。在本发明的一些优选实施方案中，所述单细胞测序数据中基因表达数通过log2(cpm
i,j
/10+1)来进行归一化，其中cpm
i,j
表示细胞j中基因i的每百万计数。
24.在本发明的一些实施方案中，所述预测模型模块中，利用线性回归方法训练预测模型。在本发明的一些优选实施方案中，利用xgboost回归方法训练预测模型。在本发明的一些具体实施方案中，利用xgboost回归方法进行训练时，超参数的选择中，将booster设置为线性；将训练最大迭代次数设置为40，并使用迭代过程的早期停止以避免过拟合；将用于更新模型的步长收缩设置为0.5；将模型的损失函数设置为线性回归；将权重的l2正则化项增加为100；将其余的超参数设为默认值。在本发明的一些实施方案中，利用套索回归或岭回归方法训练预测模型。
25.本发明第二方面提供一种基于单细胞测序从头预测调控突变的设备，包括处理器和存储器，所述存储器以所述处理器可执行的计算机程序形式储存本发明第一方面任一所述的系统，所述处理器执行所述存储器储存的计算机程序，实现所述系统的功能。
26.本发明第三方面提供一种基于单细胞测序从头预测调控突变的介质，所述介质为
cluster-ieqtl的共定位基因在“银标准”中的富集。
38.图8示出了本发明实施例3中基于不同cnn预测基因表达水平的表现。该图展示了同一细胞亚群特异性模型对990个测试集基因的表达水平预测的表现。y轴为每组模型的预测和实际测量的基因表达水平之间的pcc，x轴为用于比较的几种cnn模型。
39.图9示出了本发明实施例3中基于不同cnn预测细胞亚群特异性基因表达的表现。该图展示了同一测试集基因在不同细胞亚群特异性模型中预测表达水平的表现。y轴为每组模型的990个测试集基因的预测和实际测量的表达水平之间的pcc，x轴为用于比较的几种cnn模型。
40.图10示出了本发明实施例3中基于sei的细胞亚群特异性预测模型的训练损失曲线。该图展示了在利用sei提取的染色质特征训练基因表达预测模型中，基于训练集和测试集评估的模型表现，y轴是损失函数的值，x轴是训练的轮次。
41.图11示出了本发明实施例3中基于不同的线性回归方法预测基因表达水平的表现。该图展示了同一细胞亚群特异性模型对990个测试集基因的表达水平预测的表现。y轴为每组模型的预测和实际测量的基因表达水平之间的pcc，x轴为用于比较的几种线性回归方法。
42.图12示出了本发明实施例3中基于预测模型的基因表达预测。该密度图显示了在动脉组织细胞亚群(a)成纤维细胞、(b，c)内皮细胞和(d)平滑肌细胞中，预测模型预测的和实际观测的基因表达之间的一致性。x轴表示对数转化的基因表达预测值，y轴表示基于单细胞表达数据得到的对数转化的基因表达水平。
43.图13示出了本发明实施例3中预测模型在20个组织中对基因表达的预测表现。该箱线图显示了不同组织的细胞亚群特异性模型的基因表达预测表现。x轴为不同组织的预测模型分组，y轴为预测和实际观测到的基因表达水平的pcc相关性。
44.图14示出了本发明实施例3中预测模型对基因表达水平高低的预测。该密度图显示了在动脉组织的几种细胞亚群中，当只考虑观察到表达的测试集基因时，预测模型对基因表达水平的预测表现。x轴为对数转化的基因表达预测值，y轴表示基于单细胞表达数据得到的对数转化的基因表达水平。
45.图15示出了本发明实施例3中预测模型对细胞亚群特异性基因表达的预测表现。该树状图基于每个模型中细胞亚群特异性基因表达的预测值与所有细胞亚群中细胞亚群特异性基因表达的观察值之间的pcc，描述了357个细胞亚群的关系。树上方不同的颜色表示每个细胞亚群对应的组织信息，中间不同的颜色表示对应的细胞谱系信息，右边显示了每个细胞亚群模型中预测和观察的基因表达水平的相关性
46.图16示出了本发明实施例3中基于预测模型的eqtl统计值预测。该图显示了在外周血组织中，预测模型预测的突变效应与eqtl统计值z-scores之间的相关性。每个点代表一个独立的ld block，x轴显示了该基因组区域中最大的对数转化的突变预测效应绝对值，y轴显示了该基因组区域中最大的对数转化的eqtl z-scores绝对值。
47.图17示出了本发明实施例3中不同突变组别在cell cluster-ieqtl的富集。该箱型图展示了基于预测模型预测值分类的两组突变位点在相应组织中cell cluster-ieqtl的富集。突变位点组2代表每个组织中预测的细胞亚群特异性调控位点，这些突变的细胞亚群模型的预测绝对值显著大于bulk水平模型的预测绝对值。突变位点组1代表每个组织中
突变位点组2之外的其它突变位点。每个点代表一个测试组织。
48.图18示出了本发明实施例3中预测模型预测精细映射突变位点的表现。图为基于预测模型预测的逻辑分类器判断一个突变位点是否为精细映射突变位点的roc曲线。左图和右图分别为模型在距离(a)tss 20kb以内和(b)tss 20kb以外的预测表现。
49.图19示出了本发明实施例4中整合分析的框架流程图。(a)通过使用cnn模型beluga，整合每个基因tss周围的染色质特征；(b)根据细胞图谱，为每个细胞亚群分别训练细胞亚群特异性的突变预测模型；(c)根据学生t检验对所有预测的突变进行排序，以获得对每个细胞亚群中特异性的突变位点；(d)根据人群中ld相关性对eqtl和突变位点预测结果进行整合，确定受预测到的突变位点顺式调控的候选基因。
50.图20示出了本发明实施例4整合分析示意图。根据ld相关性对预测的突变位点和eqtl进行整合分析。
51.图21示出了本发明实施例4中关联分析阈值测试。(a)该图展示了基于chromhmm软件对动脉组织中的六种不同染色质修饰数据进行整合的结果。每行代表具有相似染色质状态的一组基因组区域；(b)不同阈值下的eqtl位点在相应组织中ccres的富集结果。x轴表示不同eqtl组别，y轴为每组eqtl相比于匹配的对照组正好位于ccres的or。
52.图22示出了本发明实施例4中细胞亚群特异性基因调控遗传变异情况。该图统计了每种细胞亚群的特异性功能调控突变以及受其顺式调控的关联基因的数量。
53.图23示出了本发明实施例4中得到的细胞亚群特异性基因调控遗传变异情况与scatac-seq数据集的比较分析。该热图展示了每种细胞亚群中推断的特异性调控突变(y轴)在scatac-seq数据集中细胞亚群特异性开放染色质区域(x轴)的标准化后的富集分数。该图反映了推断的细胞亚群特异性调控变异和实际观察到的细胞亚群特异性调控区域之间的良好一致性。x轴和y轴的标签颜色代表了细胞亚群所属的谱系。
54.图24示出了本发明实施例4中cell cluster-ieqtl与外部数据集的比较分析。该热图展示了每种细胞亚群中cell cluster-ieqtl位点(y轴)在scatac-seq数据集中细胞亚群特异性开放染色质区域(x轴)的标准化后的富集分数。相比之下，cell cluster-ieqtl位点与实际观察到的细胞亚群特异性调控区域之间的一致性较差。x轴和y轴的标签颜色代表了细胞亚群所属的谱系。
具体实施方式
55.除非另有说明、从上下文暗示或属于现有技术的惯例，否则本技术中所有的份数和百分比都基于重量，且所用的测试和表征方法都是与本技术的提交日期同步的。在适用的情况下，本技术中涉及的任何专利、专利申请或公开的内容全部结合于此作为参考，且其等价的同族专利也引入作为参考，特别这些文献所披露的关于本领域中的相关术语的定义。如果现有技术中披露的具体术语的定义与本技术中提供的任何定义不一致，则以本技术中提供的术语定义为准。
56.为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。
57.实施例
58.以下例子在此用于示范本发明的优选实施方案。本领域内的技术人员会明白，下
述例子中披露的技术代表发明人发现的可以用于实施本发明的技术，因此可以视为实施本发明的优选方案。但是本领域内的技术人员根据本说明书应该明白，这里所公开的特定实施例可以做很多修改，仍然能得到相同的或者类似的结果，而非背离本发明的精神或范围。
59.除非另有定义，所有在此使用的技术和科学的术语，和本发明所属领域内的技术人员所通常理解的意思相同，在此公开引用及他们引用的材料都将以引用的方式被并入。
60.那些本领域内的技术人员将意识到或者通过常规试验就能了解许多这里所描述的发明的特定实施方案的许多等同技术。这些等同将被包含在权利要求书中。
61.下述实施例中的实验方法，如无特殊说明，均为常规方法。下述实施例中所用的仪器设备，如无特殊说明，均为实验室常规仪器设备；下述实施例中所用的试验材料，如无特殊说明，均为自常规生化试剂商店购买得到的。
62.实施例1单细胞测序及数据处理
63.1.单细胞测序
64.获得与gtex匹配的20个成年人体组织类型，共39例样本，详细信息如表1所示。
65.表1本实施中分析的组织及相关信息
[0066][0067]
利用微孔板测序技术(microwell-seq)(han,xiaoping et al.mapping the mouse cell atlas by microwell-seq,cell,volume 172,issue 5,1091-1107.e17)对各样本进行高通量单细胞转录组测序(scrna-seq)。对microwell-seq数据进行质控，去除rna背景污染并过滤潜在双细胞双细胞后，得到高质量的单细胞数据集。
[0068]
2.整合高质量的单细胞数据集
[0069]
在对scrna-seq数据进行背景污染矫正和潜在双细胞过滤后，发明人进一步对处理后的数据进行无监督聚类，具体步骤如下：
[0070]
(1)对于检测基因数大于1000的高质量单细胞数据，将相同组织的批次的基因表达矩阵利用merge函数合并在一起。使用r包seurat中的函数createseuratobject创建seurat对象，在该函数中，将最小细胞设置为3(min.cells＝3)，最小基因设置为200(min.features＝200)。
[0071]
(2)在创建的seurat对象meta.data中，输入hcl原文注释的详细细胞类型信息。
[0072]
(3)利用函数normalizedata(seurat包)对seurat对象中的基因表达矩阵进行归一化。其中，将归一化方法设置为lognormalize，缩放因子设置为10000。该步骤将每个基因的计数除以该细胞总的基因计数，再乘以指定的缩放因子10000，最后使用log1p进行自然对数转换。
[0073]
(4)使用函数findvariablefeatures(seurat包)寻找每个seurat对象中的高变基因。在该函数中，将方法设置为vst，将高变基因数量设置为2000。该步骤首先使用局部多项式拟合log转化的基因表达方差与均值之间的关系，然后使用观测均值和由拟合线给出的期望方差对特征值进行标准化，最后基于标准化后的特征值计算方差。
[0074]
3.数据降维处理和亚群定义
[0075]
在整合每个组织的单细胞数据集后，发明人进一步对其进行了降维和分群分析，具体步骤如下：
[0076]
(1)使用scaledata函数(seurat包)对整合的seurat对象中的基因表达进行中心化，使用runpca(seurat包)对高变基因进行主成分分析(principal component analysis，pca)，并参考elbowplot的曲线拐点选择后续使用的pca维度。
[0077]
(2)在runtsne(seurat包)函数中，通过选取的pca维度使用fit-tsne进行降维。在该函数中，将dims设置为选取的pca数。
[0078]
(3)使用findneighbors和findclusters函数(seurat包)，基于细胞间的距离进行分群。将findneighbors中k.param设置为10，维度(dims)设置为选择的pca维度；将findclusters中迭代次数设置为50，分辨率设置为1。
[0079]
(4)使用dimplot函数(seurat包)来对每个组织的分群情况进行可视化。在该函数中，将reduction设置为“fitsne”。
[0080]
(5)通过findallmarkers函数(seurat包)，利用wilcoxon秩和检验来计算每个细胞亚群的特征高表达基因(marker)。
[0081]
(6)将hcl原文中的注释投射至高质量数据的细胞亚群中。
[0082]
最后，发明人对细胞亚群进行归类，将其划分为不同的主要细胞亚群以及细胞谱系。在本实施例中，一共获得357个细胞亚群，其中每个组织中包含的细胞亚群数量详见表2。
[0083]
表2单个组织中包含的细胞亚群数量
[0084][0085]
发明人将这357个细胞亚群进一步归类为44个主要细胞类型和9个不同的谱系，具体信息见表3。
[0086]
表3高质量细胞图谱的主要细胞类型
[0087]
[0088][0089]
实施例2确定细胞亚群互作eqtl的方法
[0090]
目前，eqtl通常是基于大块样本的测序数据，并利用传统的遗传关联分析方法推测得到的。这些被推测的eqtl通常被称为标准eqtl。标准eqtl可以从gtex网站上获得。在标准eqtl中，细胞亚群特异性的基因调控作用可能被组织样本中的细胞异质性所掩盖，尤其是对于那些在样本中比例较小的细胞亚群来说。
[0091]
为了实现在细胞亚群水平上对基因表达遗传调控的探索，发明人扩展了传统的eqtl关联分析方法，提出利用单细胞表达数据来推测细胞亚群互作eqtl(cell cluster interaction eqtl，cell cluster-ieqtl)。发明人将cell cluster-ieqtl定义为依赖于细胞图谱中的细胞亚群与基因表达水平产生关联的表达数量性状位点。如图1所示，在单核苷酸多态性(single nucleotide polymorphism，snp)遗传位点上具有更高t等位基因剂量的样本中，基因表达水平与样本中某细胞亚群的比例呈正相关。此外，在含有更高该细胞亚群的组织样本中，snp基因型与基因表达之间具有更强相关性。这些观察表明，该snp位点是一个细胞亚群特异性的基因表达调控位点，因此被判定为cell cluster-ieqtl。
[0092]
本实施例中，cell cluster-ieqtl确定方法的主要步骤分为大块样本的去卷积分析、大块样本表达矩阵的归一化、样本中协变量的获取和交互效应的检测。
[0093]
1.大块样本的去卷积分析
[0094]
(1)对于scrna-seq数据中每一个群，随机抽取500个细胞(小于500个细胞的群则不抽样)，并将抽取的细胞每10个细胞作为一个假细胞(pseudo-cell)。
[0095]
(2)求取假细胞中各基因表达水平的平均值，获得假细胞-基因表达参考矩阵。
[0096]
(3)在软件cibersortx中，分析模块选择“create signature matrix”，分析模式选择“custom”，输入数据类型选择“scrna-seq”，然后将上述假细胞-基因表达参考矩阵上传至自定义输入中的“single cell reference matrix file”处，作为单细胞参考矩阵文件。此外在附加选项“single cell input options”中，将基因表达阈值“min.expression”设置为0。其他选项默认，以此创建基因表达特征矩阵。
[0097]
(4)将大块样本rna-seq测序数据上传至cibersortx作为混合文件，并利用s-mode(single cell mode)来矫正不同来源的数据之间的批次效应。
[0098]
(5)在绝对模式下对大块样本进行去卷积分析，利用基因表达特征矩阵来计算每个样本中的细胞亚群的绝对分数。绝对模式下生成的绝对分数将相对细胞分数缩放成反映混合物中每种细胞亚群的绝对比例的分数，虽然绝对分数没有直接表示为细胞占比，但它可以用于在不同细胞亚群之间直接比较所估计的细胞比例大小。保留在超过1/5的样本中绝对分数大于1的细胞亚群用于下游分析。
[0099]
同时，为避免异常值的影响，对计算的绝对分数进行逆正态变换(inversely normally transformed，int)。
[0100]
2.大块样本表达矩阵的归一化
[0101]
对大块样本rna-seq测序数据进一步地进行tmm(trimmed mean of m-values)归一化，保留在20％的样本中tpm(transcripts per million)不少于0.1以及未归一化的原始序列数不小于6的基因。
[0102]
3.样本中协变量的获取
[0103]
利用peer软件从tmm归一化之后的表达矩阵中提取能够解释基因表达变化的隐藏因子。peer推断的隐藏因子(peer factor)和人群表达数据中存在的许多混杂因素显著相关，包括技术的协变量(如文库构建批次、样本缺血时间、检测到转录本的数量
……
)和生物学个体的协变量(如年龄、性别、种族
……
)。因此，在cell cluster-ieqtl确定过程中，将peer factor作为协变量可以减少结果的假阳性，提高cell cluster-ieqtl位点的检测能力。peer因子的数量选择取决于cell cluster-ieqtl关联测试中所用到的样本数量。
[0104]
利用plink2软件，对测试样本的snp基因型进行主成分分析，并从突变矩阵中提取前五个主成分(principal components，pcs)。
[0105]
将提取的peer factor、基因型主成分、文库构建批次、样本缺血时间、检测到转录本的数量、测序平台、年龄、性别、种族等其它已知协变量合并在一起，然后检测并丢弃共线性的协变量。
[0106]
4.交互效应的检测
[0107]
在细胞图谱中的每个细胞亚群中，对于tmm归一化中保留的全部基因、以及它们转录起始位点(tss)周围1mb区域中突变频率(maf)大于0.05的突变(snp)，在tensorqtl软件中根据以下线性模型进行cell cluster-ieqtl的计算：
[0108]
tmm归一化基因表达～β1
×
snp基因型:细胞亚群比例+2
×
snp基因型+3
×
4，如图7)。总的来说，这些结果显示了本实施例推测的cell cluster-ieqtl结果在生物学上是合理的。
[0119]
实施例3从头预测细胞亚群水平的突变
[0120]
本实施例提供基于基于实施例1中scrna-seq数据构建细胞亚群特异性深度学习模型的流程图。主要步骤分为染色质状态bulk水平的预测、细胞亚群特异性基因表达谱的获取、细胞亚群特异性线性回归模型的训练和基因表达预测表现的评估。
[0121]
1.染色质状态bulk水平的预测
[0122]
利用cnn模型提取目前已知的24339个蛋白质编码、核糖体rna(ribosomal rna，rrna)和长链非编码rna(long noncoding rna,lncrna)的dna序列(基因)的染色质特征(包括组蛋白标记、转录因子结合和dnase超敏反应谱等表观数据集)。
[0123]
本实施例中，发明人选择三种cnn模型进行染色质特征提取：
[0124]
cnn模型deepsea(2015年版本)能够基于dna序列预测919个染色质特征；cnn模型beluga能够基于序列预测2002个染色质特征；cnn模型sei能够基于序列预测来自超过1300种组织和细胞亚群的21907个染色质特征。利用以上三种cnn模型，发明人分别预测了24339个基因的tss
±
20kb基因组区域的染色质特征。对于deepsea、beluga和sei三种模型，发明人分别使用1000bp、2000bp和4096bp的序列窗口大小作为输入，在每个200bp的基因组区域中预测了919、2002和21907个表观特征。因此，对于每个基因，deepsea、beluga和sei的预测分别生成了183800(919
×
200)、400400(2002
×
200)和4381400(21907
×
200)个染色质特征。
[0125]
2.细胞亚群特异性基因表达谱的获取
[0126]
1)对于24339个基因，单细胞测序数据中基因表达原始计数矩阵通过log2(cpm
i,j
/10+1)来进行归一化，其中cpm
i,j
表示细胞j中基因i的每百万计数。缺失值用0来进行填充。
[0127]
2)使用magic算法对以上得到的归一化表达矩阵进行填补，来弥补由于单细胞测序导致的矩阵稀疏性。需要说明的是，单细胞rna-seq数据通常包含许多由于原始rna扩增失败而导致的缺失值或dropout，因此低表达的基因往往捕获不到，如果不进行imputation，可能导致后续模型学习及评估不够充分。
[0128]
3)为了进一步增加测得的基因数量和umi数量，发明人对scrna-seq数据中同一个群的所有细胞的基因表达量取平均值，生成细胞亚群特异性的基因转录谱。
[0129]
3.细胞亚群特异性线性回归模型的训练
[0130]
为了捕获细胞图谱中不同细胞亚群的表观特征，发明人基于cnn模型提取的每个基因tss周围的染色质状态，利用线性回归模型来拟合细胞亚群特异性的基因表达水平。在模型训练中，将位于8号染色体上的990个基因划分为测试集，将其余基因划分为训练集，利用训练集进行模型训练。
[0131]
4.基因表达预测表现的评估
[0132]
1)对于每个细胞亚群特异性模型，通过计算在990个测试集基因中的预测表达水平和实际测量表达水平的皮尔森相关性(pearson’s correlation coefficient，pcc)，来检测模型在捕获不同基因的表达水平差异方面的表现。
[0133]
2)对于每个测试集基因，通过计算在不同scrna-seq数据细胞亚群中的预测表达水平和实际测量表达水平的pcc值，来检测模型在捕获细胞亚群特异性基因表达水平方面
的表现。
[0134]
5.基于不同cnn和线性回归方法的模型性能评估
[0135]
在构建细胞亚群特异性的突变预测模型中，为了选择适合的cnn和线性回归模型，发明人针对不同的cnn模型和不同的线性回归模型，分别进行了基因表达预测，从而评估不同的模型选择对基因表达预测的影响。
[0136]
(1)基于不同cnn的模型性能评估
[0137]
然后，发明人按照前述步骤和方法，根据不同cnn提取的染色质特征，分别拟合人类细胞图谱(human cell landscape，hcl)中不同细胞亚群的基因表达水平，由此构建了三组细胞亚群特异性的预测模型，均利用xgboost线性回归模型来拟合细胞亚群特异性的基因表达水平，在超参数的选择中，将booster设置为线性；将训练最大迭代次数设置为40，并使用迭代过程的早期停止以避免过拟合(通过设置early_stopping_rounds＝10)；将用于更新模型的步长收缩设置为0.5；将模型的损失函数设置为线性回归；将权重的l2正则化项增加为100；将其余的超参数设为默认值。最后对不同组别的模型预测表现进行评估和比较。
[0138]
结果显示，在捕获不同基因表达水平的差异方面，基于beluga的模型表现优于基于deepsea和sei的模型(如图8)。在捕获不同细胞亚群基因表达水平的差异方面，基于beluga训练的模型同样也取得了最好的表现(如图9)。虽然基于sei的模型使用了最多的染色质特征来预测基因表达水平，但其表现最差。分析发现，这是由于sei提取的染色质特征中存在大量与基因表达无关的噪音信息，加上训练样本的有限，模型出现了过拟合(如图10)。
[0139]
(2)基于不同线性回归方法的模型性能评估
[0140]
由于预测基因表达的染色质因子特征数量非常大(每个基因400400个特征)，并且用于训练的基因样本数量相对有限(24339个基因)，因此使用线性回归模型拟合基因表达时可能会遇到多重共线性和过拟合的潜在问题。为了选择适合的线性回归方法，发明人利用beluga提取染色质特征，分别利用普通最小二乘法线性回归(ordinary least-squares linear regression)、套索回归(lasso regression)、岭回归(ridge regression)以及xgboost梯度提升方法分别拟合hcl中不同细胞亚群的基因表达水平，从而构建得到四组细胞亚群特异性的预测模型。发明人对不同组别的模型预测表现进行评估和比较。
[0141]
结果显示，在几种不同的线性回归模型中，基于xgboost的预测模型表现最好，其次是基于套索回归和岭回归的模型，而基于普通最小二乘法线性回归的预测表现最差(如图11)。
[0142]
此外，发明人还对基于不同的线性回归方法训练细胞亚群特异性预测模型所需要的时间进行了比较。结果如表4所示，可以看出，xgboost在训练细胞亚群特异性预测模型的速度方面显著优于套索回归和岭回归。
[0143]
表4利用不同的回归方法训练基因表达预测模型的时间
[0144][0145]
因此，xgboost能够更好地处理该线性回归问题，并更快地为scrna-seq数据批量建立细胞亚群特异性的预测模型。
[0146]
综上所述，本实施例利用beluga模型提取染色质特征，进一步利用xgboost回归拟合基因表达水平，从而建立了细胞亚群特异性的突变预测模型。
[0147]
5.突变预测模型性能评估
[0148]
发明人分别构建了来自20个组织的357个细胞亚群的特异性突变预测模型。为了评估模型的性能，发明人首先计算了每个模型预测的测试集基因表达水平与实际测量的表达水平之间的相关性。
[0149]
结果表明，这些细胞亚群特异性模型能够有效地预测相应细胞亚群中的基因表达水平(p＜2.2e-16，图12)。总体来说，在357个细胞亚群的模型中，测试集基因的预测和实际观察到的基因表达水平之间的pcc中位数达到了0.763。而在肾脏、胃和横结肠的某些细胞亚群模型中，pcc甚至超过了0.80(图13)。为了进一步证明该模型可以准确预测基因表达水平的变化，而不仅仅是二元的转录状态(表达或不表达)，发明人只保留在每个细胞亚群中观察到表达的测试集基因，并重新计算预测和实际观察到的表达水平之间的一致性。结果表明，预测模型仍然显示出对基因表达良好的预测性能。在357个细胞亚群模型中，pcc中位数达到了0.63(图14)。
[0150]
此外，发明人还评估了预测模型在预测细胞亚群特异性基因表达方面的能力。具体而言，发明计算了每个细胞亚群中预测的基因表达水平与其他所有细胞亚群中观察到的基因表达水平相对差异的对数值，并将其作为细胞亚群特异性基因表达水平的衡量指标。随后，发明人计算了每个模型中细胞亚群特异性基因表达的预测值与所有细胞亚群中细胞亚群特异性基因表达的观察值之间的pcc，并以此对细胞亚群进行了无监督聚类。结果表明，预测模型能够准确地预测不同细胞亚群的谱系特异性基因表达，同一谱系内的基因表达预测值具有更高的一致性，而不同谱系之间的基因表达预测值差异较大(图15)。
[0151]
6.具有特异性调控功能的突变位点的预测
[0152]
首先，通过计算机模拟dna碱基替换，得到包括假定突变位点的基因序列作为备选等位基因序列。其次，计算参考基因序列(野生型基因序列)和备选等位基因序列在beluga输出的染色质特征的差值。接着，将该差值输入训练好的细胞亚群特异性的xgboost模型，以预测突变在特定细胞亚群中对基因表达的突变预测效应(干扰效应)。然后，根据单侧学生t检验(student's t test)，选择top 2000并且突变预测效应绝对值大于0.5的位点，测试特定细胞亚群中突变预测效应的绝对值是否大于其它细胞亚群，并据此来确定每种细胞亚群的特异性调控突变位点。
[0153]
发明人还先将等位基因序列的染色质特征输入到输入训练好的细胞亚群特异性的xgboost模型，以预测等位基因序列在特定细胞亚群中的表达水平。同样，根据单侧学生t检验(student's t test)，测试特定细胞亚群中等位基因的表达水平与参考基因序列在相
应细胞亚群中的表达水平的差值的绝对值是否大于其它细胞亚群，并据此来确定每种细胞亚群的特异性调控突变位点。
[0154]
7.突变预测结果的合理性评估
[0155]
为了检测突变效应预测结果的合理性，发明人首先测试了预测模型是否可以仅基于dna序列重现基于人群数据的关联分析结果。发明人将来自每个组织中所有细胞亚群模型的预测结果汇总在一起，并与标准eqtl进行比较。结果显示，在同一个连锁不平衡区域(ld block，r2≥0.6)中，eqtl统计值z-scores(效应大小β/标准误)与突变预测效应的最高绝对值表现出显著的相关性(p＜2.2e-16，图16)。在20个测试的组织中，两者的斯皮尔曼相关系数(spearman’s correlation coefficient，scc)的中位数为0.3。
[0156]
接着，发明人测试了预测模型能否仅基于dna序列识别cell cluster-ieqtl位点。需要注意的是，在一个研究的组织中，本领域技术人员通常无法确定cell cluster-ieqtl的基因调控效应实际上来源于何种细胞亚群。这是因为不同细胞亚群在测试样本中估计的比例是线性相关的，cell cluster-ieqtl的调控效应既可能来源于测试snp基因型交互效应的细胞亚群，也可能来源于与之线性相关的其他细胞亚群。尽管如此，发明人发现，在同一组织中，基于所有细胞亚群的预测模型推测的遗传调控位点在所有细胞亚群的cell cluster-ieqtl中显著富集(在18/20组织中，bonferroni矫正的p值小于0.05)。
[0157]
然后，发明人利用大块组织样本的基因表达数据为20个组织分别训练了bulk水平模型，并对突变效应进行了重新预测。对于每个组织，发明人比较了细胞亚群水平的突变预测和bulk水平的突变预测结果。发明人按照细胞亚群模型的预测效应绝对值是否大于bulk水平模型的预测效应绝对值的1.96个标准差，将突变位点分为两组。结果发现，那些在细胞亚群模型中具有明显更高预测绝对值的突变位点，同样在cell cluster-ieqtl中也具有更为显著的富集(图17)。这些结果表明，在揭示细胞亚群特异性遗传调控方面，预测模型的预测与cell cluster-ieqtl具有一致性。
[0158]
此外，发明人还测试了预测模型能否直接指出基因调控的因果突变位点。由于人群中等位基因ld的存在，与基因具有最显著关联的突变位点并不一定是真正导致该基因表达水平改变的突变。eqtl因果突变位点的鉴定可以更加准确地进行基因功能分析，促进对疾病基础生物学的了解和药物靶标的鉴定。预测模型的突变预测基于dna序列碱基替换实现，因此理论上其结果不受ld的影响。发明人使用组织中每个变异体的最高预测绝对值，与gtex组织的精细映射eqtl数据集进行比较。结果表明，预测模型可以在相互具有较高ld关联(r2≥0.6)的突变位点中，实现区分精细映射的阳性和阴性突变，其平均auroc为0.780(距离tss 20kb以内auroc＝0.782，20kb以外auroc＝0.766，图18)。因此，预测模型可以识别因果突变方面，在近端和远端调控区域中都表现出良好的性能。
[0159]
为了进一步检验预测模型的突变预测能力，发明人在人永生化红细胞系(hudep-2)中训练了预测模型，并对hbg基因tss周围1kb基因组区域的突变进行了逐一地预测。结果表明，致病突变的预测效应大小是与它们匹配的随机选择的突变集合的3.25倍(95％ ci＝1.62～6.16)。这些结果表明预测模型能够在一定程度上识别非编码致病突变，从而为其在单碱基水平的因果突变预测提供了有力的证据。
[0160]
实施例4预测受突变位点顺式调控的基因
[0161]
由于基因组调控元件可以通过三维空间结构(例如dna环状结构)的协同作用来影
响基因的转录表达，仅仅基于dna序列可能无法充分估计突变对于特定基因表达水平的影响，并且也缺乏进一步的遗传学证据支持。因此，发明人通过整合实施例2得到的cell cluster-ieqtl和实施例3预测的突变位点，确定受突变位点顺式调控的基因。如图19所示，整合分析(荟萃分析)的具体步骤如下：
[0162]
(1)获取标准eqtl，并利用实施例2的方法步骤获得cell cluster-ieqtl。为了获得标准eqtl在全基因组范围的fdr，发明人采用了与cell cluster-ieqtl多重假设检验矫正方法相同的方法。具体而言，使用bh方法对所有eqtl的基因水平的显著性进行多重假设检验矫正，以获得在全基因组范围内的fdr，选取全基因组显著性fdr＜0.4作为标准eqtl和cell cluster-ieqtl的阈值进行后续分析。
[0163]
(2)基于标准eqtl和cell cluster-ieqtl计算中使用的snp基因型矩阵，使用plink来估计遗传位点的连锁不平衡(ld)相关性。
[0164]
(3)在每个
±
1mb的测试区域中，当一种细胞亚群测试出的eqtl或cell cluster-ieqtl与同一种细胞亚群预测出的特异性调控突变位点呈高连锁不平衡(ld r2》0.8)，则将该qtl或cell cluster-ieqtl关联的基因推测为受该预测的突变位点顺式调控的基因(如图20所示)。
[0165]
在步骤(1)中，为了选择合适的标准eqtl和cell cluster-ieqtl的fdr阈值，以便进行后续的整合分析，发明人对动脉组织进行了基准测试，并将其与包括启动子和增强子在内的候选顺式调节元件(ccres)进行了比较。发明人整合了动脉组织中不同染色质组蛋白修饰的染色质免疫沉淀后测序(chromatin immunoprecipitation sequencing，chip-seq)数据，并将染色质划分为e1-e10 10个不同的染色质状态(如图21a)。发明人最终确定将基因组区域e1、e2和e4作为ccres，其中e1和e2为增强子，e4为启动子。
[0166]
接着，发明人将动脉组织的qtl或cell cluster-ieqtl分为以下四组：a)fdr＜0.05且连锁位点(r2≥0.8)预测效应的绝对值大于0.5；b)fdr＜0.4且连锁位点(r2≥0.8)预测效应的绝对值大于0.5；c)fdr＜0.05；d)fdr＜0.4。发明人分别估计了每组突变位点在e1、e2、e4基因组区域中的富集情况以及95％的置信区间。结果发现，当只考虑关联分析时，fdr＜0.05组的突变相比fdr＜0.4组的突变在ccres中具有更为显著的富集(如图21b)。然而，在联合从头预测进行考虑时，情况出现了变化，fdr＜0.4相比fdr＜0.05组具有更为显著的富集。根据该结果，发明人认为整合分析本身可以避免假阳性的发生。因此，在整合分析流程中，使用fdr＜0.4作为阈值更为适合。该阈值能够在真阳性和假阳性数量之间提供了一个合理的平衡，在保证准确度的同时，也能够推测出更多基因调控遗传突变位点。
[0167]
在本实施例中，发明人利用实施例3的方法预测13182个细胞亚群特异性功能突变位点，利用本实施例的方法预测到6181基因的表达水平受到细胞亚群特异性功能突变顺式调控的影响(图22)。
[0168]
为了测试本实施例的方法是否正确地识别了各个细胞亚群中的特异性遗传调控效应，发明人将突变位点-基因关联关系与一个外部的scatac-seq数据集(zhang,k.et al.a single-cell atlas of chromatin accessibility in the human genome.cell 184,5985-6001e5919)进行比较，并计算了各个细胞亚群中推断的调控位点在实际检测的细胞亚群特异性开放染色质区域的富集程度。结果显示，相同谱系的细胞亚群具有类似的富集情况，而基于富集结果的无监督层次聚类能够反映这些细胞亚群之间众所周知的谱系
关系(图23)。同时，还在同一谱系的细胞亚群中，基因调控位点在差异开放染色质峰中表现出特异性富集。因此，该富集分析表明，本实施例推测的细胞亚群特异性调控位点与实际观察到的细胞亚群特异性调控区域之间具有一致性。
[0169]
此外，由于cell cluster-ieqtl也是细胞亚群依赖的遗传调控位点，发明人使用类似的方法计算了这些位点在scatac-seq测量的细胞亚群特异性开放染色质区域的富集程度，并进行了无监督聚类分析。相比之下，该富集结果显示出较低的细胞谱系特异性。例如，一种细胞亚群的cell cluster-ieqtl位点除了在相同谱系的细胞亚群的差异开放染色质峰中富集外，也在其他谱系的细胞亚群中表现出明显的富集(图24)。这是因为，在一种细胞亚群中检测到的cell cluster-ieqtl可能实际上来源于其他细胞亚群中发生的基因调控。因此，相比单纯的cell cluster-ieqtl分析，本实施例的整合分析作为一个集成框架能够更好地识别感兴趣细胞亚群的特异性基因调控位点。
[0170]
在本发明中，可以利用计算机模块实现上述各步骤的功能，各个模块集合或组装和配合成系统，进一步利用存储器以处理器可执行的计算机程序形式储存系统。还可利用计算机可读存储介质存储上述系统。
[0171]
在本发明提及的所有文献都在本技术中引用作为参考，就如同每一篇文献被单独引用作为参考那样。此外应理解，在阅读了本发明的上述讲授内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本技术所附权利要求书所限定的范围。

技术特征：
1.一种基于单细胞测序从头预测调控突变的系统，其特征在于，包括以下模块：染色质特征提取模块，用于基于卷积神经网络提取基因转录起始位点
±
20kb区域的染色质特征；染色质特征存储模块，与所述染色质特征提取模块连接，用于接收并存储所述染色质特征提取模块提取得到的全基因组范围内全部基因的染色质特征；细胞亚群特异性基因表达谱获取模块，用于接收所述多细胞亚群的单细胞测序数据，生成细胞亚群特异性基因表达谱及相应的细胞亚群信息；预测模型建立模块，分别与所述染色质特征存储模块和所述细胞亚群特异性基因表达谱获取模块连接，用于针对不同细胞亚群，分别利用回归方法基于所述染色质特征和所述细胞亚群特异性基因表达谱训练出基于染色质特征预测基因表达水平的预测模型；预测模型应用模块，与所述预测模型建立模块连接，用于接收所述预测模型建立模块训练出的针对不同亚群的预测模型；所述染色质特征提取模块还直接与所述预测模型应用模块连接，用于：
①
接收特定细胞亚群的一个或多个基因序列，并提取基因的染色质特征输入到所述预测模型应用模块，相应地，所述预测模型应用模块还用于利用训练得到的针对相应细胞亚群的预测模型基于所述染色质特征预测基因表达水平；和/或
②
接收一对包含一个突变位点的等位基因，提取所述等位基因的染色质特征并计算所述等位基因的染色质特征的差值；将染色质特征的差值输入到所述预测模型应用模块，相应地，所述预测模型应用模块还用于利用训练得到的针对不同细胞亚群的预测模型基于染色质特征的差值预测所述等位基因在相应细胞亚群中的突变预测效应，若所述等位基因在某个细胞亚群中的突变预测效应的绝对值大于其它细胞亚群，则该等位基因上的突变位点为该细胞亚群的特异性调控突变位点。2.根据权利要求1所述的一种基于单细胞测序从头预测调控突变的系统，其特征在于，所述预测模型应用模块还与细胞亚群特异性基因表达谱获取模块连接，所述染色质特征提取模块还用于：
③
接收一对包含一个突变位点的一个等位基因，提取所述等位基因的染色质特征并输入到所述预测模型应用模块；相应地，所述预测模型应用模块还用于利用训练得到的针对不同细胞亚群的预测模型基于相应的染色质特征预测所述等位基因在相应细胞亚群中的表达水平，并用于计算预测得到所述等位基因在相应细胞亚群中的表达水平的差值，若所述等位基因在某个细胞亚群中得到的表达水平的差值的绝对值大于其它细胞亚群，则该等位基因上的突变位点为该细胞亚群的特异性调控突变位点。3.根据权利要求1或2所述的一种基于单细胞测序从头预测调控突变的系统，其特征在于，所述染色质特征基于beluga提取所述染色质特征，具体地，利用2000bp的序列窗口大小作为输入，在每个200bp的基因组区域中2002个表观特征。4.根据权利要求1或2所述的一种基于单细胞测序从头预测调控突变的系统，其特征在于，并用于将单细胞测序数据中基因表达数进行归一化，缺失值用0进行填充，进一步对归一化后的矩阵进行填补，最后将同一细胞亚群中的所有细胞的基因表达量取平均值，得到所述细胞亚群特异性基因表达谱。5.根据权利要求4所述的一种基于单细胞测序从头预测调控突变的系统，其特征在于，
所述单细胞测序数据中基因表达数通过log2(cpm
i,j
/10+1)来进行归一化，其中cpm
i,j
表示细胞j中基因i的每百万计数。6.根据权利要求1或2所述的一种基于单细胞测序从头预测调控突变的系统，其特征在于，所述预测模型模块中，利用l2正则化的线性回归方法训练预测模型。7.根据权利要求6所述的一种基于单细胞测序从头预测调控突变的系统，其特征在于，利用xgboost回归方法训练预测模型。8.根据权利要求7所述的一种基于单细胞测序从头预测调控突变的系统，其特征在于，利用xgboost回归方法进行训练时，超参数的选择中，将booster设置为线性；将训练最大迭代次数设置为40，并使用迭代过程的早期停止以避免过拟合；将用于更新模型的步长收缩设置为0.5；将模型的损失函数设置为线性回归；将权重的l2正则化项增加为100；将其余的超参数设为默认值。9.一种基于单细胞测序从头预测调控突变的设备，其特征在于：包括处理器和存储器，所述存储器以所述处理器可执行的计算机程序形式储存所述权利要求1-8任一项所述的系统，所述处理器执行所述存储器储存的计算机程序，实现所述系统的功能。10.一种基于单细胞测序从头预测调控突变的介质，其特征在于：所述介质为计算机可读存储介质，并储存所述权利要求1-8任一项所述的系统，计算机读取并执行所储存的系统。

技术总结
本发明公开了一种基于单细胞测序从头预测调控突变的系统、设备和介质，属于单细胞测序技术领域。所述系统利用卷积神经网络，提取基因转录起始位点周围的染色质特征，并通过回归方法，由染色质特征拟合细胞亚群特异性的基因表达，为每细胞亚群分别训练特异性的预测模型；进一步利用预测模型，预测细胞亚群特异性调控突变位点。本发明的系统能够显著提升预测结果的精准度，使其成为一种研究细胞亚群水平遗传调控的通用的方法，在细胞功能基因组学领域具有巨大的应用价值。域具有巨大的应用价值。域具有巨大的应用价值。


技术研发人员：郭国骥 肖彦宇 王晶晶 李佳琦
受保护的技术使用者：浙江大学
技术研发日：2023.05.23
技术公布日：2023/7/26