一种ATAC-seq测序数据的生物信息分析方法

摘要

本发明提供一种ATAC-seq测序数据的生物信息分析方法，包括以下步骤：对ATAC-seq测序数据进行分析与质控；将分析与质控后的数据与参考基因组进行比对分析；对单样本Peak进行检测与统计；组内一致性peak提取与分析；各处理组peak合并以及多样本聚类分析；各处理组间共有和特有peak分析；组间peak丰度差异分析。本发明围绕常规ATAC-seq测序数据进行研究，构建了一个分析内容丰富的，能满足个性化需求的常规ATAC-seq测序数据分析流程。分析结果揭示了全基因组范围内的染色质开放区域的序列信息，并且可以帮助科研工作者进一步探索找到或预测参与基因组开放性高低变化的转录因子。

权利要求

1.一种ATAC-seq测序数据的生物信息分析方法，其特征在于：包括以下步骤：

S1.对ATAC-seq测序数据进行分析与质控；

S2.将分析与质控后的数据与参考基因组进行比对分析；

S3.对单样本Peak进行检测与统计；

S4.组内一致性peak提取与分析；

S5.各处理组peak合并以及多样本聚类分析；

S6.各处理组间共有和特有peak分析；

S7.组间peak丰度差异分析。

2.根据权利要求1所述的一种ATAC-seq测序数据的生物信息分析方法，其特征在于：所述步骤S1中对ATAC-seq测序数据进行分析与质控的方法为：对下机的原始数据进行过滤，去除含有adapter的数据、N比例大于10％的数据和质量值Q≤10的碱基数占整条read的

40％以上的数据。

3.根据权利要求1所述的一种ATAC-seq测序数据的生物信息分析方法，其特征在于：所述S2中比对分析具体包括：利用比对软件Bowtie2将步骤S1获得的数据比对到参考基因组进行比对结果统计，过滤掉比对到线粒体或叶绿体上的数据，确认比对质量合格后，将比对到基因组上唯一位置的序列，进行后续生物信息分析，这些生物信息分析包括：基因组测序深度积累分布、Reads相对TSS位置的分布、Reads在染色体上的分布以及插入片段分析。

4.根据权利要求1所述的一种ATAC-seq测序数据的生物信息分析方法，其特征在于：所述S3中对单样本Peak进行检测与统计，具体包括：在全基因组范围内进行peak扫描并进行单样本peak统计，包括：单样本peak长度分布，单样本peak深度分布，单样本peak富集程度分布，单样本peak显著程度分布，单样本peak在基因功能元件上的分布以及单样本peak在染色体上的分布分析。

5.根据权利要求1所述的一种ATAC-seq测序数据的生物信息分析方法，其特征在于：所述步骤S4中组内一致性peak提取与分析，具体包括：进行IDR分析及处理组一致性peak获得，然后进行组内peak相关基因分析和组内共有peak的TF motif分析，其中组内peak相关基因分析包括：组内共有peak在基因功能元件上的分布、组内共有peak相关基因分析、组内peak相关基因进行GO富集分析和KO富集分析、组内共有peak的TF motif分析；组内共有peak的TF motif分析包括TFs motif的从头预测和已知TF motif的富集分析。

6.根据权利要求1所述的一种ATAC-seq测序数据的生物信息分析方法，其特征在于：所述步骤S5中各处理组peak合并以及多样本聚类分析，具体包括：利用DiffBind软件进行组间peak的合并，得到各个处理组间peak的并集，并计算各个peak在各个样本中的丰度；然后进行主成分分析，即将样本包含的高维度信息降维成数个维度的综合指标，进行样本间的比较；最后进行聚类分析，即计算两个样品间的皮尔斯相关系数，将相关系数以热图形式进行展示样品间的相关性。

