一种联合基因组三维结构差异鉴定和转录组基因表达水平差异分析挖掘功能基因的方法

摘要

本发明涉及生物信息技术领域，具体涉及一种联合基因组三维结构差异鉴定和转录组基因表达水平差异分析挖掘功能基因的方法。本发明提供一种基因组三维结构差异的鉴定方法，包括在Hi‑C数据的基础上进行差异AB Compartment的鉴定、差异TAD的鉴定和差异Loop的鉴定，利用该方法可实现利用原位Hi‑C数据对动植物进行不同层级的基因组三维结构差异的高效鉴定。本发明提供利用该基因组三维结构差异鉴定方法联合转录组基因表达水平差异分析挖掘功能基因的方法，利用该联合分析方法可实现从基因组三维结构与转录组两个层面进行差异分析，进而实现对动植物功能基因的高效挖掘。

说明书

技术领域

本发明涉及生物信息技术领域，具体涉及一种基因组三维结构差异的鉴定方法以及联合基因组三维结构差异鉴定和转录组基因表达水平差异分析挖掘功能基因的方法。

背景技术

基因组包含一维的DNA序列，目前已经有大量物种完成基因组测序。但是，对于基因组DNA不断缠绕折叠形成的染色体高级结构，目前仅有少数物种的研究。

基因组的三维结构主要包括从隔室(A/B Compartment)到拓扑相关结构域(TAD)，最后再到环(loop)的三级层级结构。隔室主要采用第一主成分进行鉴定，通常规定正的特征值与松散染色质状态相关联，称为A Compartment；负的特征值与压缩染色质状态相关，称为B Compartment。拓扑相关结构域(TAD)是指一个高度自关联的连续区域，并且其通过明显的边界与其相邻区域分离开来，形成一个独立的调控单位。TAD的边界具有重要作用，将边界部分删除后会导致基因调控紊乱，使得原来沉默的基因被转录，而原来应该转录的基因则被沉默。TAD作为一个调控单位，其内部的基因拥有共同的调控元件，因而存在其内部各基因的协同表达特征(为染色体上临近基因的共表达提供了依据)。Loop是指如果一对DNA片段的交互频率高于线性上相邻的片段的交互频率，那么这对片段将会形成一个染色质环。染色质环的一端连接启动子，另一端连接增强子，称为增强子-启动子环，可以实现对基因表达的远端调控。综上所述，基因组的三维结构在影响基因功能方面具有重要作用。

Hi-C技术是3C(Chromosome Conformation Capture)技术的衍生技术，它实现了全基因组范围内的染色体片段间的相互作用的捕获。Hi-C通过将空间结构临近的DNA片段进行交联，并将交联的DNA片段富集，然后进行高通量测序，对测序数据进行分析，揭示染色体片段间的交互信息，阐述基因组三维结构。

目前，分析差异三维基因组结构的软件只有juicer一款，但是，该软件只能分析差异Loop，且由于该软件是针对哺乳动物开发的，参数设置极为严格，在其它物种、尤其植物的差异Loop鉴定时差异Loop数量会特别少。目前尚没有一种较为完善的通用的差异三维基因组结构鉴定与功能基因挖掘的方法，因此亟需开发具有较高通用性的、不同层级的三维基因组结构(差异AB Compartment、差异TAD和差异Loop)差异的分析方法，并通过联合三维基因组结构差异与转录组基因表达差异进行基因功能富集分析，实现功能基因的精准挖掘。

发明内容

为解决现有技术中的技术问题，本发明的目的在于提供一种基因组三维结构差异的高效鉴定方法以及利用该基因组三维结构差异鉴定方法联合转录组基因表达水平差异分析高效挖掘功能基因的方法。

为实现上述目的，本发明通过对大量动植物和微生物来源的样本进行基因组三维结构信息的分析比较，开发了一套能够适用于各物种、不同层级的基因组三维结构鉴定(差异AB Compartment、差异TAD和差异Loop)的方法。进一步基于该鉴定方法开发了联合基因组三维结构差异鉴定和转录组基因表达水平差异分析挖掘功能基因的方法。

具体地，本发明的技术方案如下：

本发明提供一种基因组三维结构差异的鉴定方法，包括在Hi-C数据的基础上进行差异AB Compartment的鉴定、差异TAD的鉴定和差异Loop的鉴定。

优选地，在进行差异AB Compartment的鉴定、差异TAD的鉴定和差异Loop的鉴定之前，利用bwa和HiC-Pro进行Hi-C数据的比对和质控，获得有效的Hi-C数据。

