使用GenomicRanges绘制彩色核苷酸替代矩形的序列图方法

实现方案 #

你可以通过将GRanges中的序列拆分到单碱基粒度，再结合ggplot2的geom_text图层实现需求，完整代码如下：

第一步：数据预处理，拆分为单碱基粒度 #

# 加载依赖包
library(tidyverse)
library(GenomicRanges)

# 将GRanges转为数据框，新增read_id用于区分不同的序列行
gr_df <- as.data.frame(gr) %>%
  mutate(read_id = factor(row_number(), levels = rev(row_number())))

# 拆分每个序列为单个碱基，匹配每个碱基对应的基因组坐标
nt_df <- gr_df %>%
  rowwise() %>%
  reframe(
    read_id = read_id,
    abundance = abundance,
    nt = strsplit(seq, "")[[1]],
    pos = start:end
  )

第二步：绘图（两种可选样式） #

样式1：仅展示带颜色区分的碱基

# 定义常规ATCG配色方案，可根据需求自行调整
nt_colors <- c(
  "A" = "#109648",
  "T" = "#255C99",
  "C" = "#F7B32B",
  "G" = "#D62828",
  "N" = "grey50"
)

ggplot(nt_df, aes(x = pos, y = read_id, label = nt, color = nt)) +
  geom_text(size = 3, fontface = "bold") +
  scale_color_manual(values = nt_colors) +
  labs(x = "基因组位置（chr8）", y = "序列 reads", color = "碱基类型") +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )

样式2：保留原有的丰度填充背景，叠加带颜色的碱基

ggplot() +
  # 底层丰度背景层
  geom_tile(
    data = gr_df,
    aes(x = (start + end)/2, y = read_id, width = width, height = 0.8, fill = abundance),
    alpha = 0.3
  ) +
  # 上层碱基层
  geom_text(
    data = nt_df,
    aes(x = pos, y = read_id, label = nt, color = nt),
    size = 3, fontface = "bold"
  ) +
  scale_fill_viridis_c() +
  scale_color_manual(values = nt_colors) +
  labs(x = "基因组位置（chr8）", y = "序列 reads", color = "碱基类型", fill = "丰度") +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )

补充说明 #

如果后续你的数据包含正负链信息，可在拆分序列前调用Biostrings::reverseComplement()对负链序列做反向互补转换，保证碱基展示方向和基因组坐标方向一致。如果序列数量较多，可适当调小geom_text的size参数避免碱基重叠。

内容的提问来源于stack exchange，提问作者Deepak Tanwar