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

通过聚合峰分析进行可视化

既然已经找出了“WT”和“FS”条件之间的差异loop结构，就可以利用聚合峰分析（APA）来直观地展示loop结构调用的质量。APA 是一种以 Hi-C 数据中的中心loop像素为中心，展示周围接触频率矩阵的堆叠图。

计算 APA

APA 的计算过程是提取并聚合围绕 Hi-C 像素的矩阵，该像素具有特定的分辨率（res）以及在每个方向上的像素数量（buffer）。举个例子，若要提取一个 21x21 的矩阵，且分辨率为 10-kb，就需要将 res 设置为 10e3，buffer 设置为 10。那些过于靠近对角线的“短”相互作用需要被过滤掉，以防止聚合时出现错误。filterBedpe() 函数会计算出哪些相互作用会与对角线相交，并将其过滤掉。在下面的代码中，将loop结构的不同类别整合成一个列表，并针对 10-kb 的分辨率和 10 的缓冲区进行这种过滤操作：

## Assemble all, WT, and FS loops into a list
loopList <- 
  list(allLoops = loopCounts,
       wtLoops = wtLoops,
       fsLoops = fsLoops)

## Define resolution and buffer (pixels from center)
res <- 10e3
buffer <- 10

## Filter out short loop interactions
filteredLoops <- 
  lapply(X = loopList,
         FUN = filterBedpe,
         res = res,
         buffer = buffer) |>
  `names<-`(value = names(loopList))

lapply(filteredLoops, summary)

你会发现很多相互作用都靠近对角线，因此被过滤掉了。接下来的代码部分将展示如何将 calcApa() 函数应用于经过筛选的loop结构列表，从而从“WT”和“FS”条件的合并重复样本 Hi-C 文件中提取和聚合计数：

## Hi-C file paths from GEO
wtHicPath <- "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE143nnn/GSE143465/suppl/GSE143465_HEK_HiC_NUP_IDR_WT_A9_megaMap_inter_30.hic"
fsHicPath <- "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE143nnn/GSE143465/suppl/GSE143465_HEK_HiC_NUP_IDR_FS_A9_megaMap_inter_30.hic"

## Calculate APA matrices for loops from WT Hi-C data
loopApaWtHic <-
  lapply(X = filteredLoops,
         FUN = calcApa,
         hic = wtHicPath,
         norm = "KR",
         buffer = buffer)

## Calculate APA matrices for loops from FS Hi-C data
loopApaFsHic <-
  lapply(X = filteredLoops,
         FUN = calcApa,
         hic = fsHicPath,
         norm = "KR",
         buffer = buffer)

由于计算 APA 矩阵可能需要较长时间，因此已经提供了一个示例数据集，其中包含了上述代码中预先计算好的 APA 矩阵。

data("loopApaWtHic")
data("loopApaFsHic")

lapply(loopApaWtHic, dim)
lapply(loopApaFsHic, dim)

在对这些矩阵进行可视化之前，最后一步是将总和值根据每个类别中的loop结构数量进行归一化处理，这样可以将解释转变为每个loop结构的平均信号。这也有助于在可视化时使各图的尺度保持一致。

## Get the number of loops for each condition
nLoops <- lapply(filteredLoops, length)

## Divide each matrix by nLoops
loopApaWtHic <- Map("/", loopApaWtHic, nLoops)
loopApaFsHic <- Map("/", loopApaFsHic, nLoops)

使用 ggplot2 进行可视化

若想通过 ggplot2 来呈现结果，得先将矩阵转换成长格式。

## Convert matrix to long-format
long <- 
  loopApaWtHic$allLoops |>
  as.table() |>
  as.data.frame() |>
  setNames(c('rows', 'cols', 'counts'))

## Visualize with ggplot2
library(ggplot2)
ggplot(data = long,
       mapping = aes(x = rows, y = cols, fill = counts)) +
  scale_fill_distiller(palette = 'YlGnBu', direction = 1) + 
  geom_tile() +
  theme(aspect.ratio=1, axis.text.x = element_text(angle = 45, hjust=1))

另外，通过在行或列上使用 rev() 函数，可以翻转矩阵，进而改变 Hi-C 对角线的方向。

## Flip the matrix
library(ggplot2)
ggplot(data = long,
       mapping = aes(x = rev(rows), y = cols, fill = counts)) +
  scale_fill_distiller(palette = 'YlGnBu', direction = 1) + 
  geom_tile() +
  theme(aspect.ratio=1, , axis.text.x = element_text(angle = 45, hjust=1))

