鸿源韬生物

CUT&Tag生信分析报告

1.项目简介

1.1 样本信息

1.2 实验原理及流程

Hi-C（High-throughput Chromosome Conformation Capture）是一种基于高通量测序的三维基因组研究技术，用于解析染色质在细胞核内的空间组织及远距离DNA相互作用(Durand NC et al., 2016)。该方法基于染色质构象捕获（Chromosome Conformation Capture, 3C）技术，通过全基因组范围内的染色质互作检测，揭示染色体拓扑结构，如染色质环（Chromatin loops）、拓扑关联结构域（Topologically Associating Domains, TADs）和A/B染色质区室（Compartments）等。

Hi-C 实验原理

染色体作为遗传信息的载体，在细胞核内并非随机分布，而是占据特定的核区域，并形成高度有序的三维结构(LANCTOT C et al., 2007)。染色质的空间构象在细胞的生长、分化及基因调控过程中起着关键作用，其变化可能直接影响基因的表达调控及生物体的发育(SEXTON T et al., 2007)。Hi-C技术源于染色体构象捕获（Chromosome Conformation Capture, 3C）技术，是一种全基因组范围内的三维基因组学研究方法。该技术以整个细胞核为研究对象，结合高通量测序和生物信息学分析，解析染色质DNA在核空间中的互作关系(VAN BERKUM N L et al., 2010)。染色质的空间相互作用对于维持细胞正常功能至关重要，不同染色质区域之间的互作可能影响基因的转录活性、增强子-启动子调控及染色质高级结构的形成(SELVARAJ S et al., 2013)。

Hi-C数据分析可以揭示基因组不同位点之间的交互关系。通常，基因组被划分为固定大小的bin（如5kb、10kb或更大的窗口），以计算两个bin之间的交互强度。交互频率较高的区域可能代表染色质环（chromatin loops）或拓扑关联结构域（TADs），这些结构对于基因表达调控至关重要(FORTIN J P et al., 2015)。Hi-C数据能够与多种组学数据（如ChIP-Seq、转录组测序等）进行系统性整合分析，从而在基因调控网络和表观遗传调控的框架下，深入解析基因表达调控的分子机制及其在生物性状形成中的关键作用。通过多层次的数据整合，研究者能够揭示远程基因调控元件的功能、疾病相关的染色质结构重排事件，以及细胞分化过程中染色质三维构象的动态演变规律。这种整合分析方法为理解复杂生物过程中的染色质空间组织与功能调控提供了重要的研究工具和理论依据。

1.3 实验流程

Hi-C 的实验流程主要包括以下步骤：
a.染色质交联（Crosslinking）：通过甲醛（formaldehyde）处理细胞，交联细胞内染色质，使空间上相互接近的DNA片段共价结合，从而保留其天然的三维结构。
b.基因组酶切（Digestion）：使用限制性内切酶（如HindIII、MboI或DpnII）切割基因组DNA，产生粘性末端，使得原本相互接近的DNA片段被切割成小片段。
c.粘性末端填充（End Filling）与生物素标记（Biotin Labeling）：采用DNA聚合酶填充酶切产生的粘性末端，同时加入生物素标记的dCTP（Biotin-14-dCTP），以便后续富集特定片段。
d.DNA片段连接（Ligation）：通过DNA连接酶（T4 DNA ligase）在稀释条件下进行片段连接，使得原本三维空间上相互靠近的DNA片段更可能发生连接，形成嵌合片段（chimeric fragments）。
e.去交联、DNA纯化（Reverse Crosslinking & Purification）：通过高温和蛋白酶K消除甲醛交联，释放DNA，同时进行RNA酶和蛋白酶处理，纯化DNA。
f.生物素富集（Biotin Pull-Down）：通过链霉亲和素（streptavidin）磁珠特异性富集含有生物素标记的片段，从而去除未成功标记的背景DNA，提高后续测序的特异性。
g.文库构建与测序（Library Preparation & Sequencing）：采用高通量测序平台（如Illumina NovaSeq）对富集后的DNA片段进行建库和测序，通常采用双端测序（paired-end sequencing）模式，以便更好地解析DNA片段的配对关系。

1.4 分析流程

