三维基因组：multiHiCcompare 差异分析

引言

multiHiCcompare是 R 包 HiCcompare 的升级版，它能够一次性对多个 Hi-C 数据集进行联合标准化和比较。其标准化方法是基于原始 HiCcompare 策略的改进版本，并且利用了“ cyclic loess procedur”。虽然该程序计算量较大，但可以通过并行化处理来提高效率，既可以按染色体（cyclic_loess）进行并行化，也可以按染色体 - 块在基于距离的基础上（fastlo）进行并行化。简单来说，标准化是通过将数据表示为 MD（均值 - 差异）图来实现的，然后通过拟合一条 LOESS 回归曲线来消除偏差。MD 图是 MA 图的一种变体，其中 M 表示两个数据集之间相互作用频率的对数比，D 表示两个相互作用区域之间的距离。差异分析是通过将负二项分布模型拟合到 Hi-C 数据中来完成的，使用的工具是 edgeR 包，其方法与 diffHic 类似。与 diffHic 不同的是，multiHiCcompare 会先将数据按距离依赖的块分割（类似于标准化步骤），然后分别对每个块进行建模。当只有两种实验条件时，可以使用 hic_exactTest 进行差异分析；而当涉及多种实验条件时，则更推荐使用 hic_glm 函数。

输入

multiHiCcompare 接受以稀疏上三角格式输入的 Hi-C 矩阵（原始值），这些矩阵可以从 .hic、.cool 和 HiC-Pro 文件中轻松获取。使用 straw 从 .hic 文件中提取矩阵，并指定不进行标准化（NONE），从而获得原始的相互作用频率。

declare -a chrs=("chr1" "chr2" "chr3" "chr4" "chr5" "chr6" "chr7" "chr8" "chr9" "chr10" "chr11" "chr12" "chr13" "chr14" "chr15" "chr16" "chr17" "chr18" "chr19" "chr20" "chr21" "chr22" "chrX")
for chr in "${chrs[@]}"
do
    echo $chr
    straw NONE /home/HiC_Matrices/GSM3262956_D00_HiC_Rep1.hic "$chr" "$chr" BP 5000 > /home/HiC_Matrices/GSM3262956_D00_HiC_Rep1_"$chr"_NONE_5kb.txt
    straw NONE /home/HiC_Matrices/GSM3262957_D00_HiC_Rep2.hic "$chr" "$chr" BP 5000 > /home/HiC_Matrices/GSM3262957_D00_HiC_Rep2_"$chr"_NONE_5kb.txt
    straw NONE Matrices/GSM3262964_D15_HiC_Rep1.hic "$chr" "$chr" BP 5000 > /home/HiC_Matrices/GSM3262964_D15_HiC_Rep1_"$chr"_NONE_5kb.txt
    straw NONE Matrices/GSM3262965_D15_HiC_Rep2.hic "$chr" "$chr" BP 5000 > /home/HiC_Matrices/GSM3262965_D15_HiC_Rep2_"$chr"_NONE_5kb.txt
done

输出

输出文件包含三列：第一个相互作用 bin 的起始坐标、第二个相互作用 bin 的起始坐标以及它们的原始相互作用频率。

实战

接下来，可以进入 R 环境，加载所需的包以及染色体 2 的输入数据：

# R # version 3.6.2
options(scipen=999)
library(multiHiCcompare)

# Set up analysis parallelization
library(BiocParallel)
numCores <- 20
register(MulticoreParam(workers = numCores))

data("hg19_cyto")

hESC_rep1<-read.table('/home/Hi-C_Matrices/GSM3262956_D00_HiC_Rep1_chr2_NONE_5kb.txt', header=FALSE)
hESC_rep2<-read.table('/home/Hi-C_Matrices/GSM3262957_D00_HiC_Rep2_chr2_NONE_5kb.txt', header=FALSE)
cardio_rep1<-read.table('/home/Hi-C_Matrices/GSM3262964_D15_HiC_Rep1_chr2_NONE_5kb.txt', header=FALSE)
cardio_rep2<-read.table('/home/Hi-C_Matrices/GSM3262965_D15_HiC_Rep2_chr2_NONE_5kb.txt', header=FALSE)

需要注意的是，通过 BiocParallel 包中的 MulticoreParam 函数利用并行处理来加速分析，指定要使用的 numCores 核心数。此外，还加载了 hg19_cyto，其中包含了人类基因组 hg19 版本中需要排除的区域，例如着丝粒和端粒。

然后，需要在数据框中添加一个包含染色体信息的列，重新排列列的顺序，将染色体名称移到第一列，并使用 colnames 为数据框设置列名：

hESC_rep1$chr <- rep("2",nrow(hESC_rep1))
hESC_rep2$chr <- rep("2",nrow(hESC_rep2))
cardio_rep1$chr <- rep("2",nrow(cardio_rep1))
cardio_rep2$chr <- rep("2",nrow(cardio_rep2))

