一种单细胞染色质可及性图谱的跨物种分析方法、设备和介质

本发明属于单细胞染色质可及性数据分析，具体地，涉及一种单细胞染色质可及性图谱的跨物种分析方法、设备和介质。

背景技术：

1、细胞是生物体基本的结果和功能单位。不同类型的细胞共享同一个基因组序列，由细胞类型特异的转录程序控制，产生特定的细胞表型。在多细胞生物体内，这些功能特化的细胞类型往往存在于同一空间。传统的群体水平的测序技术，仅能检测得到细胞间信号的均值，这往往掩盖了单个细胞的特异信息。单细胞测序(single-cell sequencing)能够在单个细胞水平上，对基因组、转录组等遗传信息进行测序分析，从而揭示单个细胞特有的基因表达水平，以及染色质开放程度、拷贝数变异等重要信息。

2、近十几年来，单细胞测序技术蓬勃发展，大批优质的单细胞技术平台应运而生，覆盖了单细胞基因组、转录组、表观组、蛋白组和空间组等多个组学层面，整体上呈现出从单个细胞到百万细胞水平，从单组学到多组学联合分析的趋势。

3、蓬勃发展的单细胞测序技术推动了大规模的多物种全身细胞图谱的绘制。研究者们已先后完成对人、食蟹猴、小鼠等物种的单细胞转录组测序以及单细胞染色质转座酶可及性测序。物种水平单细胞图谱的绘制，揭示了细胞类型的特异性调控网络，挖掘了细胞类型广泛异质性背后的调控方式。目前跨物种图谱分析主要集中在单细胞转录组水平，通过挖掘物种间高度相似的转录调控程序，分析进化过程中细胞类型的保守性和多样性，从而得以构建跨物种细胞类型的进化层次结构树。而单细胞染色质可及性图谱的跨物种系统性比较研究尚未大规模开展，目前研究物种主要限制在人类、食蟹猴和小鼠等几个哺乳动物内部，尚未进行涉及其他脊椎动物与无脊椎动物的，多物种单细胞染色质可及性图谱数据集的联合分析研究。

4、因此，亟需利用现有生物信息学工具，开发并完善一套跨物种单细胞染色质可及性数据的集成分析策略。

技术实现思路

1、为了解决上述技术问题，本发明采用的技术方案如下：

2、本发明第一方面提供一种单细胞染色质可及性图谱的跨物种分析方法，包括以下步骤：

3、s1，单细胞染色质可及性测序数据的处理：

4、获得多个物种的单细胞染色质可及性测序数据并进行分层降维，在得到细胞分群结果后，使用谱系特异基因对细胞簇进行投射，确定细胞簇的谱系来源，将相同谱系的细胞簇汇聚在一起，得到谱系水平的细胞亚群，

5、s2，识别多物种的潜在调控元件：

6、s21，将细胞类型特异基因开放程度一致的细胞亚群进行合并，并注释成细

7、胞类型；

8、s22，以细胞类型为单位，使用macs2工具调用峰值(callpeak)，并采用迭代合并的计算方式，得到每个物种非重复的调控元件，

9、s3，调控元件的跨物种分析：

10、s31，选取所述多个物种中的一个作为参考物种，其基因组作为参考基因组，设置最小碱基匹配率为0.1，应用liftover工具与对应物种基因组转换的链文件，将其余物种的调控元件迁移到参考基因组上，能迁移到参考基因组的调控元件具有序列保守性，

11、对于参考物种，选取其余物种的任一种利用同样的方法获得所述参考物种的调控元件的序列保守性；

12、s32，将除所述参考物种之外的物种的调控元件与所述参考物种的调控元件取交集，计算重叠比例，将重叠比例大于等于0.5的调控元件定义为具有功能保守性。

13、在本发明的一些实施方案中，步骤s1中，进行分层降维前，进一步包括对物种的基因组进行窗口划分的步骤，具体根据基因组大小gs确定划分窗口的大小：