7.根据权利要求1所述的一种ATAC-seq测序数据的生物信息分析方法，其特征在于：所述步骤S6中各处理组间共有和特有peak分析，具体包括：首先，通过韦恩图分析得到不同比较组之间特有和共有的peaks；接着，对应比较组间共有或特有的peak，进行peak相关基因的分析以及peak相关基因的GO和KO富集分析；最后，对于某个处理组特有的peak，对其进行TF-motif分析，包括：TF-motif的denovo预测和已知TFs-motif的富集分析。

8.根据权利要求1所述的一种ATAC-seq测序数据的生物信息分析方法，其特征在于：所述步骤S7中组间peak丰度差异分析，具体包括：利用DiffBind软件进行差异peak统计并绘制出差异peak统计图及差异比较火山图，然后进行差异peak基因分析和差异peak相关TF-motif分析，其中差异peak基因分析包括：提取出差异peak相关基因，并对peak相关基因进行GO富集分析和KO富集分析；差异peak相关TF-motif分析是对TF-motif的denovo预测和已知TF-motif的富集分析。

说明书

一种ATAC-seq测序数据的生物信息分析方法

技术领域

[0001] 本发明涉及生物学技术领域，具体涉及一种ATAC-seq测序数据的生物信息分析方法，开展常规ATAC-seq测序数据分析流程的探索与开发。

背景技术

[0002] 染色质是遗传物质的载体。真核生物的核DNA并不是裸露的，而是与组蛋白结合形成染色体的基本结构单位核小体，核小体再经逐步地压缩折叠最终形成染色体高级结构(如人的DNA链完整展开长度约2m，经过这样的折叠就变成纳米级至微米级的染色质结构而可以存储在小小的细胞核)。而DNA的复制转录是需要将DNA的紧密结构打开，从而允许一些调控因子，如转录因子、其他调控因子等结合。这部分打开的染色质，被称为开放染色质。打开的染色质允许其他调控因子结合的特性称为染色质可及性(chromatin accessibility)。因此，染色质的可及性与转录调控密切相关。

[0003] 开放染色质的研究方法主要包括ATAC-seq以及传统的DNase-Seq、FAIRE-seq。ATAC-seq(Assay for transposase-accessible chromatin with high throughput sequencing)是2013年由斯坦福大学William J.Greenleaf和Howard Y.Chang实验室开发的用于研究染色质可及性的方法，其原理是利用Tn5转座酶容易结合在开放染色质的特性，然后对Tn5转座酶捕获到的DNA序列进行测序。目前，ATAC-seq已经成为研究开放染色质的首选方法。

[0004] ATAC-seq可以全基因组范围内检测染色质的开放程度，能够得到全基因组范围内的蛋白质可能结合的位点信息。广泛应用于转录因子结合分析、核小体定位、活性调控元件分布等，在表观遗传机制研究领域具有广阔的应用前景。

[0005] 目前，对于常规ATAC-seq测序得到的数据的分析流程还没有形成的标准。因此，利用ATAC-seq研究开放染色质亟需一种既能标准化操作，又能满足个性化需求的ATAC-seq测序数据分析方法。

发明内容

[0006] 本发明的目的在于为克服上述现有技术的不足之处而提供一种ATAC-seq测序数据的生物信息分析方法。

[0007] 为实现上述目的，本发明采取的技术方案为：

[0008] 一种ATAC-seq测序数据的生物信息分析方法，包括以下步骤：

[0009] S1.对ATAC-seq测序数据进行分析与质控；

[0010] S2.将分析与质控后的数据与参考基因组进行比对分析；

[0011] S3.对单样本富集区域(Peak)进行检测与统计；

[0012] S4.组内一致性peak提取与分析；

[0013] S5.各处理组peak合并以及多样本聚类分析；

[0014] S6.各处理组间共有和特有peak分析；

[0015] S7.组间peak丰度差异分析。

[0016] 进一步地，所述步骤S1中对ATAC-seq测序数据进行分析与质控的方法为：对下机的原始数据进行过滤，去除含有测序接头(adapter)的数据、N比例大于10％的数据和质量值Q≤10的碱基数占整条测序片段(read)的40％以上的数据。

