三维基因组：diffHic 差异分析

diffHic

diffHic是一款 R 包，几乎包含了 Hi-C 数据分析的全部流程。它从对原始数据（. fastq 格式）进行比对开始，接着进行过滤、分箱、标准化，然后运用 HiCCUPS 方法来识别 loop，利用 directionality index 来识别拓扑关联域（TADs），并且基于 edgeR的统计框架开展差异分析。diffHic 的比对过程是通过一个 Python 脚本实现的，该脚本需要借助 Bowtie2，并且采用了一种适合较长 Hi-C 读长的嵌合比对策略。该工具的用户指南非常详尽，其本身也极具灵活性，即使在高分辨率下也能良好运行，因为数据是以 hdf5 格式存储，并以 GInteractions 对象的形式加载。GInteractions 是 InteractionSet R 包中一种类似 GRanges 的对象，专门用于存储成对的基因组接触信息。不过，对于哺乳动物基因组，在 5kb 及更高分辨率下，为了防止出现内存（RAM）问题，差异分析需要按染色体逐一进行。

Example

首先，需将 .hic 矩阵转换为 GInteractions 对象。借助 HiCExplorer 命令行工具套件中的 hicConvertFormat 函数，先将 .hic 矩阵在指定分辨率下转换为 .cool 文件，再从中提取染色体 2 的矩阵，转换为 GInteractions 格式。

在 /home/Hi-C_Matrices/ 文件夹中，输入：

hicConvertFormat --matrices GSM3262956_D00_HiC_Rep1. hic GSM3262957_D00_HiC_Rep2. hic GSM3262964_D15_HiC_Rep1. hic GSM3262965_D15_HiC_Rep2. hic --outFileName GSM3262956_D00_Hi-C_Rep1. cool GSM3262957_D00_HiC_Rep2. cool GSM3262964_D15_Hi-C_Rep1. cool GSM3262965_D15_HiC_Rep2. cool --inputFormat hic --outputFormat cool --resolution 5000

除了输入（--matrices）和输出文件（--out-FileName）参数之外，还需要明确输入和输出的格式（--inputFormat 和 --outputFormat），以及从 .hic 文件中提取矩阵时所使用的分辨率（--resolution）。需要注意的是，一旦指定了分辨率参数，输出文件的名称中也会包含该分辨率值（例如，在这个例子中，对于第一个样本，文件名会变成 GSM3262956_D00_Hi-C_Rep1_5000.cool）。在进行第二轮转换时，我们需要设置与之前相同的参数，不过这次输入格式为 cool，输出格式为 ginteractions，并且只需提取我们打算分析的染色体的原始计数值（--load_raw_values）（--chromosome）。以下代码可用于递归生成所有染色体的 GInteractions 格式矩阵：

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
    hicConvertFormat --matrices GSM3262956_D00_HiC_Rep1_5000.cool GSM3262957_D00_HiC_Rep2_5000.coolGSM3262964_D15_Hi-C_Rep1_5000. cool GSM3262965_D15_HiC_Rep2_5000. cool --outFile- NameGSM3262956_D00_HiC_Rep1_"$chr"_5000.gi GSM3262957_D00_HiC_Rep2_"$chr"_5000.gi GSM3262964_D15_Hi-C_Rep1_"$chr"_5000.gi GSM3262965_D15_HiC_Rep2_"$ chr"_5000.gi --inputFormat cool --outputFormat ginteractions --load_raw_-values --chromosome $chr
done

现在可以在 R 中加载 GInteractions 格式的矩阵了。为了将其转换为与 diffHic 兼容的对象，还需要进行一些额外的步骤。此外，需要加载 diffHic 包及其部分依赖项。以每个样本的染色体 2 的矩阵为例，首先使用 read.delim 函数将其加载到 R 中，然后通过 convertToGI 函数将其转换为 GInteractions 对象：

# R # version 3.6.2
options(scipen=999)
library(diffHic)
library(edgeR)
library(csaw)
library(statmod)
bin.size <- 5000

hESC_rep1<-read.delim("/home/Hi-C_Matrices/GSM3262956_D00_HiC_Rep1_chr2_5000.gi.tsv", header=FALSE)
hESC_rep2<-read.delim("/home/Hi-C_Matrices/GSM3262957_D00_HiC_Rep2_chr2_5000.gi.tsv", header=FALSE)
cardio_rep1<-read.delim("/home/Hi-C_Matrices/GSM3262964_D15_HiC_Rep1_chr2_5000.gi.tsv", header=FALSE)
cardio_rep2<-read.delim("/home/Hi-C_Matrices/GSM3262965_D15_HiC_Rep2_chr2_5000.gi.tsv", header=FALSE)

