美文网首页生信常用分析图形绘制
生信常用分析图形绘制03 -- 富集分析圈图

生信常用分析图形绘制03 -- 富集分析圈图

作者: 生信师兄 | 来源:发表于2022-08-13 00:16 被阅读0次
    封面

    有了R语言的基础,以及ggplot2绘图基础,我们的生信常用分析图形的绘制就可以提上日程了!本系列,师兄就开始带着大家一起学习如何用R语言绘制我们自己的各种分析图吧!

    由于本系列的所有分析代码均为师兄细心整理和详细注释而成的!欢迎点赞、收藏、转发!

    您的支持是我持续更新的最大动力!

    示例数据和代码获取

    系列内容包括:

    • 各种类型的热图你学会了吗?
      • 普通热图
      • 环形热图
    • 解锁火山图真谛!
      • plot函数就能画火山图?
      • 高级函数绘制火山图--ggplot2、ggpurb
    • 经典富集分析及气泡图、柱状图绘制
      • 气泡图绘制
      • 柱状图绘制
    • 富集分析圈图
    • 富集分析弦图
    • 绘制一张可以打动审稿人的桑基图
    • 生存分析 -- KM曲线图
    • 基础PCA图
    • 云雨图
    • 韦恩图
    • 环形互作网络图
    • 相互作用网络图
    • 聚类树美化
    • 富集分析气泡图进阶版
    • mantel test相关性图
    • 词云图
    • 瀑布图
    • 森林图
    • 曼哈顿图
    • 哑铃图
    • 三线表
    • 嵌套圈图
    • 列线图
    • 蜂群图
    • 箱线图+贝塞尔曲线
    • 矩阵散点图
    • 等等,想到再继续补充!!!

    本期富集分析圈图效果展示

    image.png

    示例数据构建和展示:

    # 绘制富集分析圈图
    library(circlize)
    library(grid)
    library(graphics)
    library(ComplexHeatmap)
    
    data <- read.table("GO_enrichment_stat.txt",header = T)
    
    data_new <- data[,c(2,1)]
    
    # 通路总基因数:min -- 0、 max -- 生物的总基因数、rich表示该通路中的基因数;
    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)
    
    # 随机创建一个上下调的比例:
    ratio <- runif(nrow(data),0,1)
    # 根据这个比例计算上调和下调各自占多少:
    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
    
    # 富集分数:差异基因中有多少比例在这个通路中:
    # 我这里为了让柱状图更清楚,虚构了一个ratio,不具有实际意义!
    data_new$rich.factor <- rich_gene_num/120
    colnames(data_new)[c(1,2)] <- c("id", "category")
    
    # 随机取30个GO作图,这里大家可以根据具体情况选择合适的通路:
    dat <- data_new[sample(247,30),]
    
    dat$id <- factor(rownames(dat), levels = rownames(dat))
    
    # 查看改造后的示例数据:
    > 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)) #各二级分类的颜色和数目
    
    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)  
      } )
    
    fig1
    # 第二圈,绘制富集的基因和富集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)  
      } )
    
    # 第三圈,绘制上下调基因
    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)  
      } )
    
    fig2
    # 第四圈,绘制富集因子
    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()
    
    fig3

    往期优秀图形目录

    渐变火山图 气泡图+相关性热图
    复杂提琴图
    复杂热图 复杂散点图
    复杂热图02 甘特图 百分比柱状图 箱线图美化 弦图 mantel test图 瀑布图 曼哈顿图 KEGG富集图 哑铃图

    以上内容仅为群内部分内容,不代表全部。详细目录请看下链接!

    示例数据和代码获取

    往期文章

    1. 生信常用分析图形绘制01 -- 各种类型的热图!
    2.生信常用分析图形绘制02 -- 解锁火山图真谛!

    相关文章

      网友评论

        本文标题:生信常用分析图形绘制03 -- 富集分析圈图

        本文链接:https://www.haomeiwen.com/subject/rxafgrtx.html