14、gs≥2g，窗口大小为1000bp；

15、1.5g≤gs<2g，窗口大小为800bp；

16、1g≤gs<1.5g，窗口大小为500bp；

17、gs<1g，窗口大小为50bp。

18、染色质可及性又称染色质开放程度，反映了染色质的转录活性状态，是研究基因表达调控的重要方向，在表观遗传图谱绘制、细胞分化和发育及各类疾病的发生发展研究中具有重要的作用。染色质的可及性，是指细胞核内大分子能够与染色质中dna发生物理接触的程度，由核小体或其他染色质结合因子在染色质上的占据情况和拓扑结构决定。

19、谱系特异性基因是指在某种物种中，某些基因在某种谱系的细胞中表达量高，且在其他谱系中不表达或者少表达的基因。同样的，在某种细胞类型的细胞中表达量高，且在其他细胞类型中不表达或者少表达的基因为细胞类型特异基因。这些基因对赋予各种类型细胞特异的形态结构特征与生理功能起到重要作用。

20、控制基因表达的调控元件大致上可以分成三类：promoter、enhancer、insulatory，也就是常说的启动子、增强子、绝缘子。在分化和发育过程中，这些调控元件会在高维的染色质结构中形成不同的相互作用。

21、macs基于一种新的模型可以很好的识别转录因子结合位点，调用峰值(callpeak)是macs2最主要的一个功能，能够利用bam文件寻找峰值(peak)。

22、scatac-seq数据的稀疏性使得降维非常困难。在本发明中，将scatac-seq数据矩阵二值化，也就是只有0和1两种情况(1代表开放，0代表不开放)。即便如此，因为转座酶切割的位置还是很少，所以二值矩阵大部分区域还是0。需要注意的是，scatac-seq数据中的0既有可能表示"不开放"，也可能是"没取样到"，从生物学角度上出发，这是两种截然相反的推断。正因为数据大部分都是0，所以1比0更有信息量。scatac-seq数据的稀疏性便是来源于这种低信息量。如果直接在稀疏矩阵上使用标准的降维分析，例如主成分分析(principal component analysis，pca)，然后使用前两个主成分进行绘图，可能无法得到想要的结果，这是由于数据的稀疏性会导致细胞间在0的位置有更高的相似性。为了解决这个问题，本发明采用分层的降维策略。

23、在本发明的一些实施方案中，步骤s1中，所述分层降维的步骤如下：

24、将窗口的开放程度作为特征，使用迭代的隐语义分析进行降维，对细胞进行第一轮无监督聚类；

25、以细胞簇为单位使用macs2工具调用峰值，以p<0.05的信号峰为特征，同样使用迭代的隐语义分析进行降维，对细胞进行第二轮无监督聚类；

26、利用umap的二维嵌入进行分群可视化，获得细胞的分群结果。

27、在本发明的一些实施方案中，步骤s21中，在对细胞亚群进行合并之前，进一步包括对细胞亚群中的细胞数量进行控制的步骤：不足300个细胞的细胞亚群按原本数目，超过300细胞的细胞亚群下采样到300个细胞。

28、在本发明的一些实施方案中，步骤s22中，进一步包括对调控元件进行分类的步骤：

29、将距离转录起始位点2kb以上的定义为远端调控元件；

30、将距离转录起始位点200bp与2kb之间的定义为近端调控元件；

31、将距离转录起始位点200bp以内的定义为启动子调控元件。

32、在本发明的一些实施方案中，步骤s3中，获得调控元件的序列保守性和功能保证性后，将所述调控元件分为三类：

33、序列不保守调控元件：不能迁移到参考基因组的调控元件；

34、序列保守性调控元件：能迁移到参考基因组但不具有功能保守性的调控元件；

35、序列与功能均保守调控元件：能迁移到参考基因组并且具有功能保守性的调控元件，

