火山图（Volcano plot）是一种比较“远古”的一种散点图，广泛应用于转录组、蛋白组等组间差异分析结果的展示。图表的横轴一般展示差异倍数的变化，而纵轴表示差异分析结果的可靠性。

我这里对传统的火山图做了进一步的优化，比如通过调整点的大小突出展示感兴趣基因对应的点，并将感兴趣基因的名称以标签的方式展示出来，如下。

那么，如何绘制这般个性又好看的火山图呢？我这里主要用到的还是最常见的ggplot2包。

#设置工作目录；

setwd( "C:/Users/MHY/Desktop/火山图的绘制")

dir

#读入数据；

data<- read.csv( "protein_dia_simple.csv",header = T)

head(data)

#载入相关的R包；

library(dplyr)

library(ggplot2)

library(ggrepel)

#转成tibble便于后续使用，去掉不需要的列；

dt<- as_tibble(data[c(-2,-3,-5)])

head(dt)

#对Q值取对数；

dt$log10FDR <- -log10(dt$qvalue)

#生成显著上下调数据的分组标签；

dt$group <- case_when(dt$log2fc > 1 & dt$qvalue < 0.05 ~ "Up",

dt$log2fc < -1 & dt$qvalue < 0.05 ~ "Down",

abs(dt$log2fc) <= 1 ~ "None",

dt$qvalue >= 0.05 ~ "None")

head(dt)

#获取表达差异最显著的10个基因；

top10sig<- filter(dt,group!= "None") %>% distinct(Symbol,.keep_all = T) %>% top_n(10,abs(log2fc))

top10sig

#将差异表达Top10的基因表格拆分成up和down两部分；

up<- filter(top10sig,group== "Up")

up

down<- filter(top10sig,group== "Down")

down

#新增一列，将Top10的差异基因标记为2，其他的标记为1；

dt$size <- case_when(!(dt$id %in% top10sig$id)~ 1,

dt$id %in% top10sig$id ~ 2)

head(dt)

#提取非Top10的基因表格；

df<- filter(dt,size==1)

head(df)

#指定绘图顺序,将group列转成因子型；

df$group <- factor(df$group,

levels= c( "Up", "Down", "None"),

ordered= T)

#开始绘图，建立映射；

p0<-ggplot(data=df,aes(log2fc,log10FDR,color=group))

#添加散点；

p1<- p0+geom_point(size=1.6)

p1

#自定义半透明颜色（红绿）；

mycolor<- c( "#FF9999", "#99CC00", "gray80")

p21<- p1 + scale_colour_manual(name= "",values=alpha(mycolor,0.9))

p21

#其他配色方案：

mycolor<- c( "#FF99CC", "#99CC00", "gray80")

p22<- p1 + scale_colour_manual(name= "",values=alpha(mycolor,0.7))

p22

#继续添加Top10基因对应的点；

p2<- p22+geom_point(data=up,aes(log2fc,log10FDR),

color= "#FF9999",size=3,alpha=0.9)+

geom_point(data=down,aes(log2fc,log10FDR),

color= "#7cae00",size=3,alpha=0.9)

p2

调整横轴和纵轴绘图区域的范围，重点突出Top10基因；设置y轴范围（上下两端的空白区域设为1），修改刻度标签。

#expansion函数设置坐标轴范围两端空白区域的大小；mult为“倍数”模式，add为“加性”模式；

p3<-p2+labs(y= "-log10FDR")+

scale_y_continuous(expand=expansion(add = c(2, 0)),

limits= c(0, 40),

breaks= c(0,10,20,30,40),

label= c( "0", "10", "20", "30", "40"))+

scale_x_continuous(limits = c(-4, 4),

breaks= c(-4,-2,0,2,4),

label= c( "-4", "-2", "0", "2", "4"))

p3

接下来为top10基因对应的散点添加指引线，我们主要用到geom_text_repel和geom_label_repel这两个函数，下面是比较重要的一些参数：

nudge_x/y：数据点与相应数据标签的距离，例如1表示标签在点右/上的1个单位处，而-2.2表示标签在点左/下2.2个单位处；

direction：标签分布方向，x表水平分布，y 表示垂直分布，both 表示随机分布；

segment.size：指定线段的粗细；

point.padding：表示点周围的空余区域，决定连接线端点到到数据点中心的距离，单位为line。

#添加箭头；

set.seed(007)

p4<- p3+geom_text_repel(data=top10sig,aes(log2fc,log10FDR,label=Symbol),

force=80,color= "grey20",size=3,

point.padding = 0.5,hjust = 0.5,

arrow= arrow(length = unit(0.01, "npc"),

type= "open", ends = "last"),

segment.color= "grey20",

segment.size=0.2,

segment.alpha=0.8,

nudge_x=0,

nudge_y=1)

p4

#自定义图表主题，对图表做精细调整；

top.mar=0.2

right.mar=0.2

bottom.mar=0.2

left.mar=0.2

#隐藏纵轴，并对字体样式、坐标轴的粗细、颜色、刻度长度进行限定；

mytheme<-theme_classic+

theme(text=element_text(family = "sans",colour = "gray30",size = 12),

axis.line = element_line(size = 0.6,colour = "gray30"),

axis.ticks = element_line(size = 0.6,colour = "gray30"),

axis.ticks.length = unit(1.5,units = "mm"),

plot.margin=unit(x=c(top.mar,right.mar,bottom.mar,left.mar),

units= "inches"))

#应用自定义主题；

p4+mytheme

#添加辅助线；

p5<- p3+geom_hline(yintercept = c(-log10(0.05)),

size= 0.7,

color= "orange",

lty= "dashed")+

geom_vline(xintercept = c(-1,1),

size= 0.7,

color= "orange",

lty= "dashed")

p5

#添加其他样式的标签；

#为了方便自定义左右区域的标签，这里使用up、down两个独立的子表格；

p6<- p5+geom_label_repel(

data= up,aes(log2fc,log10FDR,label=Symbol),

nudge_x= 1,

nudge_y= 5,

color= "white",

alpha= 0.9,

point.padding = 0.5,

size= 3,

fill= "#96C93D",

segment.size = 0.5,

segment.color = "grey50",

direction= "y",

hjust= 0.5) +

geom_label_repel(

data= down,aes(log2fc,log10FDR,label=Symbol),

nudge_x= -1,

nudge_y= 3,

color= "white",

alpha= 0.9,

point.padding = 0.5,

size= 3,

fill= "#9881F5",

segment.size = 0.5,

segment.color = "grey50",

direction= "y",

hjust= 0.5)

#应用自定义主题；

p7<- p6+mytheme

p7

好啦，以上就是今天分享的全部内容，你学会了吗？

转自：基迪奥生物