CN116564410A

CN116564410A - 一种预测突变位点顺式调控基因的方法、设备和介质

Info

Publication number: CN116564410A
Application number: CN202310599596.0A
Authority: CN
Inventors: 郭国骥; 肖彦宇; 王晶晶
Original assignee: Zhejiang University ZJU
Current assignee: Zhejiang University ZJU
Priority date: 2023-05-23
Filing date: 2023-05-23
Publication date: 2023-08-08

Abstract

本发明公开了一种预测突变位点顺式调控基因的方法、设备和介质，属于单细胞测序技术领域。所述方法包括获得表达数量性状位点；在待预测突变位点的周围区域中，若一个表达数量性状位点与所述突变位点具有高连锁不平衡相关性，则将该达数量性状位点关联的基因确定为受所述突变位点顺式调控的基因。本发明的方法扩展了传统的eQTL，提出利用单细胞表达数据来推测cell cluster‑ieQTL，能够克服基因和突变之间的关联被传统大块样本中细胞组成的异质性掩盖的缺陷。另外，本发明通过预测突变位点的顺式调控基因，可以使我们更加准确地进行基因功能分析，促进对疾病基础生物学的了解和药物靶标的鉴定。

Description

一种预测突变位点顺式调控基因的方法、设备和介质

技术领域

本发明涉及单细胞测序技术领域，具体地，涉及一种预测突变位点顺式调控基因的方法、设备和介质。

背景技术

随着全基因组关联研究(genome wide association study，GWAS)的快速发展，与人类复杂疾病和性状具有显著关联的遗传位点的发现日益增加。然而，人们对这些遗传关联背后的生物学机制的理解远远落后于当前已有的数据库资源。在人类复杂疾病和性状的遗传关联位点中，只有小部分位于基因的编码区，而绝大部分则位于基因的非编码区。这意味着决定生物学性状的突变除了影响蛋白质的结构和功能，还可能通过影响基因的转录表达来发挥作用。位于编码区的突变功能相对容易解释，例如，外显子突变效应可以通过遗传密码规则进行评估。但是对于位于非编码区的突变功能的理解则更加困难，迄今为止，还没有一种统一的方法实现对其的精准解释。

表达数量性状位点(expression quantitative trait locus，eQTL)是一类能够影响基因表达水平的遗传位点，它们可以通过在全基因组尺度上对人群中的突变和基因表达水平之间的关联进行测试来鉴定。通过进一步整合eQTL与GWAS的结果，许多研究将基因表达水平作为突变与生物学表型之间的关键分子表型，以揭示GWAS发现的遗传关联位点的生物学功能。迄今为止，绝大多数的eQTL分析都是基于传统大样本转录组测序(bulk RNAsequencing，bulk RNA-seq)进行的。这种方法从百万个裂解细胞中收集RNA并进行测序，所获得的基因表达数据代表的是来自样本中多个不同细胞类型的平均值。过去的许多研究表明，基因调控在不同的细胞类型之间存在非常大异质性。因此，在传统eQTL关联分析中，遗传变异对基因表达的细胞类型特异性影响可能会被掩盖。

细胞类型互作表达数量性状位点(cell type interaction eQTL，cell type-ieQTL)是一种依赖于特定细胞类型的表达数量性状位点。这种新兴的ieQTL算法是对传统eQTL关联分析方法的扩展。该算法首先利用去卷积分析估计组织样本中各种细胞类型的比例，然后通过在线性模型中测试基因型和估计的细胞比例，以推测细胞类型特异性的基因表达调控的遗传变异效应。这种方法在未来的科研和临床应用中表现出极大的潜力，因为它具有极高的通用性和实用性。该方法不需要研究人员额外收集大量特定细胞类型的基因组与转录组数据，只需要使用现有的组织水平遗传数据库(例如GTEx数据集)以及可下载的单细胞转录图谱资源，就能通过计算ieQTL检测细胞类型水平的基因表达遗传调控，而几乎不需要额外的计算资源成本。

然而，由于在组织样本中估计的不同细胞类型的比例通常相互关联(或负相关)，ieQTL除了捕获测试细胞类型的特异性调控效应，还捕获了与之相关的细胞类型的特异性调控效应。这导致了ieQTL计算中一个非常大的局限：难以确定检测到的遗传调控效应的真正细胞类型来源(即因果细胞类型)。此外，由于基因组上的多个位点通常处于连锁不平衡(linkage disequilibrium，LD)，一个遗传位点的ieQTL统计值除了捕获该位点的等位基因调控效应，还捕获了与之连锁的其他遗传位点的等位基因调控效应，因此难以确定真正具有基因表达调控功能的突变(即因果突变)。这些问题限制了ieQTL分析的精准度，阻碍了它在功能基因组学以及临床靶点转化方面的应用。

发明内容

为了解决上述技术问题中的至少一个，本发明采用的技术方案如下：

本发明第一方面提供一种预测突变位点顺式调控基因的方法，包括以下步骤：

S1，获得全基因组水平的表达数量性状位点；

S2，在待预测突变位点的±1Mb区域中，若一个表达数量性状位点与所述突变位点的连锁不平衡相关性R2>0.8，则将该数量性状位点关联的基因确定为受所述突变位点顺式调控的基因。

在本发明的一些实施方案中，所述待预测突变位点是利用以下方法预测的：

基于卷积神经网络提取基因转录起始位点±20kb区域的染色质特征；利用多细胞亚群的单细胞测序数据，生成细胞亚群特异性基因表达谱及相应的细胞亚群信息；用于针对不同细胞亚群，分别利用回归方法，基于所述染色质特征和所述细胞亚群特异性基因表达谱训练出基于染色质特征预测基因表达水平的预测模型；计算机模拟生成一对包含一个突变位点的等位基因，提取所述等位基因的染色质特征并计算所述等位基因的染色质特征的差值；将染色质特征的差值输入到所述不同细胞亚群的预测模型中，针对不同细胞亚群的预测模型基于染色质特征的差值预测所述等位基因在相应细胞亚群中的突变预测效应，若所述等位基因在某个细胞亚群中的突变预测效应的绝对值大于其它细胞亚群，则该等位基因上的突变位点为该细胞亚群的特异性调控突变位点。

在本发明中，“细胞亚群”是比“细胞类型”更下位的一个概念。一种“细胞类型”包括一个或多个“细胞亚群”。

在本发明的一些实施方案中，所述全基因组水平的表达数量性状位点包括标准表达数量性状位点(expression quantitative trait locus，eQTL)和细胞亚群互作表达数量性状位点(cell cluster interaction eQTL，cell cluster-ieQTL)。

其中，标准eQTL是一类能够影响基因表达水平的遗传位点，它们可以通过在全基因组尺度上对人群中的突变和基因表达水平之间的关联进行测试来鉴定。绝大多数的标准eQTL分析都是基于传统大样本转录组测序(bulk RNA sequencing，bulk RNA-seq)进行的，这种方法从百万个裂解细胞中收集RNA并进行测序，所获得的基因表达数据代表的是来自样本中多个不同细胞类型的平均值。在本发明的一些具体实施方案中，所述标准eQTL从GTEx网站上获得。