可以将这一操作应用于列表中的所有矩阵，并将数据集整合成适合用 ggplot2 可视化的“整洁”数据形式。

## Define function to convert a matrix to long format
toLong <- \(x) {
  x |>
    as.table() |>
    as.data.frame() |>
    setNames(c('rows', 'cols', 'counts'))
}

## Apply function to convert all matrices to long format
apas <- list(WT = lapply(loopApaWtHic, toLong),
             FS = lapply(loopApaFsHic, toLong))

## Add loopType to each data.frame and combine
apas <- lapply(apas, \(x) do.call(rbind, Map(cbind, x, loopType = names(x))))

## Add hicMap to each data.frame and combine
apas <- do.call(rbind, Map(cbind, apas, hicMap = names(apas)))

## Reorder factors
apas$loopType <- factor(x = apas$loopType,
                        levels = c("allLoops", "wtLoops", "fsLoops"))
apas$hicMap <- factor(x = apas$hicMap,
                      levels = c("WT", "FS"))

## Visualize with ggplot2
library(ggplot2)
ggplot(data = apas,
       mapping = aes(x = rows, y = cols, fill = counts)) +
  scale_fill_distiller(palette = 'YlGnBu', direction = 1) + 
  facet_grid(hicMap~loopType, scales = "free") +
  geom_tile() +
  theme(aspect.ratio=1, axis.text.x = element_text(angle = 45, hjust=1))

使用 plotgardener 进行可视化

plotgardener 是一个基因组学绘图工具，相比 ggplot2，它提供了更大的灵活性。作为 plotgardener 生态系统的一部分，hictoolsr 提供了一个 plotApa() 函数，该函数与 plotgardener 的其他功能兼容。此外，plotApa() 可以直接作用于矩阵，无需先转换为长格式。

以下是一个使用 plotApa() 和 RColorBrewer 调色板快速可视化的示例：

library(RColorBrewer)

plotApa(apa = loopApaWtHic$allLoops,
        palette = colorRampPalette(brewer.pal(9, 'YlGnBu')))

通过提供位置信息（例如 x、y、宽度、高度等），plotgardener 会切换到多图模式，允许在 pgPage 上进行多种图形排列。可以使用 hictoolsr 和 plotgardener 中的函数来可视化所有 APA 结果：

library(plotgardener)
library(purrr)

## Initiate plotgardener page
pageCreate(width = 4.25, height = 3)

## Define shared parameters
p <- pgParams(x = 0.5,
              y = 0.5,
              width = 1,
              height = 1,
              space = 0.075,
              zrange = c(0, max(unlist(c(loopApaWtHic, loopApaFsHic)))),
              palette = colorRampPalette(brewer.pal(9, 'YlGnBu')))

## Define grid of coordinate positions
xpos <- c(p$x, p$x + p$width + p$space, p$x + (p$width + p$space)*2)
ypos <- c(p$y, p$y + p$height + p$space, p$y + (p$height + p$space)*2)

## Plot row of WT APAs
wt_plots <- 
  pmap(list(loopApaWtHic, xpos, ypos[1]), \(a, x, y) {
    plotApa(params = p, apa = a, x = x, y = y)
  })

## Plot row of FS APAs
fs_plots <- 
  pmap(list(loopApaFsHic, xpos, ypos[2]), \(a, x, y) {
    plotApa(params = p, apa = a, x = x, y = y)
  })

## Add legend
annoHeatmapLegend(plot = wt_plots[[1]],
                  x = p$x + (p$width + p$space)*3,
                  y = ypos[1],
                  width = p$space,
                  height = p$height*0.75,
                  fontcolor = 'black')

## Add text labels
plotText(label = c("All loops", "WT loops", "FS loops"),
         x = xpos + p$width / 2,
         y = ypos[1] - p$space,
         just = c('center', 'bottom'))

plotText(label = c("WT", "FS"),
         x = xpos[1] - p$space,
         y = ypos[1:2] + p$height / 2,
         rot = 90,
         just = c('center', 'bottom'))

## Remove Guides
pageGuideHide()

正如你所看到的，尽管在某些方面 plotgardener 的操作可能较为复杂，但它也带来了更大的灵活性，可以精确控制基因组数据可视化的具体位置和方式。