R 数据可视化 —— circlize 基因组初始化

作者: 名本无名 | 来源:发表于2021-06-08 08:50 被阅读0次

前言

圆形可视化广泛应用于基因组及其相关的组学领域中，能够有效地展示高维基因组学数据。

在基因组数据中，通常是根据染色体进行分类，x 轴对应于基因组上的位置，也可以是其他类型的基因组数据

circlize 提供了一些专门的基因组绘图函数，让基因组分析更加简单方便，如：

circos.genomicTrack(): 添加轨迹和图形
circos.genomicPoints(): 添加点
circos.genomicLines(): 添加线条或线段
circos.genomicRect(): 添加矩形
circos.genomicText(): 添加文本
circos.genomicLink(): 添加连接

这样函数与基础的绘制函数是类似的，只是接受的输入数据格式不同，都是基于基础的 circlize 绘图函数实现的（如 circos.track(), circos.points() 等）。

输入数据

基因组数据通常使用的是 BED 格式的文件，即前三列标识某一基因组区域：染色体、起始位置、终止位置。

circlize 提供了一个简单函数 generateRandomBed() 来创建随机的基因组数据。从人类基因组中均匀的生成基因组区域，区域的数量与染色体的大小成正比。

nr 和 nc 参数用于控制需要创建的行列的数量，可能最后生成的行数不一定与 nr 相同，fun 参数用于接受自定义的值生成函数。

> bed <- generateRandomBed()
> head(bed)
   chr   start     end      value1
1 chr1  261327  520533  0.07617057
2 chr1  596180  606938  0.81289852
3 chr1  769058 1176608  0.61876561
4 chr1 1179719 1671784  0.21739949
5 chr1 1860787 2066114 -0.01665364
6 chr1 2183578 2277911  0.01477448
> # 设置行列数量
> bed <- generateRandomBed(nr = 200, nc = 4)
> nrow(bed)
[1] 205
> # 自定义值生成函数
> bed <- generateRandomBed(
+   nc = 2, fun = function(k) sample(letters, k, replace = TRUE)
+   )
> head(bed)
   chr   start     end value1 value2
1 chr1  154420  432520      o      q
2 chr1  621080  658294      w      g
3 chr1  923320  962390      b      t
4 chr1  964699 1202322      y      v
5 chr1 1336707 1405512      r      g
6 chr1 1455202 1534223      i      a

初始化

1. 染色体条带初始化

cytoband 类型的数据是理想的输入格式，其包含染色体的长度以及染色体条带信息，能够有效标识染色体的位置

1.1 基本使用

如果是绘制人类基因组数据，可以直接使用 circos.initializeWithIdeogram() 函数，例如

circos.initializeWithIdeogram()
text(0, 0, "default", cex = 1)
circos.clear()

染色体名称显示的是纯数字，但是其内部的索引名称还是带有 chr 的字符串

> circos.info()
All your sectors:
 [1] "chr1"  "chr2"  "chr3"  "chr4"  "chr5"  "chr6"  "chr7"  "chr8"  "chr9"  "chr10" "chr11" "chr12" "chr13" "chr14"
[15] "chr15" "chr16" "chr17" "chr18" "chr19" "chr20" "chr21" "chr22" "chrX"  "chrY" 

All your tracks:
[1] 1 2

Your current sector.index is chrY
Your current track.index is 2

circos.initializeWithIdeogram() 默认使用的是 hg19 版的 cytoband 数据，可以使用 species 参数来指定为 hg18 或其他物种

circos.initializeWithIdeogram(species = "hg18")
circos.initializeWithIdeogram(species = "mm10")

会自动从网上下载对应的数据，如果提供的物种不存在，还是会从 UCSC 数据库中下载染色体信息 chromInfo 文件，有染色体长度，但是不包含条带信息

如果网络受限或者 UCSC 中没有对应的数据，则可以自己手动创建数据框，或传入存储在本地的文件

cytoband.file <- system.file(
  package = "circlize", "extdata", "cytoBand.txt"
)
circos.initializeWithIdeogram(cytoband.file)

cytoband.df <- read.table(
  cytoband.file, colClasses = c(
    "character", "numeric", "numeric", "character", "character"), 
  sep = "\t"
)
circos.initializeWithIdeogram(cytoband.df)

注意，如果是从文件中读取，需要指定每列的数据类型，并指定位置列为 numeric 类型，因为 read.table 会将数字作为 integer 类型，而基因组长度是较大的数，会导致整数溢出

circos.intializeWithIdeogram() 默认会读取 cytoband 数据中的所有信息，使用 chromosome.index 参数可以选择要显示的染色体

circos.initializeWithIdeogram(
  chromosome.index = paste0("chr", c(3,5,2,8))
)
text(0, 0, "subset of chromosomes", cex = 1)
circos.clear()