获得测序原始数据（raw data）后，首先对原始数据进行过滤，获得高质量的测序数据（clean data），将测序数据（clean data）比对到项目物种的参考基因组上，对比对结果进行鉴定出可用的相互作用（valid pair）。根据valid pair在一定窗口（分辨率）大小下构建全基因组交互矩阵，然后统计特定分辨率矩阵的交互频率、统计 Trans 和 Cis 的数量，对染色体相互作用分析；之后在不同分辨率下进行A/B compartment、TAD 分析和Loop分析。

Hi-C生物信息学分析流程

2. 数据质控

我们交付的原始数据为fastq（简称fq）格式文件的压缩包，文件名后缀通常为 “.fq.gz”。交付数据前我们会计算每个压缩文件的md5值。在您拿到数据之后，请您先校>验每个压缩文件的md5值，Linux下可以在数据目录使用“md5sum -c <*md5.txt>”命令进行校验，Windows下可使用hashmyfiles等校验工具，如发现压缩文件md5值与附在数据文件目录下的md5文档中的不一致则说明文件可能在传输的过程中被损坏。数据文件大小为文件占用磁盘空间的大小，文件的大小通常与磁盘格式、压缩比例等因素有关，与测序数据量（碱基数）的多少无对应关系，因此对应PE测序的 read1和read2两个文件大小也可能不相同。

将高通量测序得到的原始图像数据经过Base Calling 转化为序列数据，即FASTQ格式，得到最原始的测序数据文件。FASTQ 格式文件可记录所测读段（read）的碱基及其质量分数。FASTQ 格式以测序读段为单位进行存储，每条读段占 4 行，第一行是序列标识（read ID）以及相关的描述信息，以“@” 开头；第二行即为碱基序列，长度由测序策略决定；第三行以“+”开头，后面是序列标示符、描述信息，或者什么也不加； 第四行是测序质量值（phred），与第二行一一对应，phred值以ASCII码标记，对应的 ASCII 值减去33，即为第二行对应碱基的测序质量值，示例如下：

@HWI-ST1276:71:C1162ACXX:1:1101:1208:2458 1:N:0:CGATGT
NAAGAACACGTTCGGTCACCTCAGCACACTTGTGAATGTCATGGGATCCAT
+
#55???BBBBB?BA@DEEFFCFFHHFFCFFHHHHHHHFAE0ECFFD/AEHH

测序错误率用e表示， 平台测得数据的碱基质量值用Qphred表示，则有：Qphred=-10log10(e)。软件中碱基识别正确率与Phred分值之间的简明对应关系见下表：

测序Reads的错误率往往会随着测序接近尾声而升高，这是由测序过程中化学试剂的消耗造成共有的特征。

2.1 原始数据质控

Hi-C实验基于第二代测序（NGS）平台完成，采用双端测序文库构建策略。我们需要对原始测序数据进行质量评估与过滤，以确保后续分析的可靠性。首先，我们使用FastQC（version 0.12.1）(Andrews, 2010)对原始测序数据（raw data）进行全局质量分析，包括碱基质量分布（Phred score）、碱基组成平衡性（base content uniformity）、重复序列比例（duplication level）及GC含量偏差等指标，以全面评估测序质量。
我们使用Fastp（version 0.24.0）(Chen et al., 2018)对原始测序数据进行以下过滤操作。
接头序列去除：识别并切除双端reads中的接头序列；
低复杂度序列过滤：剔除含模糊碱基（N碱基占比≥10%）的reads；
动态质量修剪：通过滑动窗口法（5 bp窗口步长）评估局部序列质量，当窗口平均Phred score小于20时，执行3'端截断；
长度筛选：保留长度≥25 bp的paired-end reads，长度不足的reads及其匹配reads（R1/R2）均被排除。
原始和过滤后质控结果请详见result/01.qc文件夹，raw为原始数据质控结果，clean为过滤后质控结果。

MD_R1MD_R2MN_R1MN_R2

❮

❯

图2.1 各个样本平均测序碱基质量分数，横坐标代表150 bp长度序列中各个位置，纵坐标为该位置平均的碱基质量值Q；盒形图中间的红线表示中位数（median value）；黄色部分代表四分位距（25-75%）；上下分割线代表 90%和 10%的上下临界值；蓝色的线代表碱基质量的平均值。

MD_R1MD_R2MN_R1MN_R2

