引言

本文展示了如何利用 hictoolsr、DESeq2 和 plotgardener 包來尋找和可視化兩種生物條件下差異loop結(jié)構(gòu)的流程。首先，需要對(duì)原始數(shù)據(jù)進(jìn)行預(yù)處理，生成每個(gè)生物重復(fù)樣本和條件對(duì)應(yīng)的.hic文件，并識(shí)別出顯著的loop相互作用。之后，將不同條件下的loop結(jié)構(gòu)進(jìn)行合并，并提取每個(gè)條件和重復(fù)樣本中l(wèi)oop錨點(diǎn)之間的相互作用頻率計(jì)數(shù)。接著，通過 DESeq2 調(diào)用差異loop結(jié)構(gòu)，并以聚合峰分析圖的形式呈現(xiàn)結(jié)果。

以 Ahn 等人 2021 年發(fā)表的論文 “Phase separation drives aberrant chromatin looping and cancer development” 為例，該研究使用了 GEO 數(shù)據(jù)庫中的 GSE143465 和 GSE143465 數(shù)據(jù)。論文探討了一種罕見融合蛋白在急性髓系白血?。ˋML）中的致癌機(jī)制。這種融合蛋白 NUP98-HOXA9（NHA9）包含一個(gè) DNA 結(jié)合域和一個(gè)內(nèi)在無序區(qū)域（IDR），IDR 能夠形成相分離的凝聚體，從而導(dǎo)致三維染色質(zhì)結(jié)構(gòu)的變化和基因表達(dá)失調(diào)。為了研究相分離如何影響染色質(zhì)結(jié)構(gòu)，研究人員將 NHA9 的 IDR 中的苯丙氨酸（F）氨基酸殘基突變?yōu)榻z氨酸（S），使其失去相分離能力，并在 HEK293T 細(xì)胞中分別表達(dá)這種突變體（FS 突變體）和野生型（WT）。隨后，對(duì) WT 和 FS 細(xì)胞系進(jìn)行了 Hi-C、ChIP-seq 和 RNA-seq 實(shí)驗(yàn)，以比較染色質(zhì)結(jié)構(gòu)、NHA9 結(jié)合情況以及基因表達(dá)對(duì)相分離的響應(yīng)。

本文將展示如何在 WT 和 FS 這兩種生物條件下尋找差異loop結(jié)構(gòu)，并利用聚合峰分析圖來呈現(xiàn)結(jié)果。

預(yù)處理數(shù)據(jù)

在尋找差異loop結(jié)構(gòu)之前，需要先將原始的 .fastq 文件處理成 .hic 文件，并識(shí)別出顯著的loop相互作用。這些步驟是在 hictoolsr 包之外完成的預(yù)處理工作。以下部分將介紹如何使用現(xiàn)有工具處理 Hi-C 數(shù)據(jù)并調(diào)用loop結(jié)構(gòu)。

dietJuicer 處理 Hi-C 數(shù)據(jù)

可以使用 dietJuicer 流程將原始 Hi-C fastq 文件轉(zhuǎn)換為 .hic 格式。具體操作分為兩個(gè)層面：

1. 每個(gè)重復(fù)樣本需要單獨(dú)處理，生成各自的 .hic 文件。這是為了提取用于 DESeq 分析的計(jì)數(shù)數(shù)據(jù)。需要注意的是，測序重復(fù)樣本在計(jì)數(shù)提取時(shí)也應(yīng)進(jìn)行合并。
2. 按照 dietJuicerCore 流程的說明來處理這些樣本。 接下來，將重復(fù)樣本合并成一個(gè) “mega” .hic 文件。這個(gè)文件將用于后續(xù)的loop結(jié)構(gòu)調(diào)用。按照 dietJuicerMerge 流程的說明來創(chuàng)建合并后的 .hic 文件。

使用 SIP call loop

在使用 dietJucierMerge 創(chuàng)建合并的 .hic 文件后，可以利用 SIP（Significant Interaction Peak caller）來識(shí)別loop相互作用。

使用方法如下：

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

示例（在 UNC 的 longleaf 集群上提交作業(yè)）：

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

另一個(gè)loop結(jié)構(gòu)調(diào)用工具 HiCCUPS 也曾用于調(diào)用 Ahn 等人 2021 年 GSE143465 數(shù)據(jù)中的loop結(jié)構(gòu)。不過，由于 hictoolsr 的合并功能需要使用 SIP 生成的特定列，因此 hictoolsr 的示例數(shù)據(jù)中包含了通過 SIP 調(diào)用的loop結(jié)構(gòu)。

合并loop結(jié)構(gòu)

在兩種條件下分別調(diào)用了loop結(jié)構(gòu)，目的是找出每個(gè)數(shù)據(jù)集中獨(dú)有的loop結(jié)構(gòu)。然而，通常會(huì)發(fā)現(xiàn)兩個(gè)數(shù)據(jù)集中存在一些重復(fù)的loop結(jié)構(gòu)。將這些重復(fù)的loop結(jié)構(gòu)合并起來非常重要，這樣可以避免對(duì)重復(fù)的loop結(jié)構(gòu)進(jìn)行重復(fù)測試，并且能夠確保捕捉到所有獨(dú)特的loop結(jié)構(gòu)。mergeBedpe() 函數(shù)通過使用 DBSCAN 算法，將那些在不同條件下略有偏移的重復(fù)loop結(jié)構(gòu)合并在一起，從而實(shí)現(xiàn)了這一目標(biāo)。

## 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 計(jì)數(shù)

DESeq2 在調(diào)用差異loop結(jié)構(gòu)時(shí)，需要從重復(fù)的 .hic 文件中提取一個(gè)計(jì)數(shù)表。以下代碼展示了如何通過每個(gè)重復(fù)的 .hic 文件的 GEO 鏈接遠(yuǎn)程獲取這些計(jì)數(shù)。當(dāng)然，也可以選擇下載這些文件，并直接提供每個(gè)文件的本地路徑（鑒于網(wǎng)絡(luò)的不穩(wěn)定性，這種方法更為推薦）。

## 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')

由于提取計(jì)數(shù)可能需要一些時(shí)間，因此已將示例數(shù)據(jù)集包裝在上面顯示的代碼的預(yù)提取計(jì)數(shù)中。

data("loopCounts")

提取計(jì)數(shù)取為提供為列的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

以下代碼演示了如何使用上面創(chuàng)建的Loopcounts對(duì)象使用DESEQ2調(diào)用差異loop。首先，將計(jì)數(shù)矩陣與loop計(jì)數(shù)分離：

## Load package
library(DESeq2)

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

head(cnts)

然后，使用Counts Matrix中的列名創(chuàng)建列數(shù)據(jù)：

## 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

接下來，可以構(gòu)建DESEQ數(shù)據(jù)集，并比較“ 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 探索差異結(jié)果：

plotMA(res)

或者 PCA：

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

接下來，將 DESeq2 的差異分析結(jié)果重新整合到 loopCounts 對(duì)象中，然后根據(jù)經(jīng)過 Benjamini-Hochberg（BH）校正后的 p 值 0.01 以及 log2 倍數(shù)變化高于或低于 0 的標(biāo)準(zhǔn)，分別篩選出野生型（WT）特異性和突變型（FS）特異性的loop結(jié)構(gòu)。

## 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多平臺(tái)發(fā)布