在本发明的一些实施方案中，所述cell cluster-ieQTL的获取方法如下：

S11，大块样本基于单细胞RNA测序数据的确定细胞亚群：获得大块样本的单细胞RNA测序数据，并获得细胞亚群信息，

S12，大块样本单细胞RNA测序数据的归一化：对大块样本单细胞RNA测序数据进行TMM归一化，保留在20％的样本中TPM不少于0.1以及未归一化的原始序列数不小于6的基因，

S13，样本中协变量的获取：获得样本的技术协变量和/或生物学个体协变量进行合并，去除共线性的协变量得到合并的协变量，所述技术协变量包括但不限于文库构建批次、样本缺血时间和检测到转录本的数量，所述生物学个体协变量包括但不限于年龄、性别、种族；

S14，cell cluster-ieQTL的确定：

(1)在步骤S11中保留的每个细胞亚群中，对于步骤S12中保留的全部基因、以及它们转录起始位点±1Mb区域中突变频率大于0.05的SNP，根据以下线性模型确定cellcluster-ieQTL：

TMM归一化基因表达～β1×SNP基因型:细胞亚群比例+2×SNP基因型+3×

细胞亚群比例+合并的协变量

其中，β1是SNP基因型与细胞亚群比例的交互效应，β2是SNP基因型的效应，β3是细胞亚群比例的效应，对于每一对SNP和基因，在以上公式中利用最小二乘法估计β1的名义p值；

(2)在每一个测试区域中，对正则化基因型相关矩阵进行特征值分解，获取该区域中有效独立SNP数量，并对每个基因最显著的SNP的名义p值进行多重假设检验矫正，生成矫正的p值；

(3)在矫正的p值的基础上，对所有基因的最显著的SNP进行多重假设检验矫正，获取最显著的SNP在全基因组范围的错误发现率，选取错误发现率<0.4的SNP作为cellcluster-ieQTL。

在本发明的一些实施方案中，步骤S11中，进一步包括对混合细胞亚群的大块样本转录组数据进行去卷积分析的步骤：(1)对于细胞数不少于500个细胞的每一个细胞亚群，随机抽取500个细胞，并将每10个抽取的细胞作为一个假细胞；(2)求取假细胞中各基因表达水平的平均值，获得假细胞-基因表达参考矩阵；(3)将所述假细胞-基因表达参考矩阵上传至CIBERSORTx中，创建基因表达特征矩阵；(4)将所述单细胞RNA测序数据上传至CIBERSORTx作为混合文件，并利用S-mode来矫正不同来源的数据之间的批次效应；(5)在绝对模式下使用所述单细胞RNA测序数据进行去卷积分析，利用基因表达特征矩阵来计算每个样本中的细胞亚群的绝对分数，保留在超过1/5的样本中绝对分数大于1的细胞亚群。

在本发明的一些实施方案中，在步骤S13中，利用PEER软件从TMM归一化之后的单细胞RNA测序数据提取能够解释基因表达变化的隐藏因子，在合并时将隐藏因子纳入。

在本发明的一些实施方案中，在步骤S13中，对样本的SNP基因型进行主成分分析，并提取基因型主成分，在合并时将基因型主成分纳入。

在本发明的一些实施方案中，步骤S14中，利用Bonferroni方法对每个基因最显著的SNP的名义p值进行多重假设检验矫正，生成eigenMT矫正的p值；在eigenMT矫正的p值的基础上，利用Benjamini-Hochberg方法对所有基因的最显著的SNP进行多重假设检验矫正。

在本发明的一些实施方案中，对于标准eQTL，首先利用与细胞亚群互作表达数量性状位点相同的方法对p值进行多重假设检验矫正获得在全基因组范围内的错误发现率，选取错误发现率<0.4的标准eQTL进行后续分析。

本发明的第二方面提供一种计算机设备，包括：存储器，用于存储计算机程序；处理器，用于执行所述计算机程序时实现本发明第一方面任一所述的一种预测突变位点顺式调控基因的方法的步骤。

本发明的第三方面一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现本发明第一方面任一所述的一种预测突变位点顺式调控基因的方法的步骤。

本发明的有益效果

相对于现有技术，本发明具有以下有益效果：

本发明的方法扩展了传统的eQTL，提出利用单细胞表达数据来推测cellcluster-ieQTL。能够克服基因和突变之间的关联被传统大块样本中细胞组成的异质性掩盖的缺陷。

本发明通过预测突变位点的顺式调控基因，可以使我们更加准确地进行基因功能分析，促进对疾病基础生物学的了解和药物靶标的鉴定。

附图说明

图1示出了一个cell cluster-ieQTL示意图。图中每个点代表一个大块组织样本。点的颜色代表某个遗传位点的SNP基因型，x轴代表大块组织样本中估计的细胞亚群的比例，y轴代表受该遗传位点顺式调控的基因归一化后的表达水平。

图2示出了GTEx样本中细胞类型的异质性。该图展示了GTEx数据集中20个不同组织类型中各个HCL细胞亚群的相对比例。这些相对比例是通过去卷积分析得到的绝对分数除以每个GTEx样本中所有细胞亚群的绝对分数之和获得的。

图3示出了部分cell cluster-ieQTL。图中展示了本发明实施例2中确定的几个cell cluster-ieQTL。这些cell cluster-ieQTL显示了在(a)动脉、(b)食管和(c)心脏的GTEx样本中，非编码突变对基因表达的调控作用是如何随着去卷积分析估计的细胞亚群比例改变而改变的。x轴表示标准化后的GTEx组织中估计的细胞亚群比例，y轴表示标准化后的基因表达水平。

图4示出了本发明实施例2共定位分析结果。该曼哈顿图展示了仅与cellcluster-ieQTL共定位的几个GWAS遗传位点的示例，说明cell cluster-ieQTL能够揭示与复杂性状和疾病相关的新的基因调控位点。图中从上至下依次为GWAS、cell cluster-ieQTL和标准eQTL。

图5示出了本发明实施例2中各类GWAS中共定位基因的数量统计。该柱状图显示了GWAS与cell cluster-ieQTL以及标准eQTL共定位分析所揭示的共定位基因数量。图中x轴表示114个不同GWAS数据集的缩写，y轴代表共定位基因的数量，柱状图的颜色表示共定位信号只能通过cell cluster-ieQTL(橙色)或标准eQTL(浅蓝色)检测到，或者两种方法均能检测到(粉色)。

图6示出了基因改变与表型变异的剂量-反应曲线。该图表示基因与生物学表型之间的关系模型。当一个致病基因的功能丧失时，通常会导致孟德尔遗传疾病的严重表型。然而，当非编码突变引起该致病基因表达水平的轻微改变时，可能会导致与该基因相关的复杂性疾病的风险增加。该模型揭示了非编码基因调控元件(如cell cluster-ieQTL)在揭示基因与疾病风险之间的关系方面的潜在作用。