hESC_rep1 <- HESC_rep1[,c(4,1:3)]
hESC_rep2 <- HESC_rep2[,c(4,1:3)]
cardio_rep1 <- cardio_rep1[,c(4,1:3)]
cardio_rep2 <- cardio_rep2[,c(4,1:3)]

colnames(hESC_rep1) <- c(’chr’, ’region1’, ’region2’, ’IF’)
colnames(hESC_rep2) <- c(’chr’, ’region1’, ’region2’, ’IF’)
colnames(cardio_rep1) <- c(’chr’, ’region1’, ’region2’, ’IF’)
colnames(cardio_rep2) <- c(’chr’, ’region1’, ’region2’, ’IF’)

现在数据框已经准备好转换为一个 hicexp 对象。这是一个 S4 对象，multiHiCcompare 在整个分析过程中存储数据，并包含实验设计。同时，可以使用参数 zero.p 过滤掉几乎全为零的矩阵行（即接触频率几乎处处为零的行），通过指定参数 A.min 过滤掉代表性不足的相互作用（默认情况下，multiHiCcompare 会过滤掉原始相互作用频率小于 5 的相互作用），或者通过 remove.regions 过滤掉不感兴趣的相互作用（例如涉及着丝粒和端粒的相互作用）：

hicexp <- make_hicexp(hESC_rep1, hESC_rep2, cardio_rep1, cardio_rep2, groups = c(0, 0, 1, 1), zero.p = 0.8, A.min = 5, filter = TRUE, remove.regions = hg19_cyto)

标准化

接下来是标准化步骤，这里使用 fastlo 函数进行标准化，因为 cyclic_loess 更加耗时。

可以使用 MD_hicexp 可视化标准化后的 MD 图，这与 diffHic 展示的 MA 图非常相似：

hicexp <- fastlo(hicexp, verbose = TRUE, parallel = TRUE)

pdf("hicexp_chr2_D00vsD15_MDplot_postNormalization.pdf")
MD_hicexp(hicexp)
dev.off()

如果不想在标准化步骤中使用并行化，可以将 parallel 设置为 FALSE。

差异分析

数据标准化后，可以使用 exactTest 进行差异相互作用分析。可以使用 MD_composite 制作另一个 MD 图，以可视化两组之间的差异，并使用 pdf 函数将其保存到文件中。其中，D.range 参数可以根据距离对图进行子集划分，这对于较大的染色体非常有用：

hicexp <- hic_exactTest(hicexp, p.method = ’fdr’, parallel = TRUE)

pdf("hicexp_chr2_D00vsD15_MDplot_postDiffAnalysis.pdf")
MD_composite(hicexp, D.range =0.1)
dev.off()

结果保存

最后，可以使用 topDirs 函数提取并保存结果，指定对数变化倍数的阈值（logfc_cutoff 参数）、对数每百万计数（logcpm_cutoff 参数）、校正后的 p 值（p.adj_cutoff 参数），以及结果的格式（bed 或 pairedbed 格式；后者适用于在 JuiceBox 中进行可视化）。

FDR_threshold <- 0.2
logFC_threshold <- 0.5
logCPM_threshold <- 0.5
td <- topDirs(hicexp1, logfc_cutoff = logFC_threshold, logcpm_cutoff = logCPM_threshold, p.adj_cutoff = FDR_threshold, return_df = ’pairedbed’)
Nrow <- nrow(td)
td_bedpe <- data.frame(td[, 1:6], "name"=rep(".",  Nrow), "score"=rep(".",  Nrow), "strand1"=rep(".",  Nrow), "strand2"=rep(".",  Nrow), "color"=rep(".",  Nrow), td[, 7:11], stringsAsFactors=F)
td_bedpe[td_bedpe$logFC<0, 11] <- "0, 0, 139" # RGB color code for dark blue
td_bedpe[td_bedpe$logFC>0, 11] <- "255, 140, 0" # RGB color code for orange
write.table(td_bedpe, "multiHiCcompare_chr2_D00vsD15_5kb.bedpe", sep="\t", col.names=F, row.names=F, quote=F)

通过这些阈值，multiHiCcompare 在染色体 2 上发现了 17,237 个差异相互作用。

multiHiCcompare 还提供了许多有用的函数来可视化结果。例如，manhattan_hicexp 函数可以用来可视化 p 值（比如通过将方法参数设置为 addCLT），或者以曼哈顿图的形式展示差异相互作用的数量，从而识别沿染色体的差异相互作用区域的热点。

png("hicexp_chr2_ManhattanPlot_p-value.png")
Differential Analysis of Hi-C data 81
manhattan_hicexp(hicexp, method = "addCLT")
dev.off()

png("hicexp_chr2_ManhattanPlot_count.png")
manhattan_hicexp(hicexp, method = "count")
dev.off()