本发明中，所述差异AB Compartment的鉴定为通过获得待比较的2个样本的代表性AB Compartment鉴定结果，比较每个样本的AB Compartment在基因组上的位置，获得差异AB Compartment结果。

具体地，所述差异AB Compartment的鉴定包括如下步骤：

(1)获取待比较的2个样本的各生物学重复的AB Compartment的鉴定结果；

(2)保留每个样本中各生物学重复一致的鉴定结果，得到每个样本的代表性ABCompartment结果；

(3)通过评估基因和ATAC-seq数据是否显著富集在A Compartment，基因组转座子元件(TE)和5mC甲基化数据是否显著富集在B Compartment，确定A和B Compartment鉴定的准确性；

(4)如果每个样本的A和B Compartment的鉴定结果均准确，则对每个样本的ABCompartment在基因组上的位置进行比较，获得2个样本的AB Compartment差异结果；

优选地，上述步骤(1)中，使用HiTC获取待比较的2个样本的各生物学重复的ABCompartment的鉴定结果。

优选地，上述步骤(3)中，所述显著富集的判断标准为p<0.01。具体地，使用bedtools工具统计基因组基因、ATAC-seq数据在A Compartment的分布以及基因组TE密度、5mC甲基化数据在B Compartment的分布。使用Two-sided Wilcoxon rank sum test检测显著性，若基因、ATAC数据显著富集在ACompartment(p<0.01)，TE、5mC甲基化显著富集在BCompartment(p<0.01)，则判断该样本的A和B Compartment鉴定准确。

优选地，上述步骤(4)中，若2个待比较样本的A和B Compartment均鉴定准确，则对每个样本的AB Compartment在基因组上的位置进行比较。将2个样本的AB Compartment比较结果划分为AB，BA，AA，BB和--五种类型，其中，AA表示两个样本都是A的区域；BB表示两个样本都是B的区域；AB表示一个样本是A，在另一个样本中变为B的区域；BA表示一个样本是B，在另一个样本中变为A的区域；--表示任意一个样本由于Hi-C数据覆盖、本身基因组序列是GAP或者是生物学重复之间不一致导致不能够判断A或者B的情况。

本发明中，所述差异TAD的鉴定采用基于DI值的差异TAD的鉴定方法。

具体地，所述差异TAD的鉴定包括如下步骤：

(1)鉴定每个样本的TAD并计算每个TAD的IR值，保留IR>1的TAD结果；

(2)统计TAD边界富集的基因数量与TAD内部区域富集的基因数量，若TAD边界富集的基因数量明显高于TAD内部区域富集的基因数量，则判断TAD的鉴定准确的鉴定准确；

(3)计算每个TAD边界在每个样本中的标准化的DI值，同时对于不同样本中有重叠的TAD，按照优先保留高DI值的TAD的原则，得到一套非冗余的TAD结果；

(4)对2个样本TAD边界DI值的差异进行显著性评估，计算p值，并使用Benjamini-Hochberg方法进行p值调整；

当调整后的p值<0.1，且一个样本的DI≥20、另一个样本小于20，2个样本的DI值差距在2倍以上时，判断为差异的TAD边界；

当调整后的p值不满足p值<0.1，但一个样本的DI值超过另外一个样本的4倍时，也判断为差异的TAD边界。

优选地，上述步骤(1)中，使用TadLib鉴定每个样本的TAD，并利用HOMER软件计算每个TAD的inclusion ratio(IR)值，保留IR>1的TAD结果。

优选地，上述步骤(2)中，使用bedtools统计TAD边界富集的基因数量与TAD内部区域富集的基因数量；并通过随机模拟在基因组上选取相同个数的TAD边界，作为对照，判断模拟的TAD的基因富集情况。

上述步骤(3)、(4)中，所述DI值的计算可采用文献报道的软件(Zhang Y,Li T,Preissl S,et al.Transcriptionally active HERV-H retrotransposons demarcatetopologically associating domains in human pluripotent stem cells.Naturegenetics,2019:1-9.)，该软件可通过https://github.com/shawnzhangyx/cvdc_scripts获得。

优选地，上述步骤(4)中，计算步骤(3)得到的非冗余的TAD结果中的TAD边界在每个样本中的DI值，并采用LIMMA软件计算p值，使用Benjamini-Hochberg方法进行p值调整。