[0017] 进一步地，所述S2中比对分析具体包括：利用比对软件Bowtie2将步骤S1获得的数据比对到参考基因组进行比对结果统计，过滤掉比对到线粒体或叶绿体上的数据，确认比对质量合格后，将比对到基因组上唯一位置的序列，进行后续生物信息分析，这些生物信息分析包括：基因组测序深度积累分布、数据相对TSS位置的分布、Reads在染色体上的分布以及插入片段分析。

[0018] 进一步地，所述S3中对单样本富集区域(Peak)进行检测与统计，具体包括：在全基因组范围内进行peak扫描并进行单样本peak统计，所述单样本peak统计包括：单样本peak长度分布，单样本peak深度分布，单样本peak富集程度分布，单样本peak显著程度分布，单样本peak在基因功能元件上的分布以及单样本peak在染色体上的分布分析。

[0019] 进一步地，所述步骤S4中组内一致性peak提取与分析，具体包括：进行IDR分析及处理组一致性peak获得，然后进行组内peak相关基因分析和组内共有peak的TF motif分析，其中组内peak相关基因分析包括：组内共有peak在基因功能元件上的分布、组内共有peak相关基因分析、组内peak相关基因进行GO富集分析和KO富集分析、组内共有peak的TF motif分析；组内共有peak的TF motif分析包括：TFs motif的从头(de novo)预测和已知TF motif的富集分析。

[0020] 进一步地，所述步骤S5中各处理组peak合并以及多样本聚类分析，具体包括：利用DiffBind软件进行组间peak的合并，得到各个处理组间peak的并集，并计算各个peak在各个样本中的丰度；然后进行主成分分析(PCA)，即将样本包含的高维度信息降维成数个维度的综合指标，进行样本间的比较；最后进行聚类分析，即计算两个样品间的皮尔斯(pearson)相关系数，将相关系数以热图形式进行展示样品间的相关性。

[0021] 进一步地，所述步骤S6中各处理组间共有和特有peak分析，具体包括：首先，通过韦恩图分析得到不同比较组之间特有和共有的peaks；接着，对应比较组间共有或特有的peak，进行peak相关基因的分析以及peak相关基因的GO和KO富集分析；最后，对于某个处理组特有的peak，对其进行转录因子(TF)motif分析，包括：TF-motif的denovo预测和已知TFs-motif的富集分析。

[0022] 进一步地，所述步骤S7中组间peak丰度差异分析，具体包括：利用DiffBind软件进行差异peak统计并绘制出差异peak统计图及差异比较火山图，然后进行差异peak基因分析和差异peak相关TF-motif分析，其中差异peak基因分析包括：提取出差异peak相关基因，并对peak相关基因进行GO富集分析和KO富集分析；差异peak相关TF-motif分析是对TF-motif的denovo预测和已知TF-motif的富集分析。

[0023] 本发明的有益效果：本发明围绕常规ATAC-seq测序数据进行研究，构建了一个分析内容丰富的，能满足个性化需求的常规ATAC-seq测序数据分析流程。分析结果揭示了全基因组范围内的染色质开放区域的序列信息，并且可以帮助科研工作者进一步探索找到或预测参与基因组开放性高低变化的转录因子。分析流程条理清晰、逻辑性强，分析结果通过网页版结题报告的形式进行展现，层次分明，并且设置有超链接进行帮助说明，有利于去更加深刻地理解分析内容和操作。另外，还支持查看报告展示结果中完整的统计数据信息等。

附图说明

[0024] 图1为本发明的ATAC-seq测序数据的生物信息分析方法流程图。

[0025] 图2为本发明实施例1中样本过滤频数分布图。

[0026] 图3为本发明实施例1中样本LA-1基因组测序深度累积分布统计图。

[0027] 图4为本发明实施例1中样本LA-1Reads相对于TSS的位置分布图。