图7示出了本发明实施例2中GSEA富集分析图。该GSEA图展示了GWAS和cellcluster-ieQTL的共定位基因在“银标准”中的富集。

图8示出了本发明实施例3中基于不同CNN预测基因表达水平的表现。该图展示了同一细胞亚群特异性模型对990个测试集基因的表达水平预测的表现。y轴为每组模型的预测和实际测量的基因表达水平之间的PCC，x轴为用于比较的几种CNN模型。

图9示出了本发明实施例3中基于不同CNN预测细胞亚群特异性基因表达的表现。该图展示了同一测试集基因在不同细胞亚群特异性模型中预测表达水平的表现。y轴为每组模型的990个测试集基因的预测和实际测量的表达水平之间的PCC，x轴为用于比较的几种CNN模型。

图10示出了本发明实施例3中基于Sei的细胞亚群特异性预测模型的训练损失曲线。该图展示了在利用Sei提取的染色质特征训练基因表达预测模型中，基于训练集和测试集评估的模型表现，y轴是损失函数的值，x轴是训练的轮次。

图11示出了本发明实施例3中基于不同的线性回归方法预测基因表达水平的表现。该图展示了同一细胞亚群特异性模型对990个测试集基因的表达水平预测的表现。y轴为每组模型的预测和实际测量的基因表达水平之间的PCC，x轴为用于比较的几种线性回归方法。

图12示出了本发明实施例3中基于预测模型的基因表达预测。该密度图显示了在动脉组织细胞亚群(a)成纤维细胞、(b，c)内皮细胞和(d)平滑肌细胞中，预测模型预测的和实际观测的基因表达之间的一致性。x轴表示对数转化的基因表达预测值，y轴表示基于单细胞表达数据得到的对数转化的基因表达水平。

图13示出了本发明实施例3中预测模型在20个组织中对基因表达的预测表现。该箱线图显示了不同组织的细胞亚群特异性模型的基因表达预测表现。x轴为不同组织的预测模型分组，y轴为预测和实际观测到的基因表达水平的PCC相关性。

图14示出了本发明实施例3中预测模型对基因表达水平高低的预测。该密度图显示了在动脉组织的几种细胞亚群中，当只考虑观察到表达的测试集基因时，预测模型对基因表达水平的预测表现。x轴为对数转化的基因表达预测值，y轴表示基于单细胞表达数据得到的对数转化的基因表达水平。

图15示出了本发明实施例3中预测模型对细胞亚群特异性基因表达的预测表现。该树状图基于每个模型中细胞亚群特异性基因表达的预测值与所有细胞亚群中细胞亚群特异性基因表达的观察值之间的PCC，描述了357个细胞亚群的关系。树上方不同的颜色表示每个细胞亚群对应的组织信息，中间不同的颜色表示对应的细胞谱系信息，右边显示了每个细胞亚群模型中预测和观察的基因表达水平的相关性

图16示出了本发明实施例3中基于预测模型的eQTL统计值预测。该图显示了在外周血组织中，预测模型预测的突变效应与eQTL统计值z-scores之间的相关性。每个点代表一个独立的LD block，x轴显示了该基因组区域中最大的对数转化的突变预测效应绝对值，y轴显示了该基因组区域中最大的对数转化的eQTL z-scores绝对值。

图17示出了本发明实施例3中不同突变组别在cell cluster-ieQTL的富集。该箱型图展示了基于预测模型预测值分类的两组突变位点在相应组织中cell cluster-ieQTL的富集。突变位点组2代表每个组织中预测的细胞亚群特异性调控位点，这些突变的细胞亚群模型的预测绝对值显著大于bulk水平模型的预测绝对值。突变位点组1代表每个组织中突变位点组2之外的其它突变位点。每个点代表一个测试组织。

图18示出了本发明实施例3中预测模型预测精细映射突变位点的表现。图为基于预测模型预测的逻辑分类器判断一个突变位点是否为精细映射突变位点的ROC曲线。左图和右图分别为模型在距离(a)TSS 20kb以内和(b)TSS 20kb以外的预测表现。

图19示出了本发明实施例4中整合分析的框架流程图。(a)通过使用CNN模型Beluga，整合每个基因TSS周围的染色质特征；(b)根据细胞图谱，为每个细胞亚群分别训练细胞亚群特异性的突变预测模型；(c)根据学生T检验对所有预测的突变进行排序，以获得对每个细胞亚群中特异性的突变位点；(d)根据人群中LD相关性对eQTL和突变位点预测结果进行整合，确定受预测到的突变位点顺式调控的候选基因。

图20示出了本发明实施例4整合分析示意图。根据LD相关性对预测的突变位点和eQTL进行整合分析。

图21示出了本发明实施例4中关联分析阈值测试。(a)该图展示了基于ChromHMM软件对动脉组织中的六种不同染色质修饰数据进行整合的结果。每行代表具有相似染色质状态的一组基因组区域；(b)不同阈值下的eQTL位点在相应组织中cCREs的富集结果。x轴表示不同eQTL组别，y轴为每组eQTL相比于匹配的对照组正好位于cCREs的OR。

图22示出了本发明实施例4中细胞亚群特异性基因调控遗传变异情况。该图统计了每种细胞亚群的特异性功能调控突变以及受其顺式调控的关联基因的数量。

图23示出了本发明实施例4中得到的细胞亚群特异性基因调控遗传变异情况与scATAC-seq数据集的比较分析。该热图展示了每种细胞亚群中推断的特异性调控突变(y轴)在scATAC-seq数据集中细胞亚群特异性开放染色质区域(x轴)的标准化后的富集分数。该图反映了推断的细胞亚群特异性调控变异和实际观察到的细胞亚群特异性调控区域之间的良好一致性。x轴和y轴的标签颜色代表了细胞亚群所属的谱系。

图24示出了本发明实施例4中cell cluster-ieQTL与外部数据集的比较分析。该热图展示了每种细胞亚群中cell cluster-ieQTL位点(y轴)在scATAC-seq数据集中细胞亚群特异性开放染色质区域(x轴)的标准化后的富集分数。相比之下，cell cluster-ieQTL位点与实际观察到的细胞亚群特异性调控区域之间的一致性较差。x轴和y轴的标签颜色代表了细胞亚群所属的谱系。

具体实施方式

除非另有说明、从上下文暗示或属于现有技术的惯例，否则本申请中所有的份数和百分比都基于重量，且所用的测试和表征方法都是与本申请的提交日期同步的。在适用的情况下，本申请中涉及的任何专利、专利申请或公开的内容全部结合于此作为参考，且其等价的同族专利也引入作为参考，特别这些文献所披露的关于本领域中的相关术语的定义。如果现有技术中披露的具体术语的定义与本申请中提供的任何定义不一致，则以本申请中提供的术语定义为准。

为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。

实施例