本发明中，所述差异Loop的鉴定采用基于交互强度的差异Loop的鉴定方法。

具体地，所述差异Loop的鉴定包括如下步骤：

(1)使用Fit-Hi-C对每个样本进行Loop鉴定并保留q-value≤0.00001的Loop作为候选的Loop位点；

(2)对于有重叠的Loop，按照优先保留q-value小的结果的原则，合并不同样本有重叠的Loop，得到一套非冗余的Loop鉴定结果；

(3)基于交互强度获得2个样本的差异Loop结果。

优选地，上述步骤(1)中，q-value越低，表明Loop鉴定的准确性越高。

优选地，上述步骤(3)中，使用HOMER软件基于交互强度分析，当交互强度差异在2倍以上且FDR≤0.1时，判断为差异Loop，以此获得2个样本的差异Loop。

本发明还提供一种联合基因组三维结构差异鉴定和转录组基因表达水平差异分析的功能基因挖掘方法，其通过将转录组差异表达基因信息分别与利用所述基因组三维结构差异的鉴定方法鉴定得到的差异AB Compartment、差异TAD和差异Loop进行关联分析，分别获得转录组基因表达差异中与差异AB Compartment中相对应的基因、与差异TAD相对应的基因以及与差异Loop相对应的基因。

具体地，所述联合基因组三维结构差异鉴定和转录组基因表达水平差异分析的功能基因挖掘方法包括如下步骤：

(1)将转录组差异表达基因分别与差异AB Compartment中的AB，BA，AA，BB和--五种类型进行关联，保留AB和BA类型对应的差异表达基因，获得差异AB Compartment对应的差异表达基因；

(2)将差异TAD边界向内外各延伸2个bin，与差异表达基因进行关联，获得差异TAD对应的差异表达基因；

(3)将差异Loop两端向左右各延伸1个bin，与差异表达基因的启动子取交集，获得差异Loop对应的差异表达基因；

(4)对步骤(1)、(2)和(3)获得的差异AB Compartment、差异TAD和差异Loop对应的差异表达基因进行GO term的富集分析，筛选与研究目标相关的GO注释及其对应的基因；

(5)对步骤(1)、(2)和(3)获得的差异AB Compartment、差异TAD和差异Loop对应的差异表达基因进行KEGG富集分析，筛选与研究目标相关的代谢通路及其对应的基因。

上述步骤(2)、(3)中，bin分别为鉴定TAD设置的bin大小和鉴定Loop设置的bin大小；具体bin的大小一般采用Hi-C数据的分辨率进行设置，分辨率可以采用HiC-Pro软件获得。

优选地，上述步骤(1)、(2)、(3)中，所述关联或所述取交集采用bedtools进行。

优选地，上述步骤(3)中，所述差异表达基因的启动子设置为基因转录起始位点上游2kb。

优选地，上述步骤(4)中，对GO term富集分析结果进行显著性检验，将p≤0.01的结果作为最终GO富集的结果。

优选地，上述步骤(5)中，对KEGG富集分析结果进行显著性检验，将p≤0.05的结果作为最终KEGG富集结果。

本发明的有益效果至少包括：

本发明提供一种适用于动、植物等不同物种的通用的基因组三维结构差异的鉴定方法，利用该方法可实现利用原位Hi-C数据对动、植物及微生物进行不同层级的基因组三维结构差异的高效鉴定，包括差异AB Compartment的鉴定与全基因组展示、差异TAD的鉴定与热图展示、差异Loop的鉴定与全基因组展示。

本发明提供利用该基因组三维结构差异鉴定方法联合转录组基因表达水平差异分析进行功能基因挖掘的方法，利用该联合分析方法可实现从基因组三维结构(ABCompartment、TAD和Loop)与转录组两个层面进行差异分析，通过对鉴定到的基因组三维结构和转录组基因表达水平均存在差异的基因进行GO和KEGG富集分析，可实现对动、植物功能基因的高效挖掘。

附图说明

图1为本发明实施例1中全基因组范围不同的Compartment的变化对应差异表达基因的差异表达倍数的箱式图，其中，AB表示H01样本(正常样本)是A，在H02样本(癌症样本)中变为B的区域；BA表示H01样本(正常样本)是B，在H02样本(癌症样本)中变为A的区域；Unchanged表示A和B在两个样本中均没有发生变化的。

图2为本发明实施例1中差异TAD边界的热图展示；其中，箭头标记的边界为差异的TAD边界。

