HiChIP 数据分析: 差异 Loop 检测

差异 Loop 检测

为了识别由于热休克导致的染色质 3D 构象中的变异，我们将使用 R 包 diffloop 进行差异分析，该包实现了两种策略来评估可变 DNA Loop的显著性：负二项回归（来自 edgeR R 包）或加权 t 检验（来自 limma R 包）。

为了让 diffloop 正确执行统计分析，每个条件至少需要两个重复。diffloop 的输入文件是由 hichipper 生成的 interactions.all.mango 文件；因此，我们创建一个新文件夹，并将所有相互作用文件复制到其中：

cd $WORKDIR
mkdir Diffloop
mkdir Diffloop/input
all_dir=$(ls Hichipper/)
for dir in $all_dir; do
 cp Hichipper/$dir/*.interactions.all.mango Diffloop/input/
done

这会先创建 Diffloop/ 文件夹，然后创建 input/ 子文件夹；接着我们创建一个变量（all_dir），其中保存由 hichipper 生成的每个样本的文件夹名称；

最后，我们遍历所有文件夹，将所需文件复制到我们的 Diffloop/input/ 位置。

现在我们启动 R 并加载所需的包：

library(diffloop)
library(ggplot2)
library(ggrepel)
setwd("Diffloop")

为了读取四个分析样本的相互作用文件，我们使用loopsmake.mango函数：

bed_dir <- paste(getwd(),"input", sep="/")

sample_names <- c("Rad21_Rep1","Rad21_Rep2","HS_Rad21_Rep1","HS_Rad21_Rep2")

full <- loopsMake.mango(bed_dir, samples=sample_names)

由于我们有两套不同的 ChIP-Seq 峰值，锚点（anchors）在处理组与未处理组中是分别定义的。为了得到一份唯一的锚点区域清单，loops- Make.mango 函数还会把邻近的锚点合并（默认 gap 值为 500 bp）。用 dim 函数可以概览生成的 “loop” 对象（full），它包含 5 个不同的槽（slots）：存储唯一环锚基因组位置信息的 full@anchors（返回所有坐标的汇总视图）、记录已读入的所有样本间相互作用总数的 full@interactions（打印所有相互作用）、保存所有样本名称及分组信息的 full@colData、记录每个样本中支持各相互作用的读对数的 full@counts，以及最后记录环宽度的 full@rowData。

dim 汇总显示，full 包含 4 个样本，共有 51,022 个唯一环锚和 404,794 个已鉴定相互作用。在仅指定输入目录而未作其他设置时，import 函数会把所有样本归到同一组。因此，我们用 updateLDGroups 函数为每个样本设定正确类型：分别用 “NT” 和 “HS” 表示未处理和热休克处理样本；

注意，我们还通过 levels 参数指定了标签顺序。虽然自连环（self-ligated loops）在 hichipper 的最后步骤已被移除，但 import 函数因合并邻近环锚可能重新生成它们。因此，我们将先移除这些自连环，再过滤掉所有非显著相互作用（按 hichipper 中 Mango 实现的定义）：

full <- updateLDGroups(full, factor(c("NT","NT","HS","HS"), levels=c("NT","HS")))

full <- subsetLoops(full, full@rowData$loopWidth >= 5000)

full_qc <- mangoCorrection(full, FDR = 0.01)

如上所示，mangoCorrection 函数在 rowData 槽中新增两列，内容为 Mango 计算的 p 值 和 FDR。

作为最后一步过滤，我们要在每个组里剔除那些“在一套重复中强烈出现（PETs ≥ 5），却在另一套重复中完全缺席（PETs ＝ 0）”的Loop：

counts <- full_qc@counts

Rad21_dis <- ((counts[,"Rad21_Rep1"]>=5 & counts[,"Rad21_Rep2"]==0)|(counts[,"Rad21_Rep2"]>=5 & counts[,"Rad21_Rep1"]==0))

HS_Rad21_dis <- ((counts[,"HS_Rad21_Rep1"]>=5 & counts[,"HS_Rad21_Rep2"]==0)|(counts[,"HS_Rad21_Rep2"]>=5 & counts[,"HS_Rad21_Rep1"]==0)

qc_filt <- subsetLoops(full_qc, !(Rad21_dis | HS_Rad21_dis))

得到的 qc_filt 对象现在共包含 118,050 条唯一相互作用。为了评估重复之间以及样本之间的变异性，我们可以用以下函数计算并绘制主成分：

pcaPlot(qc_filt) + geom_text_repel(aes(label=sample_names)) +

ggtitle("PC Plot with Size Factor Correction") + theme(legend.position="none")

通过观察上图可见，未处理的 Rad21 HiChIP 重复聚在同一区域，而 HS 处理样本则表现出更大的离散度。正如预期，处理样本与未处理样本被明显分开。借助 loopMetrics 函数，我们可以评估所施加过滤步骤的影响：

一旦保留显著Loop，便可进行差异分析；由于仅有两组（NT 和 HS），可直接使用基于 edgeR 的 quickAssoc 函数：

nt_hs_Rad21_res <- quickAssoc(qc_filt)

dim(nt_hs_Rad21_res)

quickAssoc 会生成一个新的 “loop” 对象（nt_hs_Rad21_res），其 rowData 被扩展 5 列，分别给出 log2FC（按 HS vs NT 顺序计算的 Fold Change）、log2CPM、线性回归值、以及每条Loop对应的 p-value 和 FDR。topLoops 函数可进一步按 FDR 对 nt_hs_Rad21_res 进行子集化。本例以 1% FDR 为阈，得到 6477 条显著差异Loop，其中 4696 条在热休克后增强，1781 条减弱：

nt_hs_Rad21_res_sig <- topLoops(nt_hs_Rad21_res, FDR=0.01)

dim(nt_hs_Rad21_res_sig)

最后，我们将这些显著差异Loop以制表符分隔格式写出；summary 函数先把 “loop” 对象转成 data.frame，再用 write.table 保存：

write.table(summary(nt_hs_Rad21_res_sig), "diffLoop_sig_Rad21.txt", row.names=F, sep="\t", quote=F)

未完待续，欢迎关注！