以下例子在此用于示范本发明的优选实施方案。本领域内的技术人员会明白，下述例子中披露的技术代表发明人发现的可以用于实施本发明的技术，因此可以视为实施本发明的优选方案。但是本领域内的技术人员根据本说明书应该明白，这里所公开的特定实施例可以做很多修改，仍然能得到相同的或者类似的结果，而非背离本发明的精神或范围。

除非另有定义，所有在此使用的技术和科学的术语，和本发明所属领域内的技术人员所通常理解的意思相同，在此公开引用及他们引用的材料都将以引用的方式被并入。

那些本领域内的技术人员将意识到或者通过常规试验就能了解许多这里所描述的发明的特定实施方案的许多等同技术。这些等同将被包含在权利要求书中。

下述实施例中的实验方法，如无特殊说明，均为常规方法。下述实施例中所用的仪器设备，如无特殊说明，均为实验室常规仪器设备；下述实施例中所用的试验材料，如无特殊说明，均为自常规生化试剂商店购买得到的。

实施例1单细胞测序及数据处理

1.单细胞测序

获得与GTEx匹配的20个成年人体组织类型，共39例样本，详细信息如表1所示。

表1本实施中分析的组织及相关信息

利用微孔板测序技术(Microwell-seq)(Han,Xiaoping et al.Mapping theMouse Cell Atlas by Microwell-Seq,Cell,Volume 172,Issue 5,1091-1107.e17)对各样本进行高通量单细胞转录组测序(scRNA-seq)。对Microwell-seq数据进行质控，去除RNA背景污染并过滤潜在双细胞双细胞后，得到高质量的单细胞数据集。

2.整合高质量的单细胞数据集

在对scRNA-seq数据进行背景污染矫正和潜在双细胞过滤后，发明人进一步对处理后的数据进行无监督聚类，具体步骤如下：

(1)对于检测基因数大于1000的高质量单细胞数据，将相同组织的批次的基因表达矩阵利用merge函数合并在一起。使用R包Seurat中的函数CreateSeuratObject创建Seurat对象，在该函数中，将最小细胞设置为3(min.cells＝3)，最小基因设置为200(min.features＝200)。

(2)在创建的Seurat对象meta.data中，输入HCL原文注释的详细细胞类型信息。

(3)利用函数NormalizeData(Seurat包)对Seurat对象中的基因表达矩阵进行归一化。其中，将归一化方法设置为LogNormalize，缩放因子设置为10000。该步骤将每个基因的计数除以该细胞总的基因计数，再乘以指定的缩放因子10000，最后使用log1p进行自然对数转换。

(4)使用函数FindVariableFeatures(Seurat包)寻找每个Seurat对象中的高变基因。在该函数中，将方法设置为vst，将高变基因数量设置为2000。该步骤首先使用局部多项式拟合log转化的基因表达方差与均值之间的关系，然后使用观测均值和由拟合线给出的期望方差对特征值进行标准化，最后基于标准化后的特征值计算方差。

3.数据降维处理和亚群定义

在整合每个组织的单细胞数据集后，发明人进一步对其进行了降维和分群分析，具体步骤如下：

(1)使用ScaleData函数(Seurat包)对整合的Seurat对象中的基因表达进行中心化，使用RunPCA(Seurat包)对高变基因进行主成分分析(Principal component analysis，PCA)，并参考ElbowPlot的曲线拐点选择后续使用的PCA维度。

(2)在RunTSNE(Seurat包)函数中，通过选取的PCA维度使用fit-tsne进行降维。在该函数中，将dims设置为选取的PCA数。

(3)使用FindNeighbors和FindClusters函数(Seurat包)，基于细胞间的距离进行分群。将FindNeighbors中k.param设置为10，维度(dims)设置为选择的PCA维度；将FindClusters中迭代次数设置为50，分辨率设置为1。

(4)使用DimPlot函数(Seurat包)来对每个组织的分群情况进行可视化。在该函数中，将reduction设置为“FItSNE”。

(5)通过FindAllMarkers函数(Seurat包)，利用Wilcoxon秩和检验来计算每个细胞亚群的特征高表达基因(marker)。

(6)将HCL原文中的注释投射至高质量数据的细胞亚群中。

最后，发明人对细胞亚群进行归类，将其划分为不同的主要细胞亚群以及细胞谱系。在本实施例中，一共获得357个细胞亚群，其中每个组织中包含的细胞亚群数量详见表2。

表2单个组织中包含的细胞亚群数量

发明人将这357个细胞亚群进一步归类为44个主要细胞类型和9个不同的谱系，具体信息见表3。

表3高质量细胞图谱的主要细胞类型

实施例2确定细胞亚群互作eQTL的方法

目前，eQTL通常是基于大块样本的测序数据，并利用传统的遗传关联分析方法推测得到的。这些被推测的eQTL通常被称为标准eQTL。标准eQTL可以从GTEx网站上获得。在标准eQTL中，细胞亚群特异性的基因调控作用可能被组织样本中的细胞异质性所掩盖，尤其是对于那些在样本中比例较小的细胞亚群来说。

为了实现在细胞亚群水平上对基因表达遗传调控的探索，发明人扩展了传统的eQTL关联分析方法，提出利用单细胞表达数据来推测细胞亚群互作eQTL(cell clusterinteraction eQTL，cell cluster-ieQTL)。发明人将cell cluster-ieQTL定义为依赖于细胞图谱中的细胞亚群与基因表达水平产生关联的表达数量性状位点。如图1所示，在单核苷酸多态性(single nucleotide polymorphism，SNP)遗传位点上具有更高T等位基因剂量的样本中，基因表达水平与样本中某细胞亚群的比例呈正相关。此外，在含有更高该细胞亚群的组织样本中，SNP基因型与基因表达之间具有更强相关性。这些观察表明，该SNP位点是一个细胞亚群特异性的基因表达调控位点，因此被判定为cell cluster-ieQTL。

本实施例中，cell cluster-ieQTL确定方法的主要步骤分为大块样本的去卷积分析、大块样本表达矩阵的归一化、样本中协变量的获取和交互效应的检测。

1.大块样本的去卷积分析

(1)对于scRNA-seq数据中每一个群，随机抽取500个细胞(小于500个细胞的群则不抽样)，并将抽取的细胞每10个细胞作为一个假细胞(pseudo-cell)。

(2)求取假细胞中各基因表达水平的平均值，获得假细胞-基因表达参考矩阵。

(3)在软件CIBERSORTx中，分析模块选择“Create Signature Matrix”，分析模式选择“Custom”，输入数据类型选择“scRNA-seq”，然后将上述假细胞-基因表达参考矩阵上传至自定义输入中的“Single cell reference matrix file”处，作为单细胞参考矩阵文件。此外在附加选项“Single Cell Input Options”中，将基因表达阈值“Min.Expression”设置为0。其他选项默认，以此创建基因表达特征矩阵。