[0028] 图5为本发明实施例1中样本LA-1TSS周围信号分布图。

[0029] 图6为本发明实施例1中插入片段的分布图。

[0030] 图7为本发明实施例1中样本LA-1peak长度分布图。

[0031] 图8为本发明实施例1中样本LA-1peak深度分布图。

[0032] 图9为本发明实施例1中样本LA-1peak富集倍数分布图。

[0033] 图10为本发明实施例1中样本LA-1peak显著程度分布图。

[0034] 图11为本发明实施例1中样本LA-1peak在基因功能元件上的分布饼图。

[0035] 图12为本发明实施例1中样本LA-1peak在染色体上分布图。

[0036] 图13为本发明实施例1中样本LA的共有peak在基因功能元件上的分布图。

[0037] 图14为本发明实施例1中样本LA的peak相关基因GO富集分类柱状图。

[0038] 图15为本发明实施例1中样本LA的peak相关基因KO富集分析气泡图。

[0039] 图16为本发明实施例1中MEME软件预测出的样本LA中显著的motif序列图。

[0040] 图17为本发明实施例1中各样本的motif富集气泡图。

[0041] 图18为本发明实施例1中各样本间相关性分析热图。

[0042] 图19为本发明实施例1中样本NC与LA组间peak韦恩图。

[0043] 图20为本发明实施例1中样本NC与LA间差异peak比较火山图。

[0044] 图21为本发明实施例1中样本NC与LA间差异peak的相关基因GO富集分类柱状图。

[0045] 图22为本发明实施例1中样本NC与LA间差异peak的相关基因KO富集柱形图。

具体实施方式

[0046] 以下通过特定的具体实例并结合附图说明本发明的实施方式，本领域技术人员可由本说明书所揭示的内容轻易地了解本发明的其它优点与功效。本发明亦可通过其它不同的具体实例加以施行或应用，本说明书中的各项细节亦可基于不同观点与应用，在不背离本发明的精神下进行各种修饰与变更。

[0047] 实施例1

[0048] 本实施例提供一种基于老年黄斑变性的ATAC-seq测序数据的生物信息分析方法，样本分为3组，NC(表示正常样本)、LA(表示晚期老年黄斑变性样本)、EA(表示早期老年黄斑变性样本)，每组样本都包含1、2两个生物学重复样本。具体流程如图1所示，该方法包括如下步骤：

[0049] 步骤S1：对老年黄斑变性的ATAC-seq测序数据进行分析与质控，为了保证数据质量，首先对下机的原始数据进行质控和过滤，过滤掉含有adapter的数据、含N比例大于10％的数据和低质量的数据(质量值Q≤10的碱基数占整条read的40％以上)，得到高质量的clean数据用于后续生物信息分析，reads过滤信息如下表1所示。并且，还会将样本的数据过滤信息进行统计，以统计表格的形式进行呈现，或者通过绘制样本过滤频数分布图或过滤分析比例图，图2为样本过滤分析频数分布图。另外，也会对过滤前后的碱基信息进行统计，并且绘制样本过滤前后的碱基组成分布图，更加直观地展现数据质量情况。

[0050] 表1 reads过滤信息统计表

[0051]

[0052] 注：ID表示样本名称；Clean Reads Num表示Clean Reads1+Clean Reads2的数目；HQ Clean Reads Num(％)表示HQ Clean Reads1+HQ Clean Reads2的数目；Read length表示测序片段长度；Adapter表示测序接头；Low quality表示低质量的reads数目。

[0053] 步骤S2:将分析与质控后的数据与参考基因组进行比对分析，具体步骤如下：

[0054] 步骤S2.1：比对基因组统计，利用比对软件Bowtie2将数据比对到参考基因组，过滤掉比对到线粒体或叶绿体上的数据，最终将比对到基因组上唯一位置的数据(唯一比对数据)用于后续的信息分析。

