三维基因组：Loop结构 差异分析(1)

引言

本文展示了如何利用 hictoolsr、DESeq2 和 plotgardener 包来寻找和可视化两种生物条件下差异loop结构的流程。首先，需要对原始数据进行预处理，生成每个生物重复样本和条件对应的 .hic 文件，并识别出显著的loop相互作用。之后，将不同条件下的loop结构进行合并，并提取每个条件和重复样本中loop锚点之间的相互作用频率计数。接着，通过 DESeq2 调用差异loop结构，并以聚合峰分析图的形式呈现结果。

以 Ahn 等人 2021 年发表的论文 “Phase separation drives aberrant chromatin looping and cancer development” 为例，该研究使用了 GEO 数据库中的 GSE143465 和 GSE143465 数据。论文探讨了一种罕见融合蛋白在急性髓系白血病（AML）中的致癌机制。这种融合蛋白 NUP98-HOXA9（NHA9）包含一个 DNA 结合域和一个内在无序区域（IDR），IDR 能够形成相分离的凝聚体，从而导致三维染色质结构的变化和基因表达失调。为了研究相分离如何影响染色质结构，研究人员将 NHA9 的 IDR 中的苯丙氨酸（F）氨基酸残基突变为丝氨酸（S），使其失去相分离能力，并在 HEK293T 细胞中分别表达这种突变体（FS 突变体）和野生型（WT）。随后，对 WT 和 FS 细胞系进行了 Hi-C、ChIP-seq 和 RNA-seq 实验，以比较染色质结构、NHA9 结合情况以及基因表达对相分离的响应。

本文将展示如何在 WT 和 FS 这两种生物条件下寻找差异loop结构，并利用聚合峰分析图来呈现结果。

预处理数据

在寻找差异loop结构之前，需要先将原始的 .fastq 文件处理成 .hic 文件，并识别出显著的loop相互作用。这些步骤是在 hictoolsr 包之外完成的预处理工作。以下部分将介绍如何使用现有工具处理 Hi-C 数据并调用loop结构。

dietJuicer 处理 Hi-C 数据

可以使用 dietJuicer 流程将原始 Hi-C fastq 文件转换为 .hic 格式。具体操作分为两个层面：

1. 每个重复样本需要单独处理，生成各自的 .hic 文件。这是为了提取用于 DESeq 分析的计数数据。需要注意的是，测序重复样本在计数提取时也应进行合并。
2. 按照 dietJuicerCore 流程的说明来处理这些样本。 接下来，将重复样本合并成一个 “mega” .hic 文件。这个文件将用于后续的loop结构调用。按照 dietJuicerMerge 流程的说明来创建合并后的 .hic 文件。

使用 SIP call loop

在使用 dietJucierMerge 创建合并的 .hic 文件后，可以利用 SIP（Significant Interaction Peak caller）来识别loop相互作用。

使用方法如下：

java -jar SIP_HiC.jar hic <hicFile> <chrSizeFile> <Output> <juicerToolsPath> [options]

示例（在 UNC 的 longleaf 集群上提交作业）：

java -jar SIP_HiC_v1.6.1.jar hic file.hic hg19_chromSizes_filt.txt outdir juicer_tools.jar -g 2.0 -t 2000 -fdr 0.05

另一个loop结构调用工具 HiCCUPS 也曾用于调用 Ahn 等人 2021 年 GSE143465 数据中的loop结构。不过，由于 hictoolsr 的合并功能需要使用 SIP 生成的特定列，因此 hictoolsr 的示例数据中包含了通过 SIP 调用的loop结构。

合并loop结构

在两种条件下分别调用了loop结构，目的是找出每个数据集中独有的loop结构。然而，通常会发现两个数据集中存在一些重复的loop结构。将这些重复的loop结构合并起来非常重要，这样可以避免对重复的loop结构进行重复测试，并且能够确保捕捉到所有独特的loop结构。mergeBedpe() 函数通过使用 DBSCAN 算法，将那些在不同条件下略有偏移的重复loop结构合并在一起，从而实现了这一目标。

## Load packages
library(hictoolsr)
library(dbscan)

## Define WT and FS loop file paths
wt_loops <- system.file("extdata/WT_5kbLoops.txt", package = "hictoolsr")
fs_loops <- system.file("extdata/FS_5kbLoops.txt", package = "hictoolsr")

## Merge loops and convert to GInteractions
loops <- 
  mergeBedpe(bedpeFiles = c(wt_loops, fs_loops), res = 10e3) |>
  as_ginteractions()

head(loops)

提取 Hi-C 计数