❮

❯

图2.2 各个样本碱基平衡性，图中四条线代表A T C G在每个位置平均含量。理论上，A和T应该相等，G和C应该相等，且4种碱基平行且接近分布。正常情况下四种碱基的出现频率应该是接近的，而且没有位置差异。因此好的样本中四条线应该平行且接近。当部分位置碱基的比例出现 bias 时，即四条线波动较大时可能存在测序数据或者文库污染。如果所有位置的碱基比例一致的表现出bias 时，即四条线平行但分开，往往代表文库有 bias (建库过程或本身特点)，或者是测序中的系统误差。测序刚开始由于测序仪状态不稳定，在15bp之前很可能出现波动。

MD_R1MD_R2MN_R1MN_R2

❮

❯

图2.3 各个样本重复序列水平，测序深度越高，越容易产生一定程度的重复（duplication），这属于正常的现象。但如果duplication 的程度很高，就提示我们可能有 bias 的存在（如建库过程中由于 PCR 扩增引起的duplication）。横坐标为 reads 重复的次数，纵坐标为重复次数对应的 reads 占 unique reads 的比例，以unique reads 的总数作为 100%。这里，我们仅对文件前 2000000 个reads 进行统计：对长度小于75bp 的reads 将其截短为 50bp，用于统计重复。

2.2 过滤后数据质控

这里展示Fastp过滤后的数据质控结果，图片内容与上面raw data类似。

MD_R1MD_R2MN_R1MN_R2

❮

❯

图2.4 各个样本平均测序碱基质量分数，横坐标代表150 bp长度序列中各个位置，纵坐标为该位置平均的碱基质量值Q；盒形图中间的红线表示中位数（median value）；黄色部分代表四分位距（25-75%）；上下分割线代表 90%和 10%的上下临界值；蓝色的线代表碱基质量的平均值。

MD_R1MD_R2MN_R1MN_R2

❮

❯

图2.5 各个样本碱基平衡性，图中四条线代表A T C G在每个位置平均含量。理论上，A和T应该相等，G和C应该相等，且4种碱基平行且接近分布。正常情况下四种碱基的出现频率应该是接近的，而且没有位置差异。因此好的样本中四条线应该平行且接近。当部分位置碱基的比例出现 bias 时，即四条线在某些位置波动较大时，可能测序数据或者文库存在污染。当所有位置的碱基比例一致的表现出bias 时，即四条线平行但分开，往往代表文库有 bias (建库过程或本身特点)，或者是测序中的系统误差。一般测序的时候，刚开始测序仪状态不稳定，在15bp之前很可能出现波动。

MD_R1MD_R2MN_R1MN_R2

❮

❯

图2.6 各个样本重复序列水平，测序深度越高，越容易产生一定程度的重复（duplication），这属于正常的现象。但如果duplication 的程度很高，就提示我们可能有 bias 的存在（如建库过程中由于 PCR 扩增引起的duplication）。横坐标为 reads 重复的次数，纵坐标为重复次数对应的 reads 占 unique reads 的比例，以unique reads 的总数作为 100%。这里，我们仅对文件前 2000000 个reads 进行统计：对长度小于75bp 的reads 将其截短为 50bp，用于统计重复。

2.3 数据过滤结果统计

我们对数据过滤结果进行统计，如下表所示：

表 2.1数据过滤结果统计：
Sample：样品名称；
Raw_Total_Reads/Clean_Total_Reads：过滤前后样本总reads数量，单位为百万；
Raw_Total_Bases/Clean_Total_Bases：过滤前后样本总碱基数量，单位为百万；
Raw_Q20_Rate/Clean_Q20_Rate：过滤前后样本Q20碱基比例；
Raw_Q30_Rate/Clean_Q30_Rate：过滤前后样本Q30碱基比例；
Raw_GC_Content/Clean_GC_Content:过滤前后样本GC含量。

3. 比对参考基因组