[0055] 步骤S2.2：基因组测序深度累积分布，以步骤S2.1得到的唯一比对序列为分析对象，分析其在参考基因组上的覆盖分布，统计基因组位点的深度信息，得到基因组上测序深度统计结果，以LA-1的基因组测序深度累积分布图为例，如图3所示。

[0056] 步骤S2.3：reads相对TSS位置分布，利用deep Tools软件对转录起始位点(TSS)上下游2kb区域的所有数据进行统计，设置大小为50bp的窗口，计算每个窗口内的平均数据深度，绘制出数据相对于TSS位置的分布折线图，同时也可以绘制TSS周围信号分布热图，以样本LA-1为例，其分布折线图、热图分别如图4、图5所示。

[0057] 步骤S2.4：reads在染色体上的分布，将每个样本中唯一比对、去重复后的数据比对上基因组上各个染色体(分正负链)的密度进行统计，具体步骤是用滑动窗口为10K，计算窗口内部比对到碱基位置上的数据的中位数，并转化为log2。

[0058] 步骤S2.5：插入片段分析，一般地，大多数DNA都不包含或只包含一个核小体短片段，同时也有一些包含多个核小体的长片段，在含量上会呈现明显的片段分布特征。利用ATACseqQC软件绘制各样品中插入片段分布图，如图6所示。其中，横轴表示插入片段长度，纵轴表示该长度的插入片段占所有插入片段的比例。

[0059] 步骤S3:单样本peak检测与统计，主要是利用MACS分析软件进行全基因组peak扫描。MACS分析的原理是首先将所有的数据都向3'方向延伸到插入片段的长度，在基因组上进行滑窗，计算该窗口的dynamicλ，λ的计算公式为：λlocal＝λBG(λBG是指背景区域上的数据数目)，得到一定长度的建模区域；根据此区域中所有唯一比对数据的情况，使用Poisson分布模型进行检验，计算该区域的显著性P值，并进行FDR校正，若高于阈值，则认为该区域是一个peak。阈值为q-value<0.05，并对peak在基因组上的位置信息，peak区域序列信息等进行分析，筛选出peak相关基因。主要包括对单个样本peak统计，单个样本peak长度、深度、富集倍数及程度分布，单样本peak在基因功能元件上的分布以及单样本peak在染色体上的分布等，其中以样本LA-1为例，其peak长度分布图、深度分布图、富集倍数分布图、显著程度分布图、在基因功能元件上的分布饼图、在染色体上的分布图分别为图7、8、9、10、11、12。

[0060] 步骤S4:组内一致性peak提取与分析，具体步骤如下：

[0061] 步骤S4.1，IDR分析以及处理组一致性peak的获得：ATAC-seq要求提供≥2个生物学重复的样本，需要针对各个生物学重复样本检测出的peak进行IDR分析，得到各处理组内重复样本间存在重叠的共有peak(common peaks)。进一步地，将IDR分析得到的IDR值＜0.05归类为高度一致的共有peak(Common Filter peaks)，一致性peaks统计结果如表2所示。

[0062] 表2一致性peaks结果统计表

[0063]

[0064]

[0065] 注：ID表示样本名称；Commonpeaks表示共有的peaks数目；CommonFilterPeaks表示通过IDR阈值得到的共有peaks数目；FilterRatio表示通过IDR阈值的共有的peaks数目占比；NC表示正常样本；LA表示晚期老年黄斑变性的人组织样本；EA表示早期老年黄斑变性的人组织样本，下文同。

[0066] 步骤S4.2，组内共有peak在基因功能元件上的分布：统计peak在基因各功能区域的分布情况，统计peak在基因各功能区域的分布情况，按照promoter(基因上游2k以内)、5'UTR、3'UTR、Exon、Intron、Downstream(基因下游300bp以内)、Intergenic(基因上游2k以外或基因下游300bp以外)的优先级顺序，确定peak与各个功能区域的对应关系。分析结果用分布统计表、分布饼图等进行展示，组内共有peak在基因功能元件上的分布如表3所示。同时，也可绘制共有peak在基因功能元件上的分布比例图或者共有peak相对TSS距离的分布比例图，其中样本LA的共有peak在基因功能元件上的分布饼图如图13所示。