DESeq2 在调用差异loop结构时，需要从重复的 .hic 文件中提取一个计数表。以下代码展示了如何通过每个重复的 .hic 文件的 GEO 链接远程获取这些计数。当然，也可以选择下载这些文件，并直接提供每个文件的本地路径（鉴于网络的不稳定性，这种方法更为推荐）。

## Hi-C file paths from GEO
hicFiles <- 
  c("https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM4259nnn/GSM4259896/suppl/GSM4259896_HEK_HiC_NUP_IDR_WT_A9_1_1_inter_30.hic",
    "https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM4259nnn/GSM4259897/suppl/GSM4259897_HEK_HiC_NUP_IDR_WT_A9_1_2_inter_30.hic",
    "https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM4259nnn/GSM4259898/suppl/GSM4259898_HEK_HiC_NUP_IDR_WT_A9_2_1_inter_30.hic",
    "https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM4259nnn/GSM4259899/suppl/GSM4259899_HEK_HiC_NUP_IDR_WT_A9_2_2_inter_30.hic",
    "https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM4259nnn/GSM4259900/suppl/GSM4259900_HEK_HiC_NUP_IDR_FS_A9_1_1_inter_30.hic",
    "https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM4259nnn/GSM4259901/suppl/GSM4259901_HEK_HiC_NUP_IDR_FS_A9_1_2_inter_30.hic",
    "https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM4259nnn/GSM4259902/suppl/GSM4259902_HEK_HiC_NUP_IDR_FS_A9_2_1_inter_30.hic",
    "https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM4259nnn/GSM4259903/suppl/GSM4259903_HEK_HiC_NUP_IDR_FS_A9_2_2_inter_30.hic")

## Extract Hi-C counts between loop pixels
loopCounts <- extractCounts(bedpe = loops,
              hic = hicFiles,
              chroms = c(1:22, "X"),
              res = 10e3,
              norm = 'NONE',
              matrix = 'observed')

由于提取计数可能需要一些时间，因此已将示例数据集包装在上面显示的代码的预提取计数中。

data("loopCounts")

提取计数取为提供为列的Filepath的Basename。让简化列名：

## Load package
library(InteractionSet)

## Simplify column names
colnames(mcols(loopCounts)) <- 
  gsub(pattern = "GSM.*_IDR_(WT|FS)_A9_(1|2)_(1|2)_.*", 
       replacement = "\\1_\\2_\\3",
       x = colnames(mcols(loopCounts)))

head(loopCounts)

差异 loop calling

以下代码演示了如何使用上面创建的Loopcounts对象使用DESEQ2调用差异loop。首先，将计数矩阵与loop计数分离：

## Load package
library(DESeq2)

## Isolate count matrix
cnts <- 
  mcols(loopCounts)[grep("WT|FS", colnames(mcols(loopCounts)))] |>
  as.matrix()

head(cnts)

然后，使用Counts Matrix中的列名创建列数据：

## Create colData from column names
colData <- 
  do.call(rbind, strsplit(x = colnames(cnts), split = "_")) |>
  as.data.frame(stringsAsFactors = TRUE) |>
  `colnames<-`(value = c("condition", "biorep", "techrep"))

colData

接下来，可以构建DESEQ数据集，并比较“ WT”和“ FS”条件之间的差异loop：

## Build DESeq data set
dds <- 
  DESeqDataSetFromMatrix(countData = cnts,
                         colData = colData,
                         design = ~techrep + biorep + condition)

## Run DEseq analysis
res <-
  DESeq(dds) |>
  lfcShrink(coef = "condition_WT_vs_FS", type="apeglm")

summary(res)  

可以使用 MA-plot 探索差异结果：

plotMA(res)

或者 PCA：

plotPCA(vst(dds), intgroup = "condition") + ggplot2::theme(aspect.ratio = 1)

接下来，将 DESeq2 的差异分析结果重新整合到 loopCounts 对象中，然后根据经过 Benjamini-Hochberg（BH）校正后的 p 值 0.01 以及 log2 倍数变化高于或低于 0 的标准，分别筛选出野生型（WT）特异性和突变型（FS）特异性的loop结构。

## Attach DESeq2 results
mcols(loopCounts) <- cbind(mcols(loopCounts), res)

## Separate WT/FS-specific loops
wtLoops <- loopCounts[loopCounts$padj <= 0.01 &
                         loopCounts$log2FoldChange > 0]

fsLoops <- loopCounts[loopCounts$padj <= 0.01 &
                         loopCounts$log2FoldChange < 0]

summary(wtLoops)
summary(fsLoops)

本文由mdnice多平台发布