我们将各样品过滤后的clean data的reads与参考基因组进行比对，获取Reads在参考基因组上的定位信息，这里使用的软件是集成在HiCPro(version 3.1)(Servant N et al., 2015)中的Bowtie2(version 2.4.5)(Langmead B. et al., 2018)。具体流程如下：分别将Reads1和Reads2与参考基因组进行单端比对；对于未能比对的reads，若检测到ligation-site（酶切后文库的连接位点），则截取ligation-site之后的部分进行二次比对。最后，合并两次比对结果，筛选出两端均唯一比对到基因组特定位置的Reads Pair。去掉未比对、自环、再连等非有效数据，以及去掉重复数据，之后根据valid reads pair在一定窗口（分辨率）大小下构建全基因组交互矩阵。

HiCPro比对流程

3.1 比对参考基因组情况

表 3.1 Bowtie2比对情况：
1.sample:样本名称;
2.Total_pairs_processed：输入的总pairs数量;
3.Unmapped_pairs：未比对上的pairs;
4.Low_qual_pairs：比对质量MAPQ低于10pairs;
5.Unique_paired_alignments：比对在基因组中唯一区域的pairs;
6.Multiple_pairs_alignments：比对在多个区域的pairs;
7.Pairs_with_singleton：双端数据中只有一端数据比对上的pairs;
8.Low_qual_singleton：低比对质量且单端比对;
9.Unique_singleton_alignments：单端比对中唯一比对pairs;
10.Multiple_singleton_alignments：单端比对中多重比对pairs;
11.Reported_pairs：可用的pairs数量;

3.2 分子类型统计

由于Hi-C建库的特殊性，Hi-C建库会产生不同的分子类型（BELTON J M et al., 2012），包括自环（Self Circle）、边缘悬挂（Dangling End）、有效互作对（Valid Pair）。不同的分子类型，read1和read2的方向存在差异，Valid Pairs呈现四种类型的比对方向，即FF(两个正向)，RR(两个反向),FR和RF，但read1和read2来自不同的酶切片段，Dangling End呈现FR，Self Circle呈现RF，此两种分子类型reads序列来自同一个酶切片段。

表 3.1 文库pairs类型：
1.sample:样本名称;
2.Valid_interaction_pairs：合格的pairs数量;
3.Valid_interaction_pairs_F/R：F代表前向链，R代表反向链，来自各个链的互作对数量;
4.Dangling_end_pairs：边缘悬挂，检测到的连接点的末端，但未找到与之连接的另一个末端的pairs;Religation_pairs：重连接对，由于DNA断裂后重新连接产生的pairs;
5.Self_Cycle_pairs：自环对，检测到的连接点与自身连接成环的pairs;
6.Single-end_pairs：单端对，只有一个末端与其他连接点连接的pairs;
7.Filtered_pairs：过滤对，数据分析过程中被过滤的pairs;
8.Dumped_pairs：丢弃对，实验或数据处理过程中被丢弃的pairs。

plotHiCFragment_MDplotHiCFragment_MN

❮

❯

图3.1 各样本文库pair类型统计。3个柱状图分别是所有比对的片段（All Pairs）、有效 3C 配对（Valid 3C Pairs）和无效 3C 配对（Invalid 3C Pairs）。其中，All Pairs 代表所有测序得到的 Hi-C 片段，部分成功比对到参考基因组，部分因低质量或无匹配区域未能比对。Valid 3C Pairs 是符合 Hi-C 实验要求的有效配对，通常经过质量控制和筛选，而 Invalid 3C Pairs 代表未能通过筛选的片段，可能由于比对错误、低质量数据或不符合 Hi-C 片段要求。不同颜色表示数据的不同分类，如唯一比对、多重比对、未比对和低质量比对等。

3.3 顺反作用统计

Cis表示Paired reads位于同一条染色体上的顺式相互作用（intra-chromosomal interactions），而Trans表示Paired reads位于不同染色体上的反式相互作用（inter-chromosomal interactions）。顺式相互作用可进一步分为长距互作（两端距离大于等于20 kb）和短距互作（两端距离小于20 kb）。在全基因组交互矩阵中，顺式相互作用表现为沿对角线的染色体内交互信号较强，而染色体间交互信号较弱。这种现象部分归因于染色体在细胞核中的物理分离，占据独立的核区域。基于此特征，可通过顺式与反式交互的比例初步评估数据质量。若矩阵中存在随机背景连接等噪音，可能导致顺式与反式交互比例接近。

