【R语言作图】富集分析结果直方图

前言

最近小Q在做自然选择分析，分析完之后简单粗暴的对候选基因做了富集分析，并做了展示，比起气泡图，我模仿了另一种作图方式，显示效果更佳。所以想在此分享一下如何用R语言画富集分析示意图（非气泡图）。

画图思路

利用ggplot2+grid包进行画图，采用分面的思想作图。

画图代码

#导入所需R包
library(ggplot2)
library(grid)

#--------------------------------------
# 数据读入与预处理
#--------------------------------------
options(stringsAsFactors = FALSE)
dat<-read.table("exampleData2.txt",header = TRUE,sep = "\t")
#---------------
# 输入文件格式：Term Inputnumber Backgroundnumber PValue label
# Inputnumber: 你的基因列表和Term基因列表重叠的基因数目
# label：是你想特别强调的Trem，1表示特别强调，0表示不特别强调
                              Term Inputnumber Backgroundnumber      PValue label
                  Plasma omega-6           10              689 0.000465508     0
        Verbal declarative memory           9              340 0.000670600     0
              Phospholipid levels           8              170 0.001023562     1
                      Temperament          12              158 0.002094460     0
 Red blood cell fatty acid levels          11               43 0.003176734     1
         Oleic acid plasma levels           7               25 0.003248211     1
                 Aortic root size           7              183 0.003264350     0
#---------------
# 按照GO的pval排序，从小到大
dat<-dat[order(dat$PValue),]
# 给排序后的Term编号，用作X-axis
group<-c(1:nrow(dat))
dat<-cbind(group,dat)
# 将label设为因子变量，便于fill
dat$label<-factor(dat$label)
#--------------------------------------
# 作图
#--------------------------------------
# 通用参数
# 画图区域外的空白框设置
mg_l <- margin(0, 0, 0, 0, 'lines') # 方向: 上右下左

# 先分别用ggplot2画好每一个子图
# logP
mg_1 <- unit(c(0, 0, 0.3, 0), "lines")  # 方向: 上右下左
# 之所以求logP而不是-logP，是为了让bar是从顶部往下绘制，最终呈现图中向左绘制的效果
# 因为logP是负值，所以X tick的label需要自定义重新设置label
limit_1 <- c(min(log10(dat$PValue))-0.5, 0)
break_1 <- seq(from = 0, to = min(log10(dat$PValue)), by = -1)
label_1 <- seq(from = 0, to = abs(min(log10(dat$PValue))), by = 1)
# 之所以在这里第一个label设为空，是为了后面拼图时原点只有一个0，而不是两个0，更好看些
# 在geneCount的label里第一个label不会设置为空
label_1[1]<-""
# 这里想要强调的条目设置成红色，更醒目
# 这里颜色的fill = label
fillCol<-c(rgb(244,153,30,max=255),rgb(255,0,0,max=255))
fillCol_brk<-c(0,1)
# limits = c(0.5, (nrow(dat) + 0.5)), breaks = 1:nrow(dat)
# 为了后面统一各个子图上term位置的一致性
p_logP <- ggplot(dat) +
  geom_bar(aes(x=group, y=log10(PValue),fill=label), stat="identity", position="dodge",width = 0.5) +
  scale_y_continuous(limits = limit_1, breaks = break_1,labels = label_1,expand = expand_scale()) +
  scale_x_continuous(limits = c(0.5, (nrow(dat) + 0.5)), breaks = 1:nrow(dat),labels = NULL, expand = expand_scale(), position = 'top') +
  scale_fill_manual(values = fillCol,breaks=fillCol_brk)+
  theme(plot.margin = mg_1, axis.text = element_text(margin = mg_l),axis.ticks.y = element_blank()) +
  xlab(NULL) +
  ylab("-Log10(Pvalue)") +
  coord_flip() + 
  guides(fill = FALSE)

# GeneCount
mg_2 <- unit(c(0, 0, 0.3, 0), "lines")  # 方向: 上右下左
max_geneNum<-max(dat$Inputnumber)
limit_2 <- c(0,max_geneNum+1)
dig_temp<-nchar(as.character(max_geneNum))
by <- round(max_geneNum/(5 * 10^(dig_temp - 2))) * 10^(dig_temp - 2)
break_2 <- seq(from = 0, to = max_geneNum, by = by)
# 这里设置label是为了更美观，label没有小数点
label_2 <- seq(from = 0, to = max_geneNum, by = by)

p_geneCount <- ggplot(dat) +
  geom_bar(aes(x=group, y=Inputnumber), fill="#696969",stat="identity", position="dodge",width = 0.5) +
  geom_text(aes(x=group, y=Inputnumber,label=Inputnumber),hjust=-0.5)+
  scale_y_continuous(limits = limit_2,breaks = break_2,labels = label_2,expand = expand_scale()) +
  scale_x_continuous(limits = c(0.5, (nrow(dat) + 0.5)), breaks = 1:nrow(dat),labels = NULL, expand = expand_scale(), position = 'top') +
  theme(plot.margin = mg_2, axis.text = element_text(margin = mg_l),axis.ticks.y = element_blank()) +
  xlab(NULL) +
  ylab("No. of genes") +
  coord_flip() + 
  guides(fill = FALSE)

# Function term
text_x <- rep(1, nrow(dat))
text_y <- 1:nrow(dat)
text_lab<-dat$Term
mg_3 <- unit(c(0, 0.3, 2.3, 0), "lines") # 方向: 上右下左

p_GoName <- ggplot() + 
  geom_text(aes(x = text_x, y= text_y, label = text_lab), hjust=1,size = 4) +
  scale_x_continuous(limits = c(0,1), breaks = NULL, expand = expand_scale()) +
  scale_y_continuous(limits = c(0.5, (nrow(dat) + 0.5)), breaks = NULL, expand = expand_scale()) +
  labs(x = NULL, y = NULL) + 
  theme(plot.margin = mg_3,
        axis.text = element_text(margin = mg_l),
        panel.grid.major =element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) 
# 分面画图
# 写一个简易函数调用viewport
vplayout <- function(x, y){
  viewport(layout.pos.row = x, layout.pos.col = y)
}

grid.newpage()  ##新建页面
pushViewport(viewport(layout = grid.layout(1, 3)))
# goTerm logP GeneNum
print(p_GoName, vp = vplayout(1, 1))
print(p_logP, vp = vplayout(1, 2))
print(p_geneCount, vp = vplayout(1, 3))

最终效果图：

转载请标明出处和作者 ^+^

撰文 & 编辑：VickieQ
校对：HCLO4 & 花毛