图3为本发明实施例1中Loop的互作热图；其中，A为正常样本，B为癌症样本，圈圈代表Loop所在的位置。

图4为本发明实施例2中基因组三维结构差异对应转录组差异表达基因的GO富集分析；其中，纵坐标为GO term注释信息；横坐标为-log10(P-value)，P-value为Fisher进行显著性检验得到的p值；柱形图中的值为注释到此GO term上的基因数目。

图5为本发明实施例2中基因组三维结构差异对应转录组差异表达基因的KEGG富集分析；其中，横坐标为富集因子，表示关联到的基因中位于该pathway条目的比例与整个基因组全部注释基因位于该pathway条目的比例的比值；纵坐标为-log10(P-value)，P-value为Fisher进行显著性检验得到的p值。

具体实施方式

下面将结合实施例对本发明的优选实施方式进行详细说明。需要理解的是以下实施例的给出仅是为了起到说明的目的，并不是用于对本发明的范围进行限制。本领域的技术人员在不背离本发明的宗旨和精神的情况下，可以对本发明进行各种修改和替换。

下述实施例中所使用的实验方法如无特殊说明，均为常规方法。

下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。

实施例1基因组三维结构差异的鉴定方法

本实施例以脊索瘤和正常样本(每个样本各2个生物学重复，2个样本分别命名为H01和H02)的Hi-C数据作为分析对象，提供一种基因组三维结构差异的鉴定方法，具体如下：

1、有效Hi-C数据的获得

利用bwa和HiC-Pro对Hi-C数据进行比对和质控，仅保留有效的Hi-C数据，最终获得bin大小为100000、20000和10000bp的交互矩阵。

2、差异AB Compartment的鉴定

(1)使用HiTC、采用bin＝100000，获得每个样本的AB Compartment的鉴定结果；

(2)对每个样本的生物学重复单独进行A和BCompartment的鉴定，保留各生物学重复之间一致的鉴定结果；

(3)使用bedtools工具统计基因、ATAC-seq、基因组TE密度、5mC甲基化数据在A和BCompartment的分布，并使用Two-sided Wilcoxon rank sum test检测显著性；若基因、ATAC数据显著富集在A Compartment(p<0.01)且基因组TE密度、5mC甲基化数据显著富集在B Compartment(p<0.01)，则A和BCompartment的鉴定结果准确；

(4)若每个样本的A和BCompartment的鉴定结果均准确，使用Compare_AB_from_different_sample.pl对每个样本的AB Compartment在基因组上的位置进行比较；将2个样本AB Compartment的比较结果划分为AB，BA，AA，BB和--五种类型，其中，AA表示两个样本都是A的区域；BB表示两个样本都是B的区域；AB表示H01样本是A，在H02样本中变为B的区域；BA表示H01样本是B，在H02样本中变为A的区域；--表示任意一个样本由于Hi-C数据覆盖、本身基因组序列是GAP或者是生物学重复之间不一致导致不能够判断A或者B的情况。

3、差异TAD的鉴定

(1)使用TadLib、采用bin＝20000，鉴定每个样本的TAD；利用HOMER软件计算每个TAD的IR值，保留IR>1的TAD结果；

(2)使用bedtools统计TAD边界和随机模拟制作的TAD边界的基因数量，观察两者的差别程度，若TAD边界富集的基因数量明显高于TAD内部区域，则判断TAD鉴定结果正确；

(3)使用produce_DI_Score.pl脚本获得差异的TAD结果。

produce_DI_Score.pl脚本可通过商业途径购买获得，可以实现下述功能：计算每个TAD边界在每个样本中的标准化的DI值，同时对于不同样本中有重叠的TAD，按照优先保留高DI值的TAD的原则，获得非冗余的TAD结果；对2个样本TAD边界DI值的差异进行显著性评估，计算p值，并使用Benjamini-Hochberg方法进行p值调整；当调整后的p值<0.1，且一个样本的DI≥20、另一个样本小于20，2个样本的DI值差距在2倍以上时，判断为差异的TAD边界；当调整后的p值不满足p值<0.1，且一个样本的的DI值超过另外一个样本的4倍时，也判断为差异的TAD边界。

4、差异Loop的鉴定

(1)对每个样本进行Loop鉴定并保留q-value≤0.00001的Loop作为候选的Loop位点；

(2)使用Fit-Hi-C软件自带的merge-filter.sh脚本合并各样本有重叠的Loop，得到一套非冗余的Loop鉴定结果；