表 3.2 顺反作用统计：
1.sample:样本名称;
2.valid_interaction:合格的pairs数量;
3.valid_interaction_rmdup:去重重复后可用pairs数量;
4.trans_interaction:染色质之间反式互作数量;
5.cis_interaction:染色质内部顺式互作数量;
6.cis_shortRange:短距离（小于20kb）顺式互作数量;
7.cis_longRange:长距离顺式互作数量。

plotHiCContactRanges_MDplotHiCContactRanges_MN

❮

❯

图3.2 各样本顺反作用统计。左侧Duplicates代表被去除的重复pairs，ValidInteractions代表可用的pairs。右侧代表顺式/反式作用比例。

3.4 样本相关性

我们计算各文库之间pearson相关性，并绘制成heatmap热图：

图3.3 各文库之间pearson相关性热图

3.5 比对可视化

全基因组范围染色体相互作用结果保存为hic格式文件(位于report/result/02.map文件夹中)，hic文件是压缩的⼆进制⽂件，无法直接作为文本打开查看。客户可以使用软件Juicebox对hic等文件进行可视化浏览。对于bedpe/bed/bam等格式可以使用IGV浏览，IGV浏览器使用方法可参考我们提供的使用说明文档IGV快速上手

4. 染色质互作分析

4.1 染色质互作频率

普遍而言，Hi-C实验获得的染色质相互作用频率（chromatin interaction frequencies, IFs）在一个给定范围内随着距离的增大而衰减，一般使用交互衰减指数（Interaction decay exponents, IDEs）来描述IFs衰减的斜率陡度。染色质相互作用强度随距离衰减图如下:

图4.1 相互作用强度随距离衰减曲线。如果出现末端翘起代表超过了单条染色体长度，末端部分可以忽视。

4.2 染色体间相互作用

一般情况下，绝大多数的位点间交互分布在同一染色体内部，占比较少的染色体间交互表明染色体并不是随机分布在细胞核内。因此我们进行了染色体间相互作用分析，以发掘不同染色体间交互的基本情况。我们采用Heatmap展示全基因组染色体之间的交互（Valid Reads Pair数取log值），横纵坐标表示不同的染色体；交互强度随颜色深度逐渐增强。各样本全基因组互作图谱请见report/result/03.matplot/wholegenome文件夹。

MD_150000_kr_wgmatMN_150000_kr_wgmat

❮

❯

图4.2 全基因组互作图谱

为了更直观的展示反式相互作用，我们用绘制弦图展示染色体间的交互，下图展示了交互强度排名前1000的bin pair之间的交互。Circos图外围表示基因组上不同染色体；Circos内的曲线表示交互强度前1000 bin pair的位置，各样本弦图请见report/result/03.matplot/ciscir文件夹。

cis_circos_MDcis_circos_MN

❮

❯

图4.3 染色体间反式相互作用弦图

4.3 染色体内部相互作用

我在这里采用Heatmap展示各个染色体内部之间的顺式相互作用。下面只展示部分染色体的热图，各样本各染色体热图请见report/result/03.matplot/chrmat文件夹。

MD_40000_kr_1MD_40000_kr_11MN_40000_kr_1MN_40000_kr_11

❮

❯

图4.4 部分染色体内部相互作用热图

5. A/B Compartments

Lieberman-Aiden等在人的Hi-C研究中发现，在Mb级别分辨率下基因组可以划分为两类compartments，称作A/B compartments，Loci之间的交互很多发生在相同的compartments(LIEBERMAN-AIDEN E et al., 2009)。研究发现A compartments和开放染色质区域联系比较紧密，B compartments和封闭染色质区域联系比较紧密。A/B compartments具有细胞特异性，在不同的组织细胞内能够发生转换，这种转换和基因表达调控有一定关系。我们使用HiCExplorer(version 3.7.2)(RAMIREZ F et al., 2018)进行A/B compartments鉴定。

5.1 样本A/B Compartments分析