(4)将大块样本RNA-seq测序数据上传至CIBERSORTx作为混合文件，并利用S-mode(single cell mode)来矫正不同来源的数据之间的批次效应。

(5)在绝对模式下对大块样本进行去卷积分析，利用基因表达特征矩阵来计算每个样本中的细胞亚群的绝对分数。绝对模式下生成的绝对分数将相对细胞分数缩放成反映混合物中每种细胞亚群的绝对比例的分数，虽然绝对分数没有直接表示为细胞占比，但它可以用于在不同细胞亚群之间直接比较所估计的细胞比例大小。保留在超过1/5的样本中绝对分数大于1的细胞亚群用于下游分析。

同时，为避免异常值的影响，对计算的绝对分数进行逆正态变换(inverselynormally transformed，INT)。

2.大块样本表达矩阵的归一化

对大块样本RNA-seq测序数据进一步地进行TMM(Trimmed Mean of M-values)归一化，保留在20％的样本中TPM(Transcripts per million)不少于0.1以及未归一化的原始序列数不小于6的基因。

3.样本中协变量的获取

利用PEER软件从TMM归一化之后的表达矩阵中提取能够解释基因表达变化的隐藏因子。PEER推断的隐藏因子(PEER factor)和人群表达数据中存在的许多混杂因素显著相关，包括技术的协变量(如文库构建批次、样本缺血时间、检测到转录本的数量……)和生物学个体的协变量(如年龄、性别、种族……)。因此，在cell cluster-ieQTL确定过程中，将PEER factor作为协变量可以减少结果的假阳性，提高cell cluster-ieQTL位点的检测能力。PEER因子的数量选择取决于cell cluster-ieQTL关联测试中所用到的样本数量。

利用PLINK2软件，对测试样本的SNP基因型进行主成分分析，并从突变矩阵中提取前五个主成分(principal components，PCs)。

将提取的PEER factor、基因型主成分、文库构建批次、样本缺血时间、检测到转录本的数量、测序平台、年龄、性别、种族等其它已知协变量合并在一起，然后检测并丢弃共线性的协变量。

4.交互效应的检测

在细胞图谱中的每个细胞亚群中，对于TMM归一化中保留的全部基因、以及它们转录起始位点(TSS)周围1Mb区域中突变频率(MAF)大于0.05的突变(SNP)，在tensorqtl软件中根据以下线性模型进行cell cluster-ieQTL的计算：

细胞亚群比例+合并的协变量

其中，β1是SNP基因型与细胞亚群比例的交互效应，β2是SNP基因型的效应，β3是细胞亚群比例的效应。对于每一对SNP和基因，在以上公式中利用最小二乘法估计β1的名义p值(nominal p-values)。

在每一个ieQTL测试区域(cis-window)中，利用软件eigenMT对正则化基因型相关矩阵进行特征值分解(SNP基因型从GTEx中获得)，获取该区域中有效独立的突变数量。并利用Bonferroni方法对每个基因最显著的ieQTL关联(top ieQTL)的名义p值进行多重假设检验矫正，生成eigenMT矫正的p值(eigenMT-adjusted p-value，pval_emt)。

在eigenMT矫正的p值的基础上，利用Benjamini-Hochberg(BH)方法对所有测试基因的top ieQTL进行多重假设检验矫正，获取cell cluster-ieQTL在全基因组范围的错误发现率(false discovery rate，FDR)。

结果发现，在同一种组织类型的GTEx样本中细胞亚群的比例存在显著差异(如图2)。本实施例确定出大量细胞亚群特异性的遗传调控位点，即cell cluster-ieQTL，这些调控位点对基因表达的影响随着样本中估计的细胞亚群比例变化而变化(如图3)。对于所有细胞亚群，该分析共推测出了214403个cell cluster-ieQTL(FDR＜0.4)。

5.cell cluster-ieQTL推测结果合理性检验

为了评估推测结果的合理性，发明人收集了每个组织中FDR＜0.4的所有cellcluster-ieQTL，并将这些cell cluster-ieQTL位点分别与GTEx数据库中相应组织的标准eQTL以及细胞亚群依赖性eQTL进行比较。

结果发现，与标准eQTL相比，cell cluster-ieQTL与细胞亚群依赖性eQTL之间具有更多的重叠，20个组织中OR中位数为3.769(双侧Fisher精确检验)。该结果符合预期，因为标准eQTL分析可以检测到在各种细胞亚群中都表现出相似基因表达调控的保守性位点，而cell cluster-ieQTL分析则侧重于检测在组织样本中的特定细胞亚群中表现出基因表达调控的位点。

然后，发明人分别对cell cluster-ieQTL和标准eQTL与GWAS数据集进行了全面的贝叶斯共定位分析。在多个组织中的研究中发现，使用cell cluster-ieQTL与GWAS进行的共定位分析能够揭示那些在使用标准eQTL分析时无法检测到的共定位信号(如图4)。在这些共定位区域中，cell cluster-ieQTL与GWAS被认为具有相同的因果突变。同时，cellcluster-ieQTL关联的基因(即ieGenes)被认为是共定位基因(colocalized genes)，其表达水平也被认为是遗传变异与GWAS性状中间的关键分子表型。在所有的GWAS数据集中，cell cluster-ieQTL揭示了大量新颖的共定位基因，这些基因可能与复杂性状和疾病具有因果关联(PP.H4＞0.85，如图5)。

为了进一步验证分析结果的合理性，发明人使用了GWAS性状因果基因的“银标准”数据集，将其与cell cluster-ieQTL揭示的共定位基因进行比较。需要注意的是，“银标准”数据集是基于孟德尔单基因疾病的OMIM数据库建立的，而GWAS性状是由多个基因共同作用导致。发明人依据孟德尔单基因疾病基因收集GWAS性状因果基因的理由是：在基因功能完全丧失将导致孟德尔疾病严重表型的情况下，由非编码调控变异引起的基因表达水平轻微改变则可能导致较为温和的相关生物学表型的变化，进而引起GWAS性状发生的风险(如图6)。发明人发现，具有更高PP.H4的GWAS性状-基因对在“银标准”数据集中显著富集(p＜1e-4，如图7)。总的来说，这些结果显示了本实施例推测的cell cluster-ieQTL结果在生物学上是合理的。

实施例3从头预测细胞亚群水平的突变

本实施例提供基于基于实施例1中scRNA-seq数据构建细胞亚群特异性深度学习模型的流程图。主要步骤分为染色质状态bulk水平的预测、细胞亚群特异性基因表达谱的获取、细胞亚群特异性线性回归模型的训练和基因表达预测表现的评估。

1.染色质状态bulk水平的预测

利用CNN模型提取目前已知的24339个蛋白质编码、核糖体RNA(ribosomal RNA，rRNA)和长链非编码RNA(long noncoding RNA,lncRNA)的DNA序列(基因)的染色质特征(包括组蛋白标记、转录因子结合和DNase超敏反应谱等表观数据集)。