如果没有相应的物种，且使用了 chromInfo 文件时，会包含很多的短的 contigs，使用 chromosome.index 可以删除一些不需要的 contigs

1.2 预定义轨迹

使用 circos.initializeWithIdeogram() 初始化圆形图，会创建两个轨迹，一个用于包含轴和染色体名称，另一个轨迹用于绘制 ideogram，使用 plotType 可以控制需要绘制的轨迹

par(mfcol = c(1, 2))
circos.initializeWithIdeogram(plotType = c("axis", "labels"))
text(0, 0, "plotType = c('axis', 'labels')", cex = 1)
circos.clear()

circos.initializeWithIdeogram(plotType = NULL)
text(0, 0, "plotType = NULL", cex = 1)
circos.clear()

1.3 其他设置

与常规的圆形图类似，可以使用 circos.par() 函数来控制圆形布局

par(mfcol = c(1, 2))
circos.par("start.degree" = 90)
circos.initializeWithIdeogram()
circos.clear()
text(0, 0, "'start.degree' = 90", cex = 1)

circos.par("gap.degree" = rep(c(2, 4), 12))
circos.initializeWithIdeogram()
circos.clear()
text(0, 0, "'gap.degree' = rep(c(2, 4), 12)", cex = 1)

2. 自定义染色体轨迹

circos.initializeWithIdeogram() 函数默认会初始化圆形布局，并添加两个轨迹，通过设置 plotType = NULL，只创建布局而不添加轨迹，我们就可以添加自定义图形样式

例如，我们为染色体设置不同的颜色，并将染色体名称放置在单元格内部

circos.initializeWithIdeogram(plotType = NULL)
circos.track(
  ylim = c(0, 1), 
  panel.fun = function(x, y) {
    chr = CELL_META$sector.index
    xlim = CELL_META$xlim
    ylim = CELL_META$ylim
    circos.rect(xlim[1], 0, xlim[2], 1, col = rand_color(1))
    circos.text(
      mean(xlim), mean(ylim), chr, cex = 0.7,
      col = "white", facing = "inside",
      niceFacing = TRUE
    )
  }, 
  track.height = 0.15, bg.border = NA
)
circos.clear()

3. 常规基因组类别初始化

染色体只是一种特殊的基因组分类，使用 circos.genomicInitialize() 可以初始化任意基因组分类的圆形布局，circos.initializeWithIdeogram() 函数也是基于 circos.genomicInitialize() 实现的。

circos.genomicInitialize() 函数也是接受数据框型的输入数据，数据必须至少包含三列，第一列代表基因组分类，后两列为每种分类在基因组中的顺序

例如，基因的位置信息

df <- data.frame(
  name  = c("TP53",  "TP63",    "TP73"),
  start = c(7565097, 189349205, 3569084),
  end   = c(7590856, 189615068, 3652765))
circos.genomicInitialize(df)

并不是说一个基因只能记录为一行，也可以是多行，比如基因的转录本信息。我们读取包中自带的示例数据，包含 TP53、TP63、TP73 三个基因的信息

> tp_family <- readRDS(system.file(
+   package = "circlize", "extdata", "tp_family_df.rds")
+   )
> head(tp_family)
  gene   start     end        transcript exon
1 TP53 7565097 7565332 ENST00000413465.2    7
2 TP53 7577499 7577608 ENST00000413465.2    6
3 TP53 7578177 7578289 ENST00000413465.2    5
4 TP53 7578371 7578554 ENST00000413465.2    4
5 TP53 7579312 7579590 ENST00000413465.2    3
6 TP53 7579700 7579721 ENST00000413465.2    2

绘制填充色来标识三个基因

circos.genomicInitialize(tp_family)
circos.track(
  ylim = c(0, 1), 
  bg.col = c("#1f78b480", "#33a02c80", "#e31a1c80"), 
  bg.border = NA, track.height = 0.05
)

绘制基因的转录本

n <- max(tapply(
  tp_family$transcript, tp_family$gene, 
  function(x) length(unique(x)))
)
circos.genomicTrack(
  tp_family, ylim = c(0.5, n + 0.5), 
  panel.fun = function(region, value, ...) {
    all_tx = unique(value$transcript)
    for(i in seq_along(all_tx)) {
      l = value$transcript == all_tx[i]
      # 对于每个转录本
      current_tx_start = min(region[l, 1])
      current_tx_end = max(region[l, 2])
      circos.lines(
        c(current_tx_start, current_tx_end),
        c(n - i + 1, n - i + 1), col = "#CCCCCC"
      )
      circos.genomicRect(
        region[l, , drop = FALSE],
        ytop = n - i + 1 + 0.4,
        ybottom = n - i + 1 - 0.4,
        col = "orange",
        border = NA
      )
    }
  }, 
  bg.border = NA, track.height = 0.4
)
circos.clear()