对于Hi-C实验得到的染色体各位点间（Bin pair）的交互（Valid reads对数）矩阵，我们称作Raw 矩阵，也称作Observed矩阵。HiC实验文库构建的过程中，由于酶切片段间随机性连接的存在，Observed矩阵通常被认为含有背景噪音。要较为真实反映染色体各位点间（Bin pair）的空间交互强弱，Observed矩阵需进行Observed/Expected（O/E）标准化处理。之后对O/E标准化后的矩阵进行Pearson相关性分析，会发现热图存在明暗相间的两种区域，它们对应的基因组区域就是两种compartment。我们对Pearson Correlation矩阵进行分析，得到染色体各位置（Bin）的第一特征向量（First eigenvector）和第二特征向量（Second eigenvector）。一般Bin的First eigenvector的正负两类标签分别代表了A/B compartments的划分。如果First eigenvector不能进行划分，可借助Second eigenvector进行区分。Pearson Correlation 热图与第一特征向量展示图请见report/result/04.ab/pearson文件夹。

MD_40000_11_pearsonMD_40000_1_pearsonMN_40000_11_pearsonMN_40000_1_pearson

❮

❯

图5.1 部分染色体Pearson Correlation 热图与第一特征向量展示图。PC1中标示的正值代表一类compartments、负值代表另一类compartments。结合两类compartments的基因密度、基因表达或表观组学等数据，一般把基因密度、基因表达高或染色质开放性高的一类compartment定义为A compartments，另一类定义为B compartments。

5.2 样本A/B Compartments转换

我们根据前面鉴定到的AB Compartments结果，比较不同样本之间A/B的差异，如果染色体某区域在对照组中是“A”，而在实验组中是“B”，那么我们就称之为发生了“A到B的转变”，转变区域信息记录在bed文件中，请见report/result/04.ab/diffab文件夹。

ABtrans_pie_MD_MN

❮

❯

图5.2 样本之间AB转变区域统计饼图。AA:在实验组和对照组都是A的区域;BB:在实验组和对照组都是B的区域;AtoB:在对照组中为A，在实验组中为B的区域;BtoA:在对照组中为B，在实验组中为A的区域。

ABtrans_chr_MD_MN

❮

❯

图5.3 转变区域在各染色体分布。nodiff:在实验组和对照组都是A或都是B的区域;AtoB:在对照组中为A，在实验组中为B的区域;BtoA:在对照组中为B，在实验组中为A的区域。

6. TAD分析

6.1 各样本TAD分析

在40kb左右分辨率的Hi-C交互热图中，我们能够观察到热图中呈现明显的三角形结。这些三角形的结构被命名为拓扑相关结构域（Topologically associated domains，TAD）。这些结构域之间具有明显的边界，研究表明CTCF和Cohesin等结构蛋白对于TAD边界结构的形成具有重要作用。研究发现TAD结构在不同时空下（组织、发育阶段等）具有一定的保守性，同时也存在一些动态变化。TAD内部之间相互作用比内部与外部相比更加活跃，TAD的边界可能对内外部互作在空间有阻碍作用，当TAD边界变动时，该区域基因表达可能收到影响。我们使用HiCExplorer进行TAD鉴定。

图6.1 TAD结构示意图

在分析过程中我们选用了Insulation的算法来鉴定TAD的边界。Insulation的算法原理是:我们按照一定的分辨率将基因组划分成n个大小相同的bin，对于每个bin，我们以该bin为中心划定一个（500kb-3500kb）左右的范围作为“社区”，计算该bin与“社区内”其他bin交互的总数作为其“社区活跃度”（community contacts），为了衡量该bin在整条染色体的社区活跃水平, 会计算所有bin的“社区活跃度“的均值（mean），并通过公式log2（community contacts/mean）计算该 bin的insulation score（绝缘分数）。在TAD的边界，这些区域的insulation score较低，而TAD内部的区域互作程度较高，因此它们的insulation score 较高。

图6.2 绝缘分数计算示意图

我们使用HiCExplorer算法进行TAD鉴定，结果请见report/result/05.tad文件夹。
其中*domains.bed记录各个TAD位置，这是一组不重叠的基因组坐标文件;*_boundaries.bed/gff记录了TAD边界的位置，基因组坐标与鉴定的分辨率相对应，gff文件内容包含了bed内容并增加了p值和delta值；*_score.bedgraph是记录了TAD分离分数的可视化文件。chrplot文件夹下为各染色体TAD分布以及分离分数。

MD_11_tadMD_1_tadMN_11_tadMN_1_tad

❮

❯

图6.3 部分染色体TAD分布以及分离分数，全部染色体结果请见chrplot文件夹。下方折线图为Separation score（分离分数，即insulation score），上方染色体互作热图中TAD边界用黑线框标注。

