karyoploteR：基于 BigWig 绘制多层信号图 + 自定义基因组支持

本教程将以 食蟹猕猴（Macaca fascicularis）的 snATAC-seq 测序数据中提取的基因组可及性信号（以 BigWig 格式 存储的 Peak 强度数据）作为基础输入，借助 karyoploteR 包构建一个定制化的 Genome Track 图谱，其中叠加了 SNP 注释层 和 PeakLink 互作层。

以下是本教程创建的信号图示例：

📌 本教程将涵盖以下关键知识点：

1. 自定义基因组配置的构建
   学会如何使用自定义的参考基因组（如从 UCSC 下载的 GTF 注释）在 karyoploteR 中创建适用于各类生物（如食蟹猕猴）的基因组信息，包含染色体结构和基因注释。

2. BigWig 信号的平滑可视化
   掌握如何替代 kpPlotBigWig 函数，通过 kpArea 等方式绘制经过平滑处理的 ATAC-seq 等 bigWig 信号区域，实现更柔和、直观的图形呈现。

3. 深入理解并运用 karyoploteR 绘图语法
   熟悉 karyoploteR 中各类函数（如 plotKaryotype, kpPlotGenes, kpAddLabels, kpPlotLinks, kpAxis 等）的基本用法与参数设置，构建结构层次分明、内容丰富的基因组轨迹图谱（Genome Track Plot）。

教程依赖软件包：Signac, dplyr, karyoploteR, rtracklayer, plyranges

1、创建基因组注释

在 karyoploteR 中绘制包含 基因注释 的基因组轨迹图时，依赖的是 TxDb 格式 的注释对象。由于我们通常从公共数据库（如 Ensembl 或 NCBI）下载的是 GTF 文件，因此需首先通过 GenomicFeatures::makeTxDbFromGFF 函数将 GTF 文件转换为 TxDb 对象。

然而需要注意的是，TxDb 并不保留 GTF 文件中的 gene_name 字段，因此如果我们希望在图中显示更直观的基因名称（而非仅有的 gene_id），就需要单独提取并构建一张 gene_id 与 gene_name 的映射表，作为后续注释和可视化的补充。

🧾 代码详解：

### 从GTF 创建 TxDb 注释对象
macfas5.txdb <- GenomicFeatures::makeTxDbFromGFF("Macaca_fascicularis_5.0.91.2.2_withchrname.gtf", format = "gtf")

### 从GTF 创建 gene_id - gene_name 对照表
macfas5.gtf <- rtracklayer::import("Macaca_fascicularis_5.0.91.2.2_withchrname.gtf")
macfas5.gene_map <- data.frame(
    gene_id = mcols(macfas5.gtf)$gene_id,
    gene_name = mcols(macfas5.gtf)$gene_name
) %>% distinct(gene_id,.keep_all = T) %>% 
    mutate(gene_name = if_else(is.na(gene_name), gene_id, gene_name)) %>%
    tibble::remove_rownames() %>% tibble::column_to_rownames("gene_id")

macfas5.gene_map %>% head
# -------------------------------------------------------------------------
gene_name
ENSMFAG00000044637              PGBD2
ENSMFAG00000011984                 U6
ENSMFAG00000039056             ZNF692
ENSMFAG00000030010             ZNF672
ENSMFAG00000002737            SH3BP5L
ENSMFAG00000005691 ENSMFAG00000005691

2、示例轨道数据创建

我们将构建一组用于示范的 轨道可视化数据集，
本示例主要用到 SNP数据、PeakLinks数据和Peaks数据三类数据：

1. 创建并扩展一个基因组区域窗口，例如选取第4号染色体的某个区段

mark.regs <- StringToGRanges("chr4-65401801-65402311") 
### 将该区域向左右各扩展 250 kb，以便形成一个更大的展示窗口用于可视化上下游特征。
zoom.wind <- Extend(mark.regs, 250000, 250000) %>% suppressWarnings() ### 
# -------------------------------------------------------------------------
> zoom.wind
GRanges object with 1 range and 0 metadata columns:
      seqnames            ranges strand
                    
  [1]     chr4 65151801-65652311      *
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

2. SNP 数据：模拟分布于指定基因组区域内的若干个单核苷酸多态性位点（SNP），并为每个位点分配对应的 rsID 与统计显著性 P 值，便于后续在图谱中添加显著性标注或颜色分层；

snp_data <- createRandomRegions(genome = zoom.wind, 
                                nregions = 100000,
                                length.mean = 1,
                                length.sd = 0) %>% 
    mutate(rsid = paste0("rs", sample(1e7:1e8, n(), replace = FALSE)),
           pval = runif(n(), min = 1e-10, max = 1),
           log.pval = -log10(pval)) 
# -------------------------------------------------------------------------
> snp_data
GRanges object with 100000 ranges and 3 metadata columns:
           seqnames    ranges strand |        rsid      pval  log.pval
                  |   
       [1]     chr4  55511005      * |  rs47333361 0.1209773  0.917296
       [2]     chr4  45814413      * |  rs24937342 0.4275664  0.368996