4. 放大基因组

基因组区域的放大方式类似于前面介绍的，也是将需要放大的区域的数据提取出来，并设置不同的分类，然后添加到输入数据中

extend_chromosomes <- function(bed, chromosome, prefix = "zoom_") {
  zoom_bed = bed[bed[[1]] %in% chromosome, , drop = FALSE]
  zoom_bed[[1]] = paste0(prefix, zoom_bed[[1]])
  rbind(bed, zoom_bed)
}

我们使用 read.cytoband() 函数从 UCSC 中下载并读取 cytoband 数据

cytoband <- read.cytoband()
cytoband_df <- cytoband$df
chromosome <- cytoband$chromosome

xrange <- c(cytoband$chr.len, cytoband$chr.len[c("chr1", "chr2")])
normal_chr_index <- 1:24
zoomed_chr_index <- 25:26

# 设置宽度
sector.width <- c(
  xrange[normal_chr_index] / sum(xrange[normal_chr_index]), 
  xrange[zoomed_chr_index] / sum(xrange[zoomed_chr_index])
)

绘制图形

circos.par(start.degree = 90)
circos.initializeWithIdeogram(
  extend_chromosomes(cytoband_df, c("chr1", "chr2")), 
  sector.width = sector.width)

添加一个新的轨迹

bed <- generateRandomBed(500)
circos.genomicTrack(
  extend_chromosomes(bed, c("chr1", "chr2")),
  panel.fun = function(region, value, ...) {
    circos.genomicPoints(
      region, value, pch = 16, 
      cex = 0.3, col = 'blue')
    }
)

添加连接

circos.link(
  "chr1", get.cell.meta.data("cell.xlim", sector.index = "chr1"),
  "zoom_chr1", get.cell.meta.data("cell.xlim", sector.index = "zoom_chr1"),
  col = "#33a02c20", border = NA)
circos.clear()

5. 合并两个基因组

在某些情况下，可能想要将多个基因组绘制在同一个圆形图中，例如，我们可以将人类与小鼠的基因组组合起来

首先，获取两个基因组的 cytoband 数据

human_cytoband <- read.cytoband(species = "hg19")$df
mouse_cytoband <- read.cytoband(species = "mm10")$df

注意，要区分两个基因组的染色体名称，给它们加上一个前缀

human_cytoband[ ,1] <- paste0("human_", human_cytoband[, 1])
mouse_cytoband[ ,1] <- paste0("mouse_", mouse_cytoband[, 1])

然后，合并两个数据

cytoband <- rbind(human_cytoband, mouse_cytoband)

> head(cytoband)
          V1       V2       V3     V4     V5
1 human_chr1        0  2300000 p36.33   gneg
2 human_chr1  2300000  5400000 p36.32 gpos25
3 human_chr1  5400000  7200000 p36.31   gneg
4 human_chr1  7200000  9200000 p36.23 gpos25
5 human_chr1  9200000 12700000 p36.22   gneg
6 human_chr1 12700000 16200000 p36.21 gpos50

通过设置 chromosome.index 参数的值，让两个基因组的 1 号染色体放置在相邻的位置

chromosome.index <- c(
  paste0("human_chr", c(1:22, "X", "Y")), 
  rev(paste0("mouse_chr", c(1:19, "X", "Y")))
)
circos.initializeWithIdeogram(
  cytoband, chromosome.index = chromosome.index
)
circos.clear()

染色体数据太多了，使得图像比较臃肿，为了让图形更加美观，我们关闭染色体名称和轴刻度的显示，只简单的显示染色体号，并添加染色体分组颜色和间距

# 两组之间 5 度间距
circos.par(
  gap.after = c(rep(1, 23), 5, rep(1, 20), 5)
)
circos.initializeWithIdeogram(
  cytoband, plotType = NULL, 
  # 染色体顺序
  chromosome.index = chromosome.index
)
# 添加染色体号
circos.track(
  ylim = c(0, 1), track.height = mm_h(1), 
  panel.fun = function(x, y) {
    circos.text(
      CELL_META$xcenter,
      CELL_META$ylim[2] + mm_y(2),
      gsub(".*chr", "", CELL_META$sector.index),
      cex = 0.6,
      niceFacing = TRUE
    )
  }, 
  cell.padding = c(0, 0, 0, 0), bg.border = NA
)

highlight.chromosome(
  paste0("human_chr", c(1:22, "X", "Y")), 
  col = "#66c2a5", track.index = 1
)
highlight.chromosome(
  paste0("mouse_chr", c(1:19, "X", "Y")), 
  col = "#fc8d62", track.index = 1
)
# 绘制 ideogram
circos.genomicIdeogram(cytoband)
circos.clear()