6.2 样本间TAD差异

差异TAD分析原理是对于两个样本的TAD分别进行TAD左边界、TAD右边界和TAD内部Wilcoxon rank-sum test，3个区域p值全部小于0.01则认为两个TAD是不同的。完整结果请见report/result/05.tad/difftad文件夹。其中*accepted.diff_tad代表两个样本相同的TAD，*rejected.diff_tad代表样本特有的TAD，chrplot文件夹为各染色体差异TAD分布。

MDvsMN_11_tadMDvsMN_1_tad

❮

❯

图6.4 部分染色体样本特异性TAD分布。

7. Loop分析

7.1 各样本Loop分析

染色质环（Chromatin loop）是二维基因组坐标上并不相邻但交互信号强烈的染色质空间结构。我们使用Juicer(version 2.20)(Durand NC et al., 2016)中HiCCUPS (Hi-C Computational Unbiased Peak Search)功能对鉴定样本染色质环并比较样本之间差异loop。它的工作原理是检查Hi-C接触矩阵中的每个像素，并将像素中的接触数量与像素周围一系列区域中的接触数量进行比较，loop两端锚点区域接触频率应该显著高于四周区域。
各样本鉴定loop结果请见report/result/06.loop文件夹，*_loops.bedpe是各样本记录loop信息可视化文件，可在基因组浏览器IGV中可视化。

图7.1 HiCCUPS原理

表7.1 各样本bedpe文件。各列代表信息： 1.chr：第一个锚点所在的染色体编号。
2.x1：第一个锚点的起始基因组坐标。
3.x2：第一个锚点的结束基因组坐标。
4.chr2：第二个锚点所在的染色体编号（通常与chr相同）。
5.y1：第二个锚点的起始基因组坐标。
6.y2：第二个锚点的结束基因组坐标。
7.name：Loop的名称或标识符（缺失时用.表示）。
8.score：Loop的置信度评分（缺失时用.表示）。
9.strand1：第一个锚点所在DNA链的方向（缺失时用.表示）。
10.strand2：第二个锚点所在DNA链的方向（缺失时用.表示）。
11.color：可视化Loop时使用的RGB颜色代码。
12.observed：实际观测到的染色质相互作用次数。
13.expectedBL：基于基线模型（Baseline）计算的预期相互作用次数。
14.expectedDonut：基于环形模型（Donut）计算的预期相互作用次数。
15.expectedH：基于水平结构模型（Horizontal）计算的预期相互作用次数。
16.expectedV：基于垂直结构模型（Vertical）计算的预期相互作用次数。
17.fdrBL：基线模型下计算的错误发现率（FDR），评估Loop显著性。
18.fdrDonut：环形模型下计算的错误发现率（FDR）。
19.fdrH：水平结构模型下计算的错误发现率（FDR）。
20.fdrV：垂直结构模型下计算的错误发现率（FDR）。
21.numCollapsed：合并到当前Loop中的邻近区域数量。
22.centroid1：第一个锚点的质心坐标（中点位置）。
23.centroid2：第二个锚点的质心坐标（中点位置）。
24.radius：Loop覆盖的半径（以碱基对bp为单位）。

7.2 差异Loop分析

我们使用Juicer中的HiCCUPS功能对样本进行差异Loop分析。样本之间特异性loop鉴定结果请见report/result/06.loop/diffloop，A_B_A.bedpe代表样本AB相比较时，只属于A的loop信息文件，可在基因组浏览器IGV中可视化。

表7.2 样本特异性Loop bedpe文件。

8. 差异区域基因注释

我们对差异的TAD边界区域以及差异的Loop锚点区域进行基因注释，并进行GO和KEGG富集分析。如果没有差异分析则本节内容为空。

8.1 差异区域基因注释

我们首先将差异区域划分为500bp 大小的bin，之后使用R包ChIPseeker(version 1.36)(Wang et al., 2022)对各个bin进行注释，我们统计bin在各基因功能元件分布情况，并将各个bin与基因关联。本节结果请详见位于report/result/07.anno文件夹。difftad代表差异的TAD边界区域；diffloop代表差异的Loop锚点区域。"A_B_A_PeakAnno.csv"的含义是A与B比较，A的特异性区域的基因注释文件。