circlize 之可视化基因组数据

点击上方公众号查看更多精彩内容

circlize 可用来在基因组学研究领域可视化及展示基因组之间不同染色体的关联性。

为了便于基因组学分析，circlize 特别提供了专注于绘制基因组图的功能。这些函数与基本图形函数是类似的，但是需要输入特殊格式的数据：

输入数据：

基因组数据通常存储为一个 table 形式，其中前三列定义基因组区域，随后的列是与相应区域相关联的值。每个基因组区域由三个要素组成:基因组范畴(多数为染色体)、基因组范畴上的起始位置和结束位置。这种数据结构被称为 BED 格式，广泛应用于基因组研究。

circlize 提供了一个简单的函数 generateRandomBed()，它生成随机的基因组数据。位置是由人类基因组均匀地产生的，染色体上的区域数与染色体的长度近似成正比。其中，nr 和 nc 参数控制用户需要的行数和列数。请注意，nr 与函数返回的行数不完全相同。Fun 参数是用来生成随机值的自定义函数。

set.seed(999)
bed = generateRandomBed()
head(bed)
##    chr   start     end       value1
## 1 chr1  118750  451929 -0.367702502
## 2 chr1  472114  805024 -0.001764946
## 3 chr1  807013  914103 -0.668379255
## 4 chr1 1058915 1081590  0.390291705
## 5 chr1 1194426 1341960 -1.305901261
## 6 chr1 1670402 2048723  0.340227447

生成特定行数和列数：

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  275067  349212      q      e
## 2 chr1  350892  658970      u      f
## 3 chr1  674620  875538      y      r
## 4 chr1 1053255 1115056      p      e
## 5 chr1 1127139 1545066      a      s
## 6 chr1 1864092 1973553      x      u

1、cytoband 基因组数据初始化

cytoband 数据是初始化基因组图的理想数据源。它包含染色体的长度以及所谓的染色体带注释，以帮助识别染色体上的位置。

使用 circos.initializeWithIdeogram()函数加载 cytoband 基因组数据，默认加载的是 hg19 的基因组数据：

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

查看信息：

circos.info()
## All your sectors:
##  [1] "chr1"  "chr2"  "chr3"  "chr4"  "chr5"  "chr6"  "chr7"  "chr8"  "chr9"
## [10] "chr10" "chr11" "chr12" "chr13" "chr14" "chr15" "chr16" "chr17" "chr18"
## [19] "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.clear()

也可以通过指定物种参数初始化其他物种，并自动下载对应物种的 cytoband 数据文件：

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

当你在处理稀有物种和没有可用的 cytoband 数据时，circos.initializeWithIdeogram()将尝试继续下载 UCSC 的 chromInfo 文件，其中也包含染色体的长度，但是当然，在图上没有 ideogram track。

在某些情况下，当没有互联网连接或 UCSC 上没有相应的数据可用的。可以手动构造一个包含染色体范围的 table，如果它存储在一个文件中或目录里，并发送到 circos.initializeWithIdeogram()：

# 放在目录下
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)

默认情况下，circos.initiizewithideogram()使用数据中可用的所有染色体来初始化环形图。用户可以通过指定 chromosome.index 来选择染色体的一个子集。这个参数也适用于染色体排序：

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

在初始化环形图之后，circos.initializeWithIdeogram()另外创建一个轨道，其中有基因组轴和染色体名称，并在有 ideogram 的地方创建另一个轨道(取决于 cytoband 数据是否可用)。plotType 参数用于控制添加哪种类型的轨道：

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()

我们可以用 circos.par()函数来设置其它参数，需要 circos.clear()来完成绘图：

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 时，只初始化环形布局，但不添加任何内容。这使得用户可以完全设计自己风格的染色体轨道。

在下面的例子中，我们使用不同的颜色来代表染色体，并将染色体的名称放在每个单元格的中心：

set.seed(123)
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

创建 3 个轨道：

circos.genomicInitialize(tp_family)
circos.track(ylim = c(0, 1),
    bg.col = c("#FF000040", "#00FF0040", "#0000FF40"),
    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]
            # for each transcript
            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()

你可能会注意到坐标轴的开始变成了 0KB，而不是原来的值。这只是一个轴标签的调整，以反映到每个基因开始的相对距离，而细胞中的坐标仍然使用原始值。设置 tickLabelsStartFromZero 为 FALSE 将坐标轴标签恢复为原始值。

4、缩放染色体

我们首先定义一个函数 extend.chromosome()，它将染色体子集中的数据复制到原始数据中：

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 数据。接下来，正常染色体和放大染色体的 x 范围分别归一化：

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

# normalize in normal chromsomes and zoomed chromosomes separately
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)
})

添加连接：

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 = "#00000020", border = NA)
circos.clear()

5、整合两个基因组

在某些情况下，希望在环形图中可视化多个基因组。这可以通过合成一个组合基因组来实现。在下面的例子中，结合了人类和老鼠的基因组：

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

排序，使让人类 1 号染色体接近老鼠 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()

关掉染色体名称和轴，新建一个小轨道来区分人类染色体和老鼠染色体用 highlight.chromosome()函数，对于染色体，我只写数值索引(也包括 X 和 Y)，人类和老鼠的染色体间隔 5。（circos.par("gap.after")）：

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), 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 = "red", track.index = 1)
highlight.chromosome(paste0("mouse_chr", c(1:19, "X", "Y")),
    col = "blue", track.index = 1)

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)
# note the levels of the factor controls the chromosome orders in the plot
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

对于指定的染色体：范围，使用 circos.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 = "red", track.index = 1)
highlight.chromosome(paste0("mouse_chr", c(1:19, "X", "Y")),
    col = "blue", 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 = "red", track.index = 1)
highlight.chromosome(paste0("mouse_chr", c(1:19, "X", "Y")),
    col = "blue", track.index = 1)

circos.genomicIdeogram(cytoband)

# a track of points
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, ...)
})

# links between human and mouse genomes
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, "Humanngenome")
text(0.9, 0.8, "Mousengenome")