利用ggplot2绘制功能富集气泡图

本期学习笔记的内容是在R语言中绘制功能富集泡泡图，最终能实现下面这种结果。

学习绘图之前，首先学一下怎么看这张图：

• 横轴表示差异基因占比，数值越大表示富集程度越大
• 纵轴表示富集的通路，越靠上的通路越显著
• 气泡的颜色表示pvalue的大小（越红表示pvalue的负对数越大，pvalue越小）
• 气泡大小表示差异基因个数，越大则数量越多

功能富集分析泡泡图用来展示一组基因参与哪些功能调节通路，有利于理解基因的功能和潜在意义，这次跟着生信宝典的R语言学习教程，记录一下重点内容笔记。

准备阶段

富集矩阵

• Description：通路描述
• GeneRatio：该通路的差异基因占总差异基因比例
• qvalue：富集的显著程度
• Count：该通路差异基因的个数
• Type：区分样本

整理数据

对富集矩阵进行处理，读取为表格形式，设置分隔符和行名，并输出整理后的数据。

e_data <- read.table(text=data,header = T,row.names = NULL,
                     sep = "\t",quote="")
head(e_data)

加载R包

library(plyr)
library(stringr)
library(grid)
library(ggplot2)

需要用到如上所述R包，其中ggplot2用于绘图。

单样本绘制

原始数据包含两个样本，这里先取出数据中的一个样本。

e_data_1=droplevels(e_data[e_data$Type=="Baodian_up",])

构造一个函数，使分数转化为小数。

mixedToFloat <- function(x){
  x <- sapply(x, as.character)
  is.integer  <- grepl("^-?\\d+$", x)
  is.fraction <- grepl("^-?\\d+\\/\\d+$", x)
  is.float <- grepl("^-?\\d+\\.\\d+$", x)
  is.mixed    <- grepl("^-?\\d+ \\d+\\/\\d+$", x)
  stopifnot(all(is.integer | is.fraction | is.float | is.mixed))
  
  numbers <- strsplit(x, "[ /]")
  
  ifelse(is.integer,  as.numeric(sapply(numbers, `[`, 1)),
         ifelse(is.float,    as.numeric(sapply(numbers, `[`, 1)),
                ifelse(is.fraction, as.numeric(sapply(numbers, `[`, 1)) /
                         as.numeric(sapply(numbers, `[`, 2)),
                       as.numeric(sapply(numbers, `[`, 1)) +
                         as.numeric(sapply(numbers, `[`, 2)) /
                         as.numeric(sapply(numbers, `[`, 3)))))
  
}

进行测试一下，看看刚才构造的函数有什么用处，可以发现该函数能将输入的分数转化为小数输出。

> mixedToFloat(c('2 3/4', '2/3', '11 1/4'))
[1]  2.7500000  0.6666667 11.2500000

转化数据

利用刚才构造的函数对数据进行转化，将基因的比例和计数转化为小数格式，并设置矩阵的列名；qvalue转换为负对数，并添加为新的一列。

> e_data_1$GeneRatio=mixedToFloat(e_data_1$GeneRatio)
> e_data_1$Count=mixedToFloat(e_data_1$Count)
> log_name <- paste0("negLog10_","qvalue")
> col_name_e_1 <- colnames(e_data_1)
> col_name_e_1 <- c(col_name_e_1,log_name)
> e_data_1$log_name <- log10(e_data_1$qvalue) * (-1)
> colnames(e_data_1) <- col_name_e_1

排序统计

获取数据中通路的出现次数并进行排序，获得Description和ID信息。

> e_data_1_freq <- as.data.frame(table(e_data_1$Description))
> colnames(e_data_1_freq) <- c("Description","ID")
> head(e_data_1_freq)
                                       Description ID
1                            ERK1 and ERK2 cascade  1
2   negative regulation of Notch signaling pathway  1
3                         neuron apoptotic process  1
4             positive regulation of cell adhesion  1
5 positive regulation of cytoskeleton organization  1
6                        regulation of cell growth  1

将两个矩阵合并到一起，以Description列为基准，然后对合并后的矩阵排序（按照ID、GeneRatio、负对数值为序）

e_data_2 <- merge(e_data_1,e_data_1_freq,by="Description")

e_data_3 <- e_data_2[order(e_data_2$ID,
                           e_data_2$GeneRatio,
                           e_data_2$negLog10_qvalue),]

unique删除重复的Description行，然后转换成因子，并进行排序。

t_order <- unique(e_data_3$Description)
e_data_1$Description <- factor(e_data_1$Description,
                               levels = t_order,ordered = T)

绘图阶段

设置颜色

color_1 <- c("green","red")

初始化坐标系

设置x轴和y轴，增加标签x="GeneRatio",y="GO description"

p <- ggplot(e_data_1,aes(x=GeneRatio,y=Description)) +
  labs(x="GeneRatio",y="GO description") + labs(title="")

设置点状图的大小和颜色信息，scale_color_gradient用于生成一组渐变色。

p <- p + geom_point(aes(size=Count,color=negLog10_qvalue)) +
  scale_color_gradient(low = color_1[1],high=color_1[2],name="negLog_qvalue")

对纵轴的标签进行处理，str_wrap用于设置字符串单行长度不超过60.

p <- p + scale_y_discrete(labels=function(x) str_wrap(x,width = 60))

设置绘图区域参数

初始化参数如下：

top="top"
bottom="bottom"
left="left"
right="right"
none="none"
legend_pos_par=right

使用判断语句来自动估算图形的长宽数据，设置平面大小。

uwid=0
vhig=12

if (uwid == 0 || vhig == 0) {
  x_len = length(unique(e_data_1$Description))
  if(x_len<10){
    vhig = 10
  } else if(x_len<20) {
    vhig = 10 + (x_len-10)/3
  } else if(x_len<100) {
    vhig = 13 + (x_len-20)/5
  } else {
    vhig = 40
  }
  uwid = vhig 
  if(legend_pos_par %in% c("left", "right")){
    uwid = 1.5 * uwid
  }
}

调整气泡图

使用如下代码调整图像，步骤：设置图例的位置——设置边界和背景——设置坐标轴参数——输出plot。

p <- p + theme(legend.position = legend_pos_par)
p <- p + theme(panel.grid=element_blank(),
               panel.border = element_blank(),
               legend.background = element_blank(),
               axis.line.x=element_line(size=0.4,colour = "black",linetype = "solid"),
               axis.line.y=element_line(size=0.4,colour = "black",linetype = "solid"),
               axis.ticks=element_line(size=0.4))
p

本文数据和代码笔记获取链接：https://z.jewin.love/20221104.R

参考资料：

http://www.ehbio.com/Bioinfo_R_course