有了R語言的基礎(chǔ),以及ggplot2繪圖基礎(chǔ)骗奖,我們的生信常用分析圖形的繪制就可以提上日程了醋粟!本系列,師兄就開始帶著大家一起學(xué)習(xí)如何用R語言繪制我們自己的各種分析圖吧重归!
由于本系列的所有分析代碼均為師兄細(xì)心整理和詳細(xì)注釋而成的米愿!歡迎點贊、收藏鼻吮、轉(zhuǎn)發(fā)育苟!
您的支持是我持續(xù)更新的最大動力!
系列內(nèi)容包括:
- 各種類型的熱圖你學(xué)會了嗎椎木?
- 普通熱圖
- 環(huán)形熱圖
- 解鎖火山圖真諦违柏!
- plot函數(shù)就能畫火山圖博烂?
- 高級函數(shù)繪制火山圖--ggplot2、ggpurb
- 經(jīng)典富集分析及氣泡圖漱竖、柱狀圖繪制
- 氣泡圖繪制
- 柱狀圖繪制
- 富集分析圈圖
- 富集分析弦圖
- 繪制一張可以打動審稿人的汕堇椋基圖
- 生存分析 -- KM曲線圖
- 基礎(chǔ)PCA圖
- 云雨圖
- 韋恩圖
- 環(huán)形互作網(wǎng)絡(luò)圖
- 相互作用網(wǎng)絡(luò)圖
- 聚類樹美化
- 富集分析氣泡圖進(jìn)階版
- mantel test相關(guān)性圖
- 詞云圖
- 瀑布圖
- 森林圖
- 曼哈頓圖
- 啞鈴圖
- 三線表
- 嵌套圈圖
- 列線圖
- 蜂群圖
- 箱線圖+貝塞爾曲線
- 矩陣散點圖
- 等等,想到再繼續(xù)補充b扇恰L陕省!
本期富集分析圈圖效果展示
示例數(shù)據(jù)構(gòu)建和展示:
# 繪制富集分析圈圖
library(circlize)
library(grid)
library(graphics)
library(ComplexHeatmap)
data <- read.table("GO_enrichment_stat.txt",header = T)
data_new <- data[,c(2,1)]
# 通路總基因數(shù):min -- 0万矾、 max -- 生物的總基因數(shù)悼吱、rich表示該通路中的基因數(shù);
data_new$gene_num.min <- 0
data_new$gene_num.max <- as.numeric(strsplit(data$BgRatio[1], "/")[[1]][2])
data_new$gene_num.rich <- as.numeric(unlist(lapply(data$BgRatio,
function(x) strsplit(x,"/")[[1]][1])))
# -log10 p值
data_new$"-log10Pvalue" <- -log10(data$pvalue)
# 隨機創(chuàng)建一個上下調(diào)的比例:
ratio <- runif(nrow(data),0,1)
# 根據(jù)這個比例計算上調(diào)和下調(diào)各自占多少:
rich_gene_num <- as.numeric(unlist(lapply(data$GeneRatio,
function(x) strsplit(x,"/")[[1]][1])))
data_new$up.regulated <- round(rich_gene_num*ratio)
data_new$down.regulated <- rich_gene_num - data_new$up.regulated
# 富集分?jǐn)?shù):差異基因中有多少比例在這個通路中:
# 我這里為了讓柱狀圖更清楚良狈,虛構(gòu)了一個ratio后添,不具有實際意義!
data_new$rich.factor <- rich_gene_num/120
colnames(data_new)[c(1,2)] <- c("id", "category")
# 隨機取30個GO作圖薪丁,這里大家可以根據(jù)具體情況選擇合適的通路:
dat <- data_new[sample(247,30),]
dat$id <- factor(rownames(dat), levels = rownames(dat))
# 查看改造后的示例數(shù)據(jù):
> head(dat)
id category gene_num.min gene_num.max gene_num.rich -log10Pvalue up.regulated down.regulated
76 76 BP 0 18866 151 4.038230 6 29
242 242 MF 0 18866 26 3.090989 7 3
34 34 BP 0 18866 138 5.335862 31 5
123 123 BP 0 18866 466 3.397800 45 36
211 211 CC 0 18866 50 3.718427 12 4
29 29 BP 0 18866 22 5.741276 7 5
rich.factor all.regulated up.proportion down.proportion up down
76 0.29166667 35 0.1714286 0.8285714 3234.171 18866
242 0.08333333 10 0.7000000 0.3000000 13206.200 18866
34 0.30000000 36 0.8611111 0.1388889 16245.722 18866
123 0.67500000 81 0.5555556 0.4444444 10481.111 18866
211 0.13333333 16 0.7500000 0.2500000 14149.500 18866
29 0.10000000 12 0.5833333 0.4166667 11005.167 18866
繪圖
# 第一個圈:繪制id
circle_size = unit(1, 'snpc')
circos.par(gap.degree = 0.5, start.degree = 90)
plot_data <- dat[c('id', 'gene_num.min', 'gene_num.max')]
ko_color <- c(rep('#F7CC13',10), rep('#954572',10), rep('#0796E0',9), rep('green', 6)) #各二級分類的顏色和數(shù)目
circos.genomicInitialize(plot_data, plotType = NULL, major.by = 1)
circos.track(
ylim = c(0, 1), track.height = 0.05, bg.border = NA, bg.col = ko_color,
panel.fun = function(x, y) {
ylim = get.cell.meta.data('ycenter')
xlim = get.cell.meta.data('xcenter')
sector.name = get.cell.meta.data('sector.index')
circos.axis(h = 'top', labels.cex = 0.4, labels.niceFacing = FALSE)
circos.text(xlim, ylim, sector.name, cex = 0.4, niceFacing = FALSE)
} )
# 第二圈遇西,繪制富集的基因和富集p值
plot_data <- dat[c('id', 'gene_num.min', 'gene_num.rich', '-log10Pvalue')]
label_data <- dat['gene_num.rich']
p_max <- round(max(dat$'-log10Pvalue')) + 1
colorsChoice <- colorRampPalette(c('#FF906F', '#861D30'))
color_assign <- colorRamp2(breaks = 0:p_max, col = colorsChoice(p_max + 1))
circos.genomicTrackPlotRegion(
plot_data, track.height = 0.08, bg.border = NA, stack = TRUE,
panel.fun = function(region, value, ...) {
circos.genomicRect(region, value, col = color_assign(value[[1]]), border = NA, ...)
ylim = get.cell.meta.data('ycenter')
xlim = label_data[get.cell.meta.data('sector.index'),1] / 2
sector.name = label_data[get.cell.meta.data('sector.index'),1]
circos.text(xlim, ylim, sector.name, cex = 0.4, niceFacing = FALSE)
} )
# 第三圈,繪制上下調(diào)基因
dat$all.regulated <- dat$up.regulated + dat$down.regulated
dat$up.proportion <- dat$up.regulated / dat$all.regulated
dat$down.proportion <- dat$down.regulated / dat$all.regulated
dat$up <- dat$up.proportion * dat$gene_num.max
plot_data_up <- dat[c('id', 'gene_num.min', 'up')]
names(plot_data_up) <- c('id', 'start', 'end')
plot_data_up$type <- 1
dat$down <- dat$down.proportion * dat$gene_num.max + dat$up
plot_data_down <- dat[c('id', 'up', 'down')]
names(plot_data_down) <- c('id', 'start', 'end')
plot_data_down$type <- 2
plot_data <- rbind(plot_data_up, plot_data_down)
label_data <- dat[c('up', 'down', 'up.regulated', 'down.regulated')]
color_assign <- colorRamp2(breaks = c(1, 2), col = c('red', 'blue'))
circos.genomicTrackPlotRegion(
plot_data, track.height = 0.08, bg.border = NA, stack = TRUE,
panel.fun = function(region, value, ...) {
circos.genomicRect(region, value, col = color_assign(value[[1]]), border = NA, ...)
ylim = get.cell.meta.data('cell.bottom.radius') - 0.5
xlim = label_data[get.cell.meta.data('sector.index'),1] / 2
sector.name = label_data[get.cell.meta.data('sector.index'),3]
circos.text(xlim, ylim, sector.name, cex = 0.4, niceFacing = FALSE)
xlim = (label_data[get.cell.meta.data('sector.index'),2]+label_data[get.cell.meta.data('sector.index'),1]) / 2
sector.name = label_data[get.cell.meta.data('sector.index'),4]
circos.text(xlim, ylim, sector.name, cex = 0.4, niceFacing = FALSE)
} )
# 第四圈严嗜,繪制富集因子
plot_data <- dat[c('id', 'gene_num.min', 'gene_num.max', 'rich.factor')]
label_data <- dat['category']
color_assign <- c('BP' = '#F7CC13', 'CC' = '#954572', 'MF' = '#0796E0')#各二級分類的名稱和顏色
circos.genomicTrack(
plot_data, ylim = c(0, 1), track.height = 0.3, bg.col = 'gray95', bg.border = NA,
panel.fun = function(region, value, ...) {
sector.name = get.cell.meta.data('sector.index')
circos.genomicRect(region, value, col = color_assign[label_data[sector.name,1]], border = NA, ytop.column = 1, ybottom = 0, ...)
circos.lines(c(0, max(region)), c(0.5, 0.5), col = 'gray', lwd = 0.3)
} )
category_legend <- Legend(
labels = c('BP', 'CC', 'MF'),#各二級分類的名稱
type = 'points', pch = NA, background = c('#F7CC13', '#954572', '#0796E0'), #各二級分類的顏色
labels_gp = gpar(fontsize = 8), grid_height = unit(0.5, 'cm'), grid_width = unit(0.5, 'cm'))
updown_legend <- Legend(
labels = c('Up-regulated', 'Down-regulated'),
type = 'points', pch = NA, background = c('red', 'blue'),
labels_gp = gpar(fontsize = 8), grid_height = unit(0.5, 'cm'), grid_width = unit(0.5, 'cm'))
pvalue_legend <- Legend(
col_fun = colorRamp2(round(seq(0, p_max, length.out = 6), 0),
colorRampPalette(c('#FF906F', '#861D30'))(6)),
legend_height = unit(3, 'cm'), labels_gp = gpar(fontsize = 8),
title_gp = gpar(fontsize = 9), title_position = 'topleft', title = '-Log10(Pvalue)')
lgd_list_vertical <- packLegend(category_legend, updown_legend, pvalue_legend)
grid.draw(lgd_list_vertical)
circos.clear()
往期優(yōu)秀圖形目錄
以上內(nèi)容僅為群內(nèi)部分內(nèi)容粱檀,不代表全部。詳細(xì)目錄請看下鏈接阻问!