同样地，对于不包含条带信息，只有染色体范围的输入数据，我们也可以进行组合。

首先，使用 read.chromInfo() 来获取染色体范围信息，然后合并两份数据

human_chromInfo <- read.chromInfo(species = "hg19")$df
mouse_chromInfo <- read.chromInfo(species = "mm10")$df
human_chromInfo[ ,1] <- paste0("human_", human_chromInfo[, 1])
mouse_chromInfo[ ,1] <- paste0("mouse_", mouse_chromInfo[, 1])
chromInfo <- rbind(human_chromInfo, mouse_chromInfo)
# 控制染色体的顺序
chromInfo[, 1] <- factor(chromInfo[ ,1], levels = chromosome.index)

> head(chromInfo)
         chr start       end
1 human_chr1     0 249250621
2 human_chr2     0 243199373
3 human_chr3     0 198022430
4 human_chr4     0 191154276
5 human_chr5     0 180915260
6 human_chr6     0 171115067

使用 genomicInitialize() 来初始化布局，并添加图形

circos.par(gap.after = c(rep(1, 23), 5, rep(1, 20), 5))
circos.genomicInitialize(chromInfo, plotType = NULL)
circos.track(
  ylim = c(0, 1),
  panel.fun = function(x, y) {
    circos.text(
      CELL_META$xcenter,
      CELL_META$ylim[2] + mm_y(2),
      gsub(".*chr", "", CELL_META$sector.index),
      cex = 0.6,
      niceFacing = TRUE
    )
  },
  track.height = mm_h(1),
  cell.padding = c(0, 0, 0, 0),
  bg.border = NA
)
highlight.chromosome(
  paste0("human_chr", c(1:22, "X", "Y")), 
  col = "#66c2a5", track.index = 1
)
highlight.chromosome(
  paste0("mouse_chr", c(1:19, "X", "Y")), 
  col = "#fc8d62", track.index = 1
)
# 添加空白轨迹
circos.track(ylim = c(0, 1))
circos.clear()

我们可以为图形添加更多的轨迹以及连接，让图形看起来更加的充实

# 设置间距
circos.par(gap.after = c(rep(1, 23), 5, rep(1, 20), 5))
# 初始化布局，不添加图形
circos.genomicInitialize(chromInfo, plotType = NULL)
# 添加数字染色体号
circos.track(
  ylim = c(0, 1),
  panel.fun = function(x, y) {
    circos.text(
      CELL_META$xcenter,
      CELL_META$ylim[2] + mm_y(2),
      gsub(".*chr", "", CELL_META$sector.index),
      cex = 0.6,
      niceFacing = TRUE
    )
  },
  track.height = mm_h(1),
  cell.padding = c(0, 0, 0, 0),
  bg.border = NA
)
# 添加分组颜色轨迹
highlight.chromosome(
  paste0("human_chr", c(1:22, "X", "Y")),
  col = "#66c2a5", track.index = 1
)
highlight.chromosome(
  paste0("mouse_chr", c(1:19, "X", "Y")),
  col = "#fc8d62", track.index = 1
)

# 添加 ideogram
circos.genomicIdeogram(cytoband)

# 创建随机数据
human_df <- generateRandomBed(200, species = "hg19")
mouse_df <- generateRandomBed(200, species = "mm10")
human_df[, 1] <- paste0("human_", human_df[, 1])
mouse_df[, 1] <- paste0("mouse_", mouse_df[, 1])
df <- rbind(human_df, mouse_df)
# 添加点图
circos.genomicTrack(
  df,
  panel.fun = function(region, value, ...) {
    circos.genomicPoints(region, value, col = rand_color(1), cex = 0.5, ...)
  }
)

# 添加人类与小鼠基因组之间的连接
human_mid <- data.frame(
  chr = paste0("human_chr", 1:19),
  mid = round((human_chromInfo[1:19, 2] + human_chromInfo[1:19, 3]) / 2)
)
mouse_mid <- data.frame(
  chr = paste0("mouse_chr", 1:19),
  mid = round((mouse_chromInfo[1:19, 2] + mouse_chromInfo[1:19, 3]) / 2)
)
circos.genomicLink(human_mid, mouse_mid, col = rand_color(19))
circos.clear()
# 添加注释
text(-0.9,-0.8, "Human\ngenome")
text(0.9, 0.8, "Mouse\ngenome")

完整代码：https://github.com/dxsbiocc/learn/blob/main/R/plot/genomic_human_mouse.R

R 数据可视化 —— circlize 基因组初始化

前言

输入数据

初始化

1. 染色体条带初始化

1.1 基本使用

1.2 预定义轨迹

1.3 其他设置

2. 自定义染色体轨迹

3. 常规基因组类别初始化

4. 放大基因组

5. 合并两个基因组

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

R语言绘图

R生信相关