前言
最近小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 & 花毛