# Function to convert data frame into GInteractions object,paste it into the R session to use it
convertToGI <- function(df){
    row.regions <- GRanges(df$V1, IRanges(df$V2,df$V3))
    # interaction start
    col.regions <- GRanges(df$V4, IRanges(df$V5,df$V6))
    # interaction end
    gi <- GInteractions(row.regions, col.regions)
    gi$counts <- df$V7 # Interaction frequencies
    return(gi)
}
hESC_rep1.gi <- convertToGI(hESC_rep1)
hESC_rep2.gi <- convertToGI(hESC_rep2)
cardio_rep1.gi <- convertToGI(cardio_rep1)
cardio_rep2.gi <- convertToGI(cardio_rep2)

当使用自定义函数 convertToGI 为每个样本生成染色体 2 的 GInteractions 对象之后，我们把它们合并到一个列表里，并找出所有样本中都存在的相互作用：

hic.gi <- list(hESC_rep1.gi,hESC_rep2.gi,cardio_rep1.gi,cardio_rep2.gi)
names(hic.gi)<-c("hESC_rep1.gi","hESC_rep2.gi","cardio_rep1.gi","cardio_rep2.gi")

# Get common regions between the two objects
combined <- unique(c(hic.gi$hESC_rep1.gi, hic.gi$hESC_rep2.gi,hic.gi$cardio_rep1.gi,hic.gi$cardio_rep2.gi))

# replace original regions with the common regions
replaceRegions(hic.gi$hESC_rep1.gi) <- regions(combined)
replaceRegions(hic.gi$hESC_rep2.gi) <- regions(combined)
replaceRegions(hic.gi$cardio_rep1.gi) <- regions(combined)
replaceRegions(hic.gi$cardio_rep2.gi) <- regions(combined)

# Get the matching indexes between the two objects
matched <- lapply(hic.gi, function(x) {
    match(x, combined)
})

接下来，需要创建一个矩阵，用于存储索引以及每个样本的原始接触频率计数。这个矩阵会被转换为一个 iSet 对象，该对象首先用于 diffHic 的标准化过程，然后用于差异分析环节：

# Create a count matrix (for interaction frequencies)
counts <- matrix(0, ncol = 4, nrow=length(combined)) # counts for unmatched bins set to zero

# fill in the counts for matched bins
counts[matched$hESC_rep1.gi,1] <- hic.gi$hESC_rep1.gi$counts
counts[matched$hESC_rep2.gi,2] <- hic.gi$hESC_rep2.gi$counts
counts[matched$cardio_rep1.gi,3] <- hic.gi$cardio_rep1.gi$counts
counts[matched$cardio_rep2.gi,4] <- hic.gi$cardio_rep2.gi$counts

# create the InteractionSet object
iset <- InteractionSet(counts, combined)
names(assays(iset))<-"counts"
interactions(iset) <- as(interactions(iset), "ReverseStrictGInteractions")
iset$totals <- colSums(assay(iset))

现在可以利用 edgeR 的 aveLogCPM 函数计算数据集的平均对数每百万计数（logCPM），并剔除低丰度的相互作用。

keep <- aveLogCPM(asDGEList(iset)) > 0
iset <- iset[keep,]

在去除低计数相互作用后，我们需要校正因样本制备过程差异而产生的趋势偏差。我们首先通过均值平均（MA）图来观察这些偏差，它们在图中表现为与丰度相关的趋势（Figure a）：

pdf("diffHic_chr2_D00vsD15_pre-norm_MAplot.pdf")
ab <- aveLogCPM(asDGEList(iset))
o <- order(ab)
adj.counts <- cpm(asDGEList(iset), log=TRUE)
mval <- adj.counts[,3]-adj.counts[,2]
smoothScatter(ab, mval, xlab="A", ylab="M", main="hESC (1) vs.cardiomyocyte (2)")
fit <- loessFit(x=ab, y=mval)
lines(ab[o], fit$fitted[o], col="red")
dev.off()

为了校正不同文库之间效率偏差的均值依赖性，我们使用 normOffsets 函数对数据集进行标准化，选择 loess 方法作为标准化方式，并且分别对靠近对角线的区域和地图的其他部分进行标准化处理：