(3)使用HOMER软件基于交互的强度，获得2个样本的差异Loop结果，具体标准为：交互强度差异在2倍以上且FDR≤0.1时。

结果分析如下：

差异AB Compartment的鉴定结果：以“AB”和“BA”分别代表不同的差异AB的区域，同时利用转录组数据进行验证，结果如图1所示，结果表明A－>B分布有376个下调表达的基因，65个上调的，B－>A分布有38个下调表达的基因，303个上调表达的基因。这与A区域倾向于基因表达(上调)，B区域倾向于基因沉默(下调)是一致的，证明差异AB鉴定的准确性。

差异TAD鉴定结果：以其中的一个热图作为示例(图2)，结果表明，在220.8-222.5Mb左右的TAD，在热图中的上三角(代表正常样本)，呈现一个完整的，无明显分界的三角形，表明形成一个的TAD，下三角(代表癌症样本)，热图中有明显的分界，形成了两个三角形，表明形成2个TAD，证明，鉴定的差异TAD是准确的。

差异Loop的鉴定结果：以其中一个Loop的互作热图作为示例(图3)。由于整个热图展示了基因组之间的互作信号，且左上三角和右下三角信号是对称的，颜色越深信号越强，代表越有可能形成Loop。热图的右上三角显示的正常样本中标记圈圈的地方是鉴定到差异的Loop所在的位置，该位置信号相对周边更强(图片是对称的，可以对比左下三角)，而癌症样本中则没有差异的Loop，证明此处鉴定的差异Loop是准确的。

实施例2联合基因组三维结构差异鉴定和转录组基因表达水平差异分析的功能基因挖掘方法

本实施例以实施例1的基因组三维结构差异鉴定结果以及正常和脊索瘤样本的转录组学基因表达差异数据作为分析对象，提供一种联合基因组三维结构差异鉴定和转录组基因表达水平差异分析的功能基因挖掘方法，具体如下：

(1)使用bedtools将将差异表达基因与差异AB Compartment中的AB，BA类型进行关联取交集，获得差异AB Compartment对应的差异表达基因信息；

(2)将差异TAD的边界向内外各延伸40000bp，然后使用bedtools与差异表达基因进行关联，获得差异TAD对应的差异表达基因；

(3)将差异Loop的两端向左右各延伸10000bp的长度，然后使用bedtools与差异表达基因的启动子(基因转录起始位点上游2kb)取交集，获得差异Loop对应的差异表达基因的。

(4)对步骤(1)、(2)和(3)获取得到的差异AB Compartment、差异TAD和差异Loop对应的差异表达基因利用go_enrichment.pl进行GO term富集分析，将p≤0.01的结果作为最终GO富集的结果，得到相应的富集的GO term和对应的基因列表；筛选与研究目标相关的GO注释及对应的基因。

(5)对步骤(1)、(2)和(3)获取得到的差异AB Compartment、差异TAD和差异Loop对应的差异表达基因利用draw_KEGG_graph.pl进行KEGG富集分析，保留p≤0.05结果作为KEGG富集分析的最终结果；筛选与研究目标相关的代谢通路及对应的基因。

结果分析：

基因组三维结构差异对应转录组差异表达基因的GO富集分析中，挑选20个富集最显著的GO term(P-value最小的20条)进行展示，结果如图4所示，其中，-log10(P-value)越大，表明关联基因在该GO term中的富集越显著。结果表明，利用联合基因组三维结构差异鉴定和转录组基因表达水平差异分析方法鉴定到了一些与细胞发育和生成相关的GOterm，这与癌细胞的快速生长也是相关的，证明利用本发明提供的基因组三维结构差异鉴定和转录组基因表达水平联合分析方法鉴定相关功能基因的可靠性。

基因组三维结构差异对应转录组差异表达基因的KEGG富集分析中，挑选20个富集最显著的pathway条目(P-value最小的20条)进行展示，结果如图5所示，富集因子(enrichment factor)较大，表示Peak区域基因在该pathway中的富集程度越明显；-log10(P-value)越大，表明关联基因在该pathway中的富集越显著。结果表明，利用联合基因组三维结构差异鉴定和转录组基因表达水平差异分析方法鉴定的差异Loop富集在cancer、PI3K-Akt signaling等与癌症有关的信号通路上，证明利用本发明提供的基因组三维结构差异鉴定和转录组基因表达水平联合分析方法鉴定相关功能基因的可靠性。

虽然，上文中已经用一般性说明、具体实施方式及试验，对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。