多基因的荧光定量pcr数据处理（R数据可视化）
多基因的荧光定量pcr数据处理（R数据可视化）如果是绘制人类基因组数据，可以直接使用 circos.initializeWithIdeogram() 函数，例如> 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(n
前言圆形可视化广泛应用于基因组及其相关的组学领域中，能够有效地展示高维基因组学数据。
在基因组数据中，通常是根据染色体进行分类，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()
与常规的圆形图类似，可以使用 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)
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()
染色体只是一种特殊的基因组分类，使用 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()
基因组区域的放大方式类似于前面介绍的，也是将需要放大的区域的数据提取出来，并设置不同的分类，然后添加到输入数据中
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()
在某些情况下，可能想要将多个基因组绘制在同一个圆形图中，例如，我们可以将人类与小鼠的基因组组合起来
首先，获取两个基因组的 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