neardiag <- filterDiag(iset, by.dist=2*bin.size)
nb.off <- matrix(0, nrow=nrow(iset), ncol=ncol(iset))
nb.off[neardiag] <- normOffsets(iset[neardiag,], method="loess", se.out=F)
nb.off[!neardiag] <- normOffsets(iset[!neardiag,], method="loess", se.out=F)
assay(iset, "offset") <- nb.off

接着，我们再次绘制 MA 图，以检查偏差是否得到了有效校正（Figure b）：

pdf("diffHic_chr2_D00vsD15_post-norm_MAplot.pdf")
adj.counts <- log2(assay(iset) + 0.5) - nb.off/log(2)
mval <- adj.counts[,3]-adj.counts[,2]
smoothScatter(ab, mval, xlab="A", ylab="M", main="hESC (1) vs.cardiomyocyte (2)")
fit <- loessFit(x=ab, y=mval)
lines(ab[o], fit$fitted[o], col="red")
dev.off()

随后，我们将 iSet 转换为可用于 edgeR 的 DGEList 对象（通过 asDGEList 函数），确定实验设计以明确要比较的条件（利用 model.matrix 函数），并通过估计 NB 分散度（estimateDisp）来模拟同一实验条件下重复样本之间的生物学变异性。这一步旨在模拟整个数据集的均值 - 分散度关系，同时通过 glmQLFit 函数（用于 QL 分散度）来捕捉每个相互作用的变异性及其分散度的估计不确定性：

iset_dgeList <- asDGEList(iset)

# Establish groups to compare
design<-model.matrix(~factor(c(1,1,2,2)))
colnames(design) <- c("Intercept","cardio")

# Estimate NB and QL dispersion
iset_dgeList <- estimateDisp(iset_dgeList, design)
fit.iset_dgeList <- glmQLFit(iset_dgeList, design, robust=-TRUE)

现在可以使用 glmQLFTest 进行差异相互作用分析，从该函数返回的基于 EdgeR 列表的对象中提取结果，根据调整后的 p 值、logFC 和 logCPM 对结果进行筛选，并最终将结果转换为 .bedpe 格式，以便在 JuiceBox 中进行可视化：

# Perform differential analysis
result.iset_dgeList <- glmQLFTest(fit.iset_dgeList, coef=2)

# Filter interactions for FDR, logFC and CPM
FDR_threshold <- 0.2
logFC_threshold <- 0.5
logCPM_threshold <- 0.5
results_df<-(topTags(result.iset_dgeList, n=nrow(result.iset_dgeList$AveLogCPM), adjust.method = "BH", sort.by = "logFC", p.value = FDR_threshold))@.Data[[1]]
results_df_logFC <- results_df[results_df$logFC< -(logFC_threshold) | results_df$logFC > logFC_threshold,]
results_df_logFC_logCPM <- results_df_logFC[results_df_logFC$logCPM < -(logCPM_threshold) | results_df_logFC$logCPM > logCPM_threshold,]
sigInts<-rownames(results_df_logFC_logCPM)

# Create a dataframe containing differential interactions coordinates, raw interaction counts for all samples and differential test statistics
dimnames(iset_dgeList@.Data[[1]])[2][[1]] <- c("hESC_rep1","hESC_rep2","cardio_rep1","cardio_rep2")

differentials_df <- cbind(as.data.frame(iset@interactions)[sigInts,], iset_dgeList@.Data[[1]][sigInts,], results_df_logFC_logCPM)
Nrow <- nrow(differentials_df)
differentials_df_bedpe <- data.frame(differentials_df[,1:3], differentials_df[,6:8], "name"=rep(".", Nrow), "score"=rep(".", Nrow), differentials_df$strand1, differentials_df$strand2, "color"=rep(".", Nrow), differentials_df[,12:20], stringsAsFactors=F)
differentials_df_bedpe[differentials_df_bedpe$logFC<0,11] <- "0,0,139" # RGB color code for dark blue
differentials_df_bedpe[differentials_df_bedpe$logFC>0,11] <- "255,140,0" # RGB color code for orange
write.table(differentials_df_bedpe, "diffHic_D00vsD15_chr2_5kb.bedpe", sep="\t", col.names=T, row.names=F, quote=F)

按照这些阈值，diffHic 在染色体 2 上发现了 850 个差异相互作用。diffHic 还包含一些用于可视化的函数，例如 plotPlaid 和 rotPlaid，它们可以展示包含显著相互作用的 Hi-C 矩阵部分，但这些函数只能提供数据的静态视图。

本文由mdnice多平台发布