HiChIP 数据分析: 鉴定 Loops

鉴定HiChIP Loops

由于 HiChIP 实验流程中包含免疫沉淀步骤，因此可以把 loops 定义为被靶向蛋白结合的区间（即 anchor regions）之间的相互作用。

因此，这些区域应以存在 ChIP peak 为特征。hichipper 要求以 ChIP peaks 信息作为定义 anchors regions 的起点；这些 peaks 既可以是用户提供的、来自独立 ChIP-Seq 实验的结果，也可以由 hichipper 利用 HiChIP 自身的一部分 reads 推断出来。

与其他基于 ChIP 的实验不同，HiChIP 协议还包括限制性内切酶的消化步骤，这会使得读段密度图谱因靠近酶切位点而出现强烈偏差，从而妨碍对代表 loop anchors 的 peaks（进而对 loops）的正确识别。

事实上，当某个 ChIP peak 位于一个较长 restriction fragment 的中间时，在 peak 附近几乎找不到支持该 peak 的 HiChIP Paired End Tags (PETs)，因为这些 PETs 主要聚集于 restriction site 附近，导致鉴定出的 loops 被低估。

为克服这一限制，hichipper 采用了一种片段感知 padding 方法：首先，通过固定窗口（默认左右各 500 bp）对 peaks 进行 padding，以考虑 peak 定义的不确定性；随后，通过将每个 peak 扩展至其重叠的 restriction fragment 或 fragments 的末端来定义 loop anchors。该方法的效果是增加了支持 anchors 之间相互作用所需的 PETs 数量。此外，hichipper 还根据互作 anchors 的基因组位置及距离对 loops 进行筛选：仅保留 intrachromosomal loops；anchors 之间的最小距离设为 5 kb，以避免纳入自连接 loops，而最大距离设为 2 Mb，因为更远距离的相互作用生物学意义较小。最后，hichipper 采用 Mango loop significance model以提高 HiChIP 数据中相互作用识别的准确性。

Mango 将基因组区域之间相互作用的可能性建模为它们基因组分离距离和读段深度的函数，随后使用 Benjamini-Hochberg 方法对观察到互作位点的概率进行多重检验校正。

hichipper 的输入是 HiC-Pro 的输出文件夹，其中包含已比对且有效的 read pairs 坐标文件。hichipper 既可通过提供样本描述文件（或 manifest 文件，.yaml）运行，也可通过在命令行中直接指定所有输入文件（call mode）运行。

样本描述文件必须包含三个变量：

1. ChIP-Seq peak 文件，或指示如何直接从 HiChIP 数据中推断 peaks
2. 用于 HiC-Pro 的 restriction fragment 坐标文件
3. HiC-Pro 输出文件夹的路径

本方案将使用样本描述文件方法。由于 manifest 文件中只能指定单个 ChIP-Seq peak 文件，因此我们将为每个 ChIP-Seq 实验准备一份 manifest 文件。

下面给出指向未处理 Rad21 ChIP-Seq peaks 的 .yaml 文件内容，该文件存储在名为 yaml 的文件夹中，并假设 hichipper 从 $WORKDIR 运行：

peaks:
- ChIP_peaks/Rad21_peaks.noBL.narrowPeak

resfrags:
- /Programs/HiC-Pro_2.11.1/annotation/DpnII_resfrag_hg19.bed

hicpro_output:
- HiC_Pro

这份 manifest 文件（Rad21.yaml）可用于两份未经处理的 Rad21 HiChIP 重复实验，使用方法如下：

cd $WORKDIR  && mkdir Hichipper

hichipper --keep-samples Rad21_Rep1 -l 50 --out RAD21_1 Rad21.
yaml

mv RAD21_1 Hichipper/

hichipper --keep-samples Rad21_Rep2 -l 50 --out RAD21_2 Rad21.
yaml

mv RAD21_2 Hichipper/

其中 --keep-samples 告诉 hichipper 仅处理指定的样本（Rad21_Rep1），而忽略在 HiC-Pro/hic_results/data 文件夹中找到的任何其他样本。

-l 指定 HiChIP 的 read 长度，而 --out 参数指示 hichipper 将所有生成的输出保存到哪里；如果该文件夹已存在，程序会报错并要求更改名称或删除现有文件夹。

在此，我们为每个样本创建一个文件夹，然后为方便起见将其移动到 Hichipper/ 文件夹中。

hichipper 会生成两个“每次运行”文件（RAD21_1.log 和一个 html QC 报告文件，RAD21_1.hichipper.qcreport.html），以及七个“每个样本”文件：

• Rad21_Rep1.stat 包含 PETs 的总体统计信息。
• Rad21_Rep1.anchors.bed 列出所有已识别 loop anchors 的基因组坐标。
• Rad21_Rep1.inter.loop_counts.bedpe 与 Rad21_Rep1.intra.loop_counts.bedpe 分别记录 inter-chromosomal loops 和 intrachromosomal loops。
• Rad21_Rep1.filt.intra.loop_counts.bedpe 包含所有 anchors 距离落在 min/max 截断值范围内的 intrachromosomal loops。
• Rad21_Rep1.interactions.all.mango 与 Rad21_Rep1.filt.intra.loop_counts.bedpe 中的 loops 完全一致，但额外附加了每条 loop 的 FDR 值，该值由 Mango实现的偏差校正算法计算。这些相互作用还以 R “loop” 对象的形式保存为 Rad21_Rep1-HiChIP.rds 文件。

在我们的示例中，没有修改任何参数。锚点之间最小/最大距离的默认值分别为 5 kb 和 2 Mb，但可以通过 --min-dist 和 --max-dist 参数分别进行修改。用户还可以修改片段感知填充方法，例如通过调节左/右峰延伸（--peak-pad）以增加可以映射到环的读段对的数量。

“Reads in Anchors” 是指唯一、已映射且两条读段都与环锚点重叠的读段对（注意，这里也计算了被 HiC-Pro 过滤掉的读段对，例如落在同一限制性片段上的读段对）；这应该能很好地估计 HiChIP 的质量。“Mapped unique quality valid” 是指存储在 HiC-Pro 保存的 .allValidPairs 文件中的去重有效 PETs。“Valid intra 5kb-2Mb” 表示在指定距离范围内参与染色体内相互作用的 PETs，是对原位连接效率的评估。最后，“Reads in Loops” 是指包含在环中的读段对，即唯一、染色体内、宽度在 5 kb 到 2 Mb 之间且锚点已映射的 PETs。如果“reads in loops”的比例较低，但“valid intra 5kb-2Mb”的比例较高，这表明 ChIP 效率相对较低。在案例中，所有百分比都与其它高质量 HiChIP 相当。

未完待续，欢迎关注！