本实施例中，发明人选择三种CNN模型进行染色质特征提取：

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)个染色质特征。

2.细胞亚群特异性基因表达谱的获取

1)对于24339个基因，单细胞测序数据中基因表达原始计数矩阵通过log2(CPM_i,j/10+1)来进行归一化，其中CPM_i,j表示细胞j中基因i的每百万计数。缺失值用0来进行填充。

2)使用MAGIC算法对以上得到的归一化表达矩阵进行填补，来弥补由于单细胞测序导致的矩阵稀疏性。需要说明的是，单细胞RNA-seq数据通常包含许多由于原始RNA扩增失败而导致的缺失值或dropout，因此低表达的基因往往捕获不到，如果不进行Imputation，可能导致后续模型学习及评估不够充分。

3)为了进一步增加测得的基因数量和UMI数量，发明人对scRNA-seq数据中同一个群的所有细胞的基因表达量取平均值，生成细胞亚群特异性的基因转录谱。

3.细胞亚群特异性线性回归模型的训练

为了捕获细胞图谱中不同细胞亚群的表观特征，发明人基于CNN模型提取的每个基因TSS周围的染色质状态，利用线性回归模型来拟合细胞亚群特异性的基因表达水平。在模型训练中，将位于8号染色体上的990个基因划分为测试集，将其余基因划分为训练集，利用训练集进行模型训练。

4.基因表达预测表现的评估

1)对于每个细胞亚群特异性模型，通过计算在990个测试集基因中的预测表达水平和实际测量表达水平的皮尔森相关性(Pearson’s correlation coefficient，PCC)，来检测模型在捕获不同基因的表达水平差异方面的表现。

2)对于每个测试集基因，通过计算在不同scRNA-seq数据细胞亚群中的预测表达水平和实际测量表达水平的PCC值，来检测模型在捕获细胞亚群特异性基因表达水平方面的表现。

5.基于不同CNN和线性回归方法的模型性能评估

在构建细胞亚群特异性的突变预测模型中，为了选择适合的CNN和线性回归模型，发明人针对不同的CNN模型和不同的线性回归模型，分别进行了基因表达预测，从而评估不同的模型选择对基因表达预测的影响。

(1)基于不同CNN的模型性能评估

然后，发明人按照前述步骤和方法，根据不同CNN提取的染色质特征，分别拟合人类细胞图谱(Human Cell Landscape，HCL)中不同细胞亚群的基因表达水平，由此构建了三组细胞亚群特异性的预测模型，均利用XGBoost线性回归模型来拟合细胞亚群特异性的基因表达水平，在超参数的选择中，将booster设置为线性；将训练最大迭代次数设置为40，并使用迭代过程的早期停止以避免过拟合(通过设置early_stopping_rounds＝10)；将用于更新模型的步长收缩设置为0.5；将模型的损失函数设置为线性回归；将权重的L2正则化项增加为100；将其余的超参数设为默认值。最后对不同组别的模型预测表现进行评估和比较。

结果显示，在捕获不同基因表达水平的差异方面，基于Beluga的模型表现优于基于DeepSEA和Sei的模型(如图8)。在捕获不同细胞亚群基因表达水平的差异方面，基于Beluga训练的模型同样也取得了最好的表现(如图9)。虽然基于Sei的模型使用了最多的染色质特征来预测基因表达水平，但其表现最差。分析发现，这是由于Sei提取的染色质特征中存在大量与基因表达无关的噪音信息，加上训练样本的有限，模型出现了过拟合(如图10)。

(2)基于不同线性回归方法的模型性能评估

由于预测基因表达的染色质因子特征数量非常大(每个基因400400个特征)，并且用于训练的基因样本数量相对有限(24339个基因)，因此使用线性回归模型拟合基因表达时可能会遇到多重共线性和过拟合的潜在问题。为了选择适合的线性回归方法，发明人利用Beluga提取染色质特征，分别利用普通最小二乘法线性回归(ordinary least-squareslinear regression)、套索回归(lasso regression)、岭回归(ridge regression)以及XGBoost梯度提升方法分别拟合HCL中不同细胞亚群的基因表达水平，从而构建得到四组细胞亚群特异性的预测模型。发明人对不同组别的模型预测表现进行评估和比较。

结果显示，在几种不同的线性回归模型中，基于XGBoost的预测模型表现最好，其次是基于套索回归和岭回归的模型，而基于普通最小二乘法线性回归的预测表现最差(如图11)。

此外，发明人还对基于不同的线性回归方法训练细胞亚群特异性预测模型所需要的时间进行了比较。结果如表4所示，可以看出，XGBoost在训练细胞亚群特异性预测模型的速度方面显著优于套索回归和岭回归。

表4利用不同的回归方法训练基因表达预测模型的时间

因此，XGBoost能够更好地处理该线性回归问题，并更快地为scRNA-seq数据批量建立细胞亚群特异性的预测模型。

综上所述，本实施例利用Beluga模型提取染色质特征，进一步利用XGBoost回归拟合基因表达水平，从而建立了细胞亚群特异性的突变预测模型。

5.突变预测模型性能评估

发明人分别构建了来自20个组织的357个细胞亚群的特异性突变预测模型。为了评估模型的性能，发明人首先计算了每个模型预测的测试集基因表达水平与实际测量的表达水平之间的相关性。

结果表明，这些细胞亚群特异性模型能够有效地预测相应细胞亚群中的基因表达水平(p＜2.2e-16，图12)。总体来说，在357个细胞亚群的模型中，测试集基因的预测和实际观察到的基因表达水平之间的PCC中位数达到了0.763。而在肾脏、胃和横结肠的某些细胞亚群模型中，PCC甚至超过了0.80(图13)。为了进一步证明该模型可以准确预测基因表达水平的变化，而不仅仅是二元的转录状态(表达或不表达)，发明人只保留在每个细胞亚群中观察到表达的测试集基因，并重新计算预测和实际观察到的表达水平之间的一致性。结果表明，预测模型仍然显示出对基因表达良好的预测性能。在357个细胞亚群模型中，PCC中位数达到了0.63(图14)。

此外，发明人还评估了预测模型在预测细胞亚群特异性基因表达方面的能力。具体而言，发明计算了每个细胞亚群中预测的基因表达水平与其他所有细胞亚群中观察到的基因表达水平相对差异的对数值，并将其作为细胞亚群特异性基因表达水平的衡量指标。随后，发明人计算了每个模型中细胞亚群特异性基因表达的预测值与所有细胞亚群中细胞亚群特异性基因表达的观察值之间的PCC，并以此对细胞亚群进行了无监督聚类。结果表明，预测模型能够准确地预测不同细胞亚群的谱系特异性基因表达，同一谱系内的基因表达预测值具有更高的一致性，而不同谱系之间的基因表达预测值差异较大(图15)。

6.具有特异性调控功能的突变位点的预测