[0067] 表3 peak在基因功能元件上的分布统计表

[0068]

[0069] 注：PeakCount表示peak区域reads count值；Promoter表示peak位于启动子区域的数目、比例；5’UTR表示peak位于5’UTR区域的数目、比例；3’UTR表示peak位于3’UTR区域的数目、比例；1st exon表示层peak位于1st exon区域的数目、比例；other exon表示peak位于other exon区域的数目、比例；1st intron表示peak位于1st intron区域的数目、比例；Other intron表示peak位于Other intron区域的数目、比例；downstream表示peak位于downstream区域的数目、比例；Dista intergenic表示peak位于Dista intergenic区域的数目、比例。

[0070] 步骤S4.3，组内共有peak相关基因分析：根据peak在基因组上的区域信息及基因注释信息，得到peak相关基因，可能会被调控发生转录。将剔除归类为Distal intergenic的peak(peak位于基因上游2k以外或基因下游300bp以外)对应的peak相关基因。

[0071] 步骤S4.4，GO富集分析:是将步骤S4.3中组内peak相关进行GO富集分析。分析结果可以用GO富集分类柱状图、GO富集气泡图和GO富集条形图等形式进行结果展示，其中样本LA的peak相关基因的GO富集分析以GO富集分类柱状图为例进行展示，如图14所示。

[0072] 步骤S4.5，KO富集分析:是将步骤S4.3中组内peak相关进行KO富集分析。分析结果可以用KO富集气泡图和KO富集条形图等形式进行结果展示，以样本LA的peak相关基因的KO富集气泡图为例如图15所示。

[0073] 步骤S4.6，组间共有peak相关的转录因子分析:主要包括TFs motif的从头(de novo)预测和已知TF motif的富集分析。ATAC-seq检测的开放区间中包括：已知TFs motif(数据库已报道)和新TFs motif(数据库未报道)。

[0074] 步骤S4.6.1，TF-motif的denovo预测：为了保证新TFs motif也可以被检索与发现，使用MEME Suite从头预测peak区域内的保守motif。用MEME软件从头预测较长motif，用Dreme软件从头预测较短motif。不同颜色代表不同碱基类型，字母的高度代表这个碱基的保守性(字母越高，其在整个位点的频率越高，越保守)。图16展示的是用MEME软件预测样本LA中显著的motif序列图。