36、利用ucsc的bigwigaverageoverbed工具与从ucsc下载的多种物种的多序列比对结果文件，对各物种三类调控元件进行phylop打分，获得三类调控元件在跨物种水平上的保守性情况。

37、在本发明的一些具体实施方案中，所述功能保守性的判断方法如下：

38、对于某个物种经跨物种转换的调控元件，使用bedtools工具将其与参考物种的调控元件取交集；然后，计算重叠区间的比例(重叠区间的长度/该物种调控元件的长度)；最后，将该物种重叠比例大于等于0.5的调控元件定义为具有功能保守性。

39、在本发明的一些实施方案中，进一步包括：

40、s4，基于同源基因的开放模式，鉴定跨物种相似谱系与细胞类型：

41、s41，使用archr的“getmatrixfromproject”函数，获得各物种的染色质开放性基因得分矩阵；

42、s42，基于ncbi-blastp的基因蛋白序列相似性结果，获得物种间直系同源和旁系同源基因；

43、s43，利用图嵌入映射识别跨物种互近邻，计算单个细胞的20个最近的跨物种邻居细胞；利用sam算法构建细胞图谱的流式结构，将单细胞染色质可及性数据投射到一个统一的主成分空间，整合跨物种流式结构到统一的流式结构图，并于图中利用umap投影所述直系同源和旁系同源基因，

44、s44，利用基因间相关性更新基因-基因间的权重，得到跨物种谱系或细胞类型的匹配得分，根据得分是否大于0.2，筛选同源谱系对或细胞类型对，

45、s45，鉴定支持跨物种细胞类型同源性的基因对：对于任意两个细胞类型，首先，根据映射得分将基因对进行降序排列；然后，取前1000基因对作为驱动两个细胞类型相似的可能基因对集合；最后，对可能基因对集合进行筛选，筛选条件为：基因对中的两个基因必须分别是两个细胞类型的高变基因；

46、s46，获得多物种谱系特异性转录因子：基于调控元件的开放程度，采用js散度的方法来评估细胞谱系特异性的调控元件；使用chipseeker将调控元件进行注释到对应的基因，并从中筛选转录因子，从而获得多物种具有谱系特异性的转录因子的集合。

47、在本发明的一些实施方案中，进一步包括：

48、s5，使用scenicplus构建跨物种基因调控网络：

49、s51，获得各物种的motif注释文件、基因组区域-模体排序文件、物种染色体大小文件和基因注释文件；

50、s52，结合多物种单细胞转录组测序数据与单细胞染色质可及性测序数据，将基因表达矩阵和调控元件矩阵作为输入，

51、使用pycistopic工具无监督地识别具有差异可及性的调控元件，

52、使用pycistarget工具预测转录因子的结合位点，富集得到转录因子对应的模体，

53、使用“run_scenicplus”函数，根据转录因子调控区域、转录因子表达和靶基因表达的一致性联合预测调节因子，根据转录因子表达和调控区域富集分数的相关性≥0.2，筛选谱系特异性转录因子的调控子。

54、使用cytoscape，任选一个物种作为参考物种，呈现谱系特异性转录因子调控网络。

55、本发明第二方面提供一种计算机设备，包括：存储器，用于存储计算机程序；处理器，用于执行所述计算机程序时实现本发明第一方面任一所述的一种单细胞染色质可及性图谱的跨物种分析方法的步骤。

56、本发明第三方面提供一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现本发明第一方面任一所述的一种单细胞染色质可及性图谱的跨物种分析方法的步骤。

57、本发明的有益效果

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

59、本发明扩展了跨物种单细胞染色质可及性比较涉及的物种范围，开发并完善一套跨物种单细胞染色质可及性数据的集成分析方法，包括不同物种内部调控元件序列与功能保守性比较、不同物种之间同源基因与保守细胞类型对的鉴定与跨物种基因调控网络的构建。

60、本发明能够实现多物种的染色质可及性数据的整合分析，全面解析转录调控与基因表达的关系，从进化角度大规模筛选细胞谱系相关的转录调控元件。