3. Peaks 数据：生成一组不重叠的基因组区域，模拟 snATAC-seq 检测到的开放染色质区域（peaks），作为信号轨迹（如 BigWig）之外的结构性标注进行叠加展示。

random_peaks <- createRandomRegions(genome = zoom.wind, nregions = 10000,
                                    length.mean = 500,length.sd = 0) %>% 
    subsetByOverlaps(., zoom.wind)
# -------------------------------------------------------------------------
> random_peaks
GRanges object with 83 ranges and 0 metadata columns:
       seqnames            ranges strand
                     
   [1]     chr4 65422201-65422700      *
   [2]     chr4 65620347-65620846      *

4. PeakLinks 数据：构建一组模拟的 chromatin interaction（染色质互作）信息，例如从 ATAC-seq 或 Hi-C 数据中推断得到的 peak-to-peak 或 peak-to-gene 连接关系。

loops.data <- data.frame(
    site.A = gsub(":","-",as.character(sample(random_peaks, 100,replace = T))),
    site.B = gsub(":","-",as.character(sample(random_peaks, 100,replace = T)))
) %>% 
    tidyr::separate(site.A,into = c("chr", "start", "end"),sep = "-") %>% 
    tidyr::separate(site.B,into = c("link.chr", "link.start", "link.end"),sep = "-") %>% 
    mutate(link.strand = "*") %>% 
    makeGRangesFromDataFrame(keep.extra.columns = T)
# -------------------------------------------------------------------------
> loops.data
GRanges object with 100 ranges and 4 metadata columns:
        seqnames            ranges strand |    link.chr  link.start    link.end link.strand
                       |    
    [1]     chr4 65282464-65282963      * |        chr4    65205043    65205542           *
    [2]     chr4 65369837-65370336      * |        chr4    65299425    65299924           *

📌说明：以上是示例从随机生成的Peaks 数据 random_peaks 中随机抽取 100 个 peak 对，分别作为互作的两端 site.A 和 site.B，根据 karyoploteR 中的说明规范化成 GRanges 格式，其中 site.A 分成 chr/start/end，site.B 分成 link.chr/link.start/link.end，并添加 link.strand 列，设为 "*" 表示不区分链方向，然后将site.A 的区域作为 GRanges site.B的位点拆分信息作为附加列 （即 loop 的另一端信息），形成染色质互作（loops）示例数据集。

3、对BigWig信号进行平滑处理

有时候我们获取的（比如从别人的数据集里面下载）原始 snATAC-seq bigWig 信号 可能会是在 单碱基或较小窗口 下生成的，因而常伴随较大锯齿波动。在较大尺度的基因组浏览窗口中，这类高分辨率信号往往不利于直观展示其整体变化趋势。因此，我们需要通过灵活滑动窗口（如移动平均）对信号进行平滑处理，以降低局部噪音、突出区域性趋势，从而提升可视化效果，便于识别潜在的调控区域并辅助功能解读。

🧾 代码详解：

bw.file <- ".bw" ### 选择导入 ATAC-seq 特征的 bigwig 文件
{
    k.used = 100 #设置滑动窗口大小，用于对信号进行平滑处理。表示在 zoom 区间内取前后共 k 个点进行平均。
    
    bw_smoothed <- ggcoverage::LoadTrackFile(track.file = bw.file, format = "bw",
                                             region = as.character(zoom.wind), 
                                             extend = 0) %>%
        group_by(seqnames) %>% # 按染色体分组，避免跨染色体做平滑处理
        # 对信号强度 (score) 进行滑动平均（rollmean）平滑处理
        # 使用 zoo 包的函数，fill = NA 表示边缘不足 k 的点将填 NA
        mutate(score = zoo::rollmean(score, k = k.used, fill = NA)) %>%
        na.omit() %>% makeGRangesFromDataFrame(keep.extra.columns = T)
}
# -------------------------------------------------------------------------
> bw_smoothed
GRanges object with 18300 ranges and 1 metadata column:
          seqnames            ranges strand |       score
                         |   
      [1]     chr4 65154571-65154580      * | 0.000681153
      [2]     chr4 65154581-65154610      * | 0.000705053

以下展示在250kb窗口下，平滑前后的信号图视觉效果对比 (添加了k = 250 的平滑)：

# 直接绘制bigWig信号-------------------------------------------------------------------------
kp <- plotKaryotype(zoom = zoom.wind, plot.type = 5, cex=1)
kpPlotBigWig(kp, data = bw.file, ymax="visible.region" )

# 绘制平滑后的 bigWig信号-------------------------------------------------------------------------
kp <- plotKaryotype(zoom = zoom.wind, plot.type = 5, cex=1)
kpArea(kp, data = bw_smoothed, y = bw_smoothed$score, 
       ymin = min(bw_smoothed$score), ymax = max(bw_smoothed$score))

4、将所有示例数据整合创建综合基因组轨迹图谱

🧾 代码详解：