首先，通过计算机模拟DNA碱基替换，得到包括假定突变位点的基因序列作为备选等位基因序列。其次，计算参考基因序列(野生型基因序列)和备选等位基因序列在Beluga输出的染色质特征的差值。接着，将该差值输入训练好的细胞亚群特异性的XGBoost模型，以预测突变在特定细胞亚群中对基因表达的突变预测效应(干扰效应)。然后，根据单侧学生T检验(Student's t test)，选择top 2000并且突变预测效应绝对值大于0.5的位点，测试特定细胞亚群中突变预测效应的绝对值是否大于其它细胞亚群，并据此来确定每种细胞亚群的特异性调控突变位点。

发明人还先将等位基因序列的染色质特征输入到输入训练好的细胞亚群特异性的XGBoost模型，以预测等位基因序列在特定细胞亚群中的表达水平。同样，根据单侧学生T检验(Student's t test)，测试特定细胞亚群中等位基因的表达水平与参考基因序列在相应细胞亚群中的表达水平的差值的绝对值是否大于其它细胞亚群，并据此来确定每种细胞亚群的特异性调控突变位点。

7.突变预测结果的合理性评估

为了检测突变效应预测结果的合理性，发明人首先测试了预测模型是否可以仅基于DNA序列重现基于人群数据的关联分析结果。发明人将来自每个组织中所有细胞亚群模型的预测结果汇总在一起，并与标准eQTL进行比较。结果显示，在同一个连锁不平衡区域(LD block，r2≥0.6)中，eQTL统计值z-scores(效应大小β/标准误)与突变预测效应的最高绝对值表现出显著的相关性(p＜2.2e-16，图16)。在20个测试的组织中，两者的斯皮尔曼相关系数(Spearman’s Correlation Coefficient，SCC)的中位数为0.3。

接着，发明人测试了预测模型能否仅基于DNA序列识别cell cluster-ieQTL位点。需要注意的是，在一个研究的组织中，本领域技术人员通常无法确定cell cluster-ieQTL的基因调控效应实际上来源于何种细胞亚群。这是因为不同细胞亚群在测试样本中估计的比例是线性相关的，cell cluster-ieQTL的调控效应既可能来源于测试SNP基因型交互效应的细胞亚群，也可能来源于与之线性相关的其他细胞亚群。尽管如此，发明人发现，在同一组织中，基于所有细胞亚群的预测模型推测的遗传调控位点在所有细胞亚群的cellcluster-ieQTL中显著富集(在18/20组织中，Bonferroni矫正的p值小于0.05)。

然后，发明人利用大块组织样本的基因表达数据为20个组织分别训练了bulk水平模型，并对突变效应进行了重新预测。对于每个组织，发明人比较了细胞亚群水平的突变预测和bulk水平的突变预测结果。发明人按照细胞亚群模型的预测效应绝对值是否大于bulk水平模型的预测效应绝对值的1.96个标准差，将突变位点分为两组。结果发现，那些在细胞亚群模型中具有明显更高预测绝对值的突变位点，同样在cell cluster-ieQTL中也具有更为显著的富集(图17)。这些结果表明，在揭示细胞亚群特异性遗传调控方面，预测模型的预测与cell cluster-ieQTL具有一致性。

此外，发明人还测试了预测模型能否直接指出基因调控的因果突变位点。由于人群中等位基因LD的存在，与基因具有最显著关联的突变位点并不一定是真正导致该基因表达水平改变的突变。eQTL因果突变位点的鉴定可以更加准确地进行基因功能分析，促进对疾病基础生物学的了解和药物靶标的鉴定。预测模型的突变预测基于DNA序列碱基替换实现，因此理论上其结果不受LD的影响。发明人使用组织中每个变异体的最高预测绝对值，与GTEx组织的精细映射eQTL数据集进行比较。结果表明，预测模型可以在相互具有较高LD关联(r2≥0.6)的突变位点中，实现区分精细映射的阳性和阴性突变，其平均AUROC为0.780(距离TSS 20kb以内AUROC＝0.782，20kb以外AUROC＝0.766，图18)。因此，预测模型可以识别因果突变方面，在近端和远端调控区域中都表现出良好的性能。

为了进一步检验预测模型的突变预测能力，发明人在人永生化红细胞系(HUDEP-2)中训练了预测模型，并对HBG基因TSS周围1kb基因组区域的突变进行了逐一地预测。结果表明，致病突变的预测效应大小是与它们匹配的随机选择的突变集合的3.25倍(95％ CI＝1.62～6.16)。这些结果表明预测模型能够在一定程度上识别非编码致病突变，从而为其在单碱基水平的因果突变预测提供了有力的证据。

实施例4预测受突变位点顺式调控的基因

由于基因组调控元件可以通过三维空间结构(例如DNA环状结构)的协同作用来影响基因的转录表达，仅仅基于DNA序列可能无法充分估计突变对于特定基因表达水平的影响，并且也缺乏进一步的遗传学证据支持。因此，发明人通过整合实施例2得到的cellcluster-ieQTL和实施例3预测的突变位点，确定受突变位点顺式调控的基因。如图19所示，整合分析(荟萃分析)的具体步骤如下：

(1)获取标准eQTL，并利用实施例2的方法步骤获得cell cluster-ieQTL。为了获得标准eQTL在全基因组范围的FDR，发明人采用了与cell cluster-ieQTL多重假设检验矫正方法相同的方法。具体而言，使用BH方法对所有eQTL的基因水平的显著性进行多重假设检验矫正，以获得在全基因组范围内的FDR，选取全基因组显著性FDR＜0.4作为标准eQTL和cell cluster-ieQTL的阈值进行后续分析。

(2)基于标准eQTL和cell cluster-ieQTL计算中使用的SNP基因型矩阵，使用PLINK来估计遗传位点的连锁不平衡(LD)相关性。

(3)在每个±1Mb的测试区域中，当一种细胞亚群测试出的eQTL或cell cluster-ieQTL与同一种细胞亚群预测出的特异性调控突变位点呈高连锁不平衡(LD R²>0.8)，则将该QTL或cell cluster-ieQTL关联的基因推测为受该预测的突变位点顺式调控的基因(如图20所示)。

在步骤(1)中，为了选择合适的标准eQTL和cell cluster-ieQTL的FDR阈值，以便进行后续的整合分析，发明人对动脉组织进行了基准测试，并将其与包括启动子和增强子在内的候选顺式调节元件(cCREs)进行了比较。发明人整合了动脉组织中不同染色质组蛋白修饰的染色质免疫沉淀后测序(chromatin immunoprecipitation sequencing，ChIP-seq)数据，并将染色质划分为E1-E10 10个不同的染色质状态(如图21a)。发明人最终确定将基因组区域E1、E2和E4作为cCREs，其中E1和E2为增强子，E4为启动子。

接着，发明人将动脉组织的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作为阈值更为适合。该阈值能够在真阳性和假阳性数量之间提供了一个合理的平衡，在保证准确度的同时，也能够推测出更多基因调控遗传突变位点。

在本实施例中，发明人利用实施例3的方法预测13182个细胞亚群特异性功能突变位点，利用本实施例的方法预测到6181基因的表达水平受到细胞亚群特异性功能突变顺式调控的影响(图22)。

为了测试本实施例的方法是否正确地识别了各个细胞亚群中的特异性遗传调控效应，发明人将突变位点-基因关联关系与一个外部的scATAC-seq数据集(Zhang,K.etal.A single-cell atlas of chromatin accessibility in the human genome.Cell184,5985-6001e5919)进行比较，并计算了各个细胞亚群中推断的调控位点在实际检测的细胞亚群特异性开放染色质区域的富集程度。结果显示，相同谱系的细胞亚群具有类似的富集情况，而基于富集结果的无监督层次聚类能够反映这些细胞亚群之间众所周知的谱系关系(图23)。同时，还在同一谱系的细胞亚群中，基因调控位点在差异开放染色质峰中表现出特异性富集。因此，该富集分析表明，本实施例推测的细胞亚群特异性调控位点与实际观察到的细胞亚群特异性调控区域之间具有一致性。

此外，由于cell cluster-ieQTL也是细胞亚群依赖的遗传调控位点，发明人使用类似的方法计算了这些位点在scATAC-seq测量的细胞亚群特异性开放染色质区域的富集程度，并进行了无监督聚类分析。相比之下，该富集结果显示出较低的细胞谱系特异性。例如，一种细胞亚群的cell cluster-ieQTL位点除了在相同谱系的细胞亚群的差异开放染色质峰中富集外，也在其他谱系的细胞亚群中表现出明显的富集(图24)。这是因为，在一种细胞亚群中检测到的cell cluster-ieQTL可能实际上来源于其他细胞亚群中发生的基因调控。因此，相比单纯的cell cluster-ieQTL分析，本实施例的整合分析作为一个集成框架能够更好地识别感兴趣细胞亚群的特异性基因调控位点。

在本发明中，可以利用计算机模块实现上述各步骤的功能，各个模块集合或组装和配合成系统，进一步利用存储器以处理器可执行的计算机程序形式储存系统。还可利用计算机可读存储介质存储上述系统。

在本发明提及的所有文献都在本申请中引用作为参考，就如同每一篇文献被单独引用作为参考那样。此外应理解，在阅读了本发明的上述讲授内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。

Claims

1.一种预测突变位点顺式调控基因的方法，其特征在于，包括以下步骤：

S1，获得全基因组水平的表达数量性状位点；

S2，在待预测突变位点的±1Mb区域中，若一个表达数量性状位点与所述突变位点的连锁不平衡相关性R²>0.8，则将该数量性状位点关联的基因确定为受所述突变位点顺式调控的基因。

2.根据权利要求1所述的一种预测突变位点顺式调控基因的方法，其特征在于，所述全基因组水平的表达数量性状位点包括标准表达数量性状位点和细胞亚群互作表达数量性状位点。

3.根据权利要求2所述的一种预测突变位点顺式调控基因的方法，其特征在于，所述细胞亚群互作表达数量位点的获取方法如下：

S11，大块样本基于单细胞RNA测序数据的确定细胞亚群：

获得大块样本的单细胞RNA测序数据，并获得细胞亚群信息，

S12，大块样本单细胞RNA测序数据的归一化：

对大块样本单细胞RNA测序数据进行TMM归一化，保留在20％的样本中TPM不少于0.1以及未归一化的原始序列数不小于6的基因，

S13，样本中协变量的获取：

获得样本的技术协变量和/或生物学个体协变量进行合并，去除共线性的协变量得到合并的协变量，所述技术协变量包括文库构建批次、样本缺血时间和检测到转录本的数量，所述生物学个体协变量包括年龄、性别、种族；

S14，细胞亚群互作表达数量位点的确定：

(1)在步骤S11中保留的每个细胞亚群中，对于步骤S12中保留的全部基因、以及它们转录起始位点±1Mb区域中突变频率大于0.05的SNP，根据以下线性模型确定细胞亚群互作表达数量位点：

TMM归一化基因表达～β1×SNP基因型：细胞亚群比例+2×SNP基因型+3×

细胞亚群比例+合并的协变量

(3)在矫正的p值的基础上，对所有基因的最显著的SNP进行多重假设检验矫正，获取最显著的SNP在全基因组范围的错误发现率，选取错误发现率<0.4的SNP作为细胞亚群互作表达数量位点。

4.根据权利要求3所述的一种预测突变位点顺式调控基因的方法，其特征在于，步骤S11中，进一步包括对混合细胞亚群的大块样本转录组数据进行去卷积分析的步骤：

(1)对于细胞数不少于500个细胞的每一个细胞亚群，随机抽取500个细胞，并将每10个抽取的细胞作为一个假细胞；

(2)求取假细胞中各基因表达水平的平均值，获得假细胞-基因表达参考矩阵；

(3)将所述假细胞-基因表达参考矩阵上传至CIBERSORTx中，创建基因表达特征矩阵；

(4)将所述单细胞RNA测序数据上传至CIBERSORTx作为混合文件，并利用S-mode来矫正不同来源的数据之间的批次效应；

(5)在绝对模式下使用所述单细胞RNA测序数据进行去卷积分析，利用基因表达特征矩阵来计算每个样本中的细胞亚群的绝对分数，保留在超过1/5的样本中绝对分数大于1的细胞亚群。

5.根据权利要求3所述的一种预测突变位点顺式调控基因的方法，其特征在于，在步骤S13中，

利用PEER软件从TMM归一化之后的单细胞RNA测序数据提取能够解释基因表达变化的隐藏因子，在合并时将隐藏因子纳入。

6.根据权利要求3或5所述的一种预测突变位点顺式调控基因的方法，其特征在于，在步骤S13中，

对样本的SNP基因型进行主成分分析，并提取基因型主成分，在合并时将基因型主成分纳入。

7.根据权利要求3所述的一种预测突变位点顺式调控基因的方法，其特征在于，步骤S14中，

利用Bonferroni方法对每个基因最显著的SNP的名义p值进行多重假设检验矫正，生成eigenMT矫正的p值；在eigenMT矫正的p值的基础上，利用Benjamini-Hochberg方法对所有基因的最显著的SNP进行多重假设检验矫正。

8.根据权利要求3或7所述的一种预测突变位点顺式调控基因的方法，其特征在于，

对于标准表达数量性状位点，首先利用与细胞亚群互作表达数量性状位点相同的方法对p值进行多重假设检验矫正获得在全基因组范围内的错误发现率，选取错误发现率<0.4标准表达数量性状位点进行后续分析。

9.一种计算机设备，其特征在于，包括：

存储器，用于存储计算机程序；

处理器，用于执行所述计算机程序时实现如权利要求1-8任一所述的一种预测突变位点顺式调控基因的方法的步骤。

10.一种计算机可读存储介质，其特征在于，

所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如权利要求1-8任一所述的一种预测突变位点顺式调控基因的方法的步骤。