[0075] 步骤S4.6.2，已知TFs-motif的富集分析：使用MEME Suite中AME(Analysis of Motif Enrichment)对各样本中已知的TF motif进行检测，并进行已知TF motif的富集程度分析。AME使用JASPAR CORE非冗余数据库作为输入序列富集的对照，用motif扫描输入的每个序列，利用Average odds score进行打分，计算每个序列的score值，对序列得分较高的进行统计学检验(默认Fisher's exact test)，得到具有显著p值(P<0.05)的motif，认为该motif是富集的TF motif。用样本TF motif数量统计表进行各样本TF motif数量展示，各样本分组TF motif富集统计表详细展示出了motif名称、别名、序列、校正后的P值等信息，其中样本NC motif的富集统计结果如表4所示。另外，还可从各样本中挑选基于P值最显著富集的前10个motif，然后取并集，用于绘制motif富集程度的气泡图，如图17所示。

[0076] 表4样本NC motif的富集统计结果

[0077]

[0078] 注：motif ID表示motif名称；motif alt ID表示motif别称；Consensus表示motif序列；adj-p-value表示校正后的p值；E-value表示表富集显著性的期望值；TP表示真阳性motif序列数目；％TP表示真阳性motif序列数占阳性motif序列数的比例；FP表示假阳性motif序列数目。

[0079] 步骤S5：各处理组peak合并以及多样本聚类分析，主要包括各种处理组peak合并以及peak丰度计算、主成分分析(PCA)和相关性热图分析，具体分析步骤如下：

[0080] 步骤S5.1，各种处理组peak合并以及peak丰度计算：利用DiffBind软件进行组间peak的合并，得到各个处理组间peak的并集，并计算各个peak在各个样本中的丰度,所有样本peak丰度统计表如5所示。

[0081] 表5所有样本peak丰度统计表

[0082]

[0083] 注：CHR表示peak所在的染色体名称；START表示peak的起始位置；END表示peak的结束位置。

[0084] 步骤S5.2，主成分分析(PCA)：利用主成分分析(PCA)，将样本包含的高维度的信息(成千上万个peak富集信息)，降维为数个维度的综合指标(主成分)，以便于进行样本间的比较，同时保证原始数据中包含信息尽可能多地被保留。

[0085] 步骤S5.3，相关性热图分析：首先计算出每两个样品间的皮尔斯(pearson)相关系数，再将这些相关系数以热图的形式来直观地展示任意两个样品之间的相关性，各样本间的相关性如图18所示。

[0086] 步骤S6，各处理组间共有特有peak分析：主要包括组间共有特有peak提取、共有特有peak基因分析及共有特有peak相关TF-motif分析。分析步骤和方法与上述步骤S4中的部分方法相似，其区别在于，本步骤首先是通过韦恩图分析得到不同比较组之间特有和共有的peaks，结果以韦恩图进行展示，直观明了，图19为NC组与LA组之间的peak分析结果。

[0087] 步骤S7，组间peak丰度差异分析：主要包括了丰度差异peak分析、差异peak基因分析以及差异peak相关TF-motif分析等内容。具体步骤包括：

[0088] 步骤S7.1，peak差异分析及比较火山图：利用DiffBind软件对各个比较组间的所有peak(两组peak的并集)差异peak分析，选择丰度差异倍数log2值>＝1，FDR<＝0.05筛选差异peak，并注释到差异peak相关基因。差异分析结果可以用统计表进行展示，表中详细列出比较组的名称，上调、下调差异peak数量，所有差异peak数量，所有差异peak信息表和显著差异peak信息表等信息。还可以以直方图的形式进行结果展示，根据各比较组的显著差异peak，可以进行火山图分析，其中样本NC与LA间差异peak比较火山图如图20所示。

[0089] 步骤S7.2,组间差异peak相关基因分析：其分析步骤与步骤S4.3相似，是对比较组间差异peak相关基因进行注释分析。

[0090] 步骤S7.3，GO富集分析：针对步骤S7.2中得到的比较组间差异peak相关基因进行GO富集分析，其中样本NC与LA差异peak相关基因的GO富集分类柱状图如图21所示。

[0091] 步骤S7.4，KO富集分析：针对步骤S7.2中得到的比较组间差异peak相关基因进行KO富集分析，其中样本NC与LA差异peak相关基因的KO富集分析条形图如图22所示。

[0092] 步骤S7.5，差异peak相关TF-motif分析:主要是对TF-motif的denovo预测和已知TF-motif的富集分析。TF-motif的denovo预测包括1)上调peak的motif denovo预测：用MEME软件，从头预测较长的motif(8-15bp)；利用Dreme软件，从头预测较短的motif(3-8bp)；2)下调peak的motif denovo预测：用MEME软件，从头预测较长的motif(8-15bp)；利用Dreme软件，从头预测较短的motif(3-8bp)。不同颜色代表不同碱基类型，字母的高度代表这个碱基的保守性(字母越高，其在整个位点的频率越高，越保守)。已知TF-motif的富集分析主要是分别对上调、下调peak的motif进行富集分析，分析结果主要包括上调/下调TF-motif数量统计表、比较组之间motif富集统计表和比较组富集气泡图等。分析结果展现形式多样，直观明了。

[0093] 以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。