volcano plot 火山图精修

首先我们必须明确今天的任务就是，把左边我们昨天画的草图精修成类似右边那样式儿的，从谷歌图片上搜的符合paper格式要求的图：

草图 VS 精修
昨天简书那张图是我从百度上随便搜的，今天想了一下，还是从谷歌上搜一张文章里确实真实存在的图吧，比较有说服力

首先我们对比一下看需要修改哪里？

• 坐标轴需要修改，一般X轴我们只要[-4,4]之间的点；
• 在最底端有很多基点，也就是-log10(p-value)=0的点不要；
• 图上的颜色太难看了，我们需要把显著的点找出来，给它一个蓝色，不显著的用灰色；
• p-value=0.05，（即Y轴：-log10(0.05)=1.30103）这根辅助线添加上；
• X轴在1和-1处的辅助线加上
  下面我一个一个修改
  首先我们先把昨天的代码跑一遍，复习下昨天的代码：

rm(list=ls())  #好习惯要养成，先清空下环境变量
volcano_plot = read.table(file="gene_exp.diff文件路径",header = TRUE)
log2_foldchange = log2(volcano_plot$value_1 / volcano_plot$value_2) #X轴

log2_foldchange[volcano_plot$value_2 == 0 ] = 0  #筛选X轴
log2_foldchange[volcano_plot$value_1 == 0 ] = 0 #筛选X轴

log10_p_value = log10(volcano_plot$p_value) * -1 #准备Y轴
#X/Y轴都准备好了，画草图
plot(x=log2_foldchange,y=log10_p_value)

昨天那么啰嗦，其实就是讲了这么几行代码

下面开始正式精修图

1. 修正X/Y轴，去掉-log10(p-value)=0的点

用xlim=c(-4,4)参数限定X轴只取[-4,4]；
用ylim=c(0,4)参数限定Y轴只取[0,4]
我们不想要最底下那一横排的点（即-log10(p-value)=0的点）。所以，下面我们对-log10(p-value)进行筛选，把等于0的点过虑掉：

log10_p_value_qc = log10_p_value[log10_p_value >= 0.001]
log2_foldchange_qc = log2_foldchange[log10_p_value >= 0.001]

用过滤后的X/Y轴再次画图，X轴取[-4,4]；Y轴取[0，4]

plot(x=log2_foldchange_qc, y=log10_p_value_qc, xlim=c(-4,4), ylim=c(0,4))

从图中我们可以看出来，最底下那一横排的点确实去掉了。

去掉-log10(p-value)=0的点

2.改变颜色

思路：所有点先改成灰色 -> 找出显著性的点 -> 将显著性的点变成蓝色

所有点先改成灰色

plot(x=log2_foldchange_qc, y=log10_p_value_qc,
     xlim=c(-4,4),ylim=c(0,4),
     col="#BCBABE",pch=16
     )

灰色的颜色编号是"#BCBABE"，所以给一个参数col="#BCBABE"所有点的颜色就变成灰色了
pch=16的意思是选择了一种点的样式，不同样式的点对应着不同的编号。

所有点先改成灰色

找出显著性的点

我们先认为同时满足这些条件的就是显著的点：

• p_value <= 0.05 ；即Y轴log10_p_value_qc >= 1.30103
• fold change 大于2，或小于二分之一 ，即X轴的绝对值大于等于1abs(log2_foldchange_qc) >= 1。
  为了方便后续我们操作，现在我们把同时满足这两个条件的基因筛选出来并写进一个向量sign_point

sign_point = (abs(log2_foldchange_qc) >= 1) & (log10_p_value_qc >= 1.30103)

将所有点的颜色编号#BCBABE(灰色)放入到向量col_point里：

> length(log2_foldchange_qc)
[1] 12769

注：length(log2_foldchange_qc)在查看，过滤之后还有几个点，即剩下多少个基因。

> col_point = rep("#BCBABE", length(log2_foldchange_qc))
> col_point
   [1] "#BCBABE" "#BCBABE" "#BCBABE" "#BCBABE" "#BCBABE" "#BCBABE" "#BCBABE" "#BCBABE"
   [9] "#BCBABE" "#BCBABE" "#BCBABE" "#BCBABE" "#BCBABE" "#BCBABE" "#BCBABE" "#BCBABE"
  [17] "#BCBABE" "#BCBABE" "#BCBABE" "#BCBABE" "#BCBABE" "#BCBABE" "#BCBABE" "#BCBABE"

把筛选出来的显著性的点sign_point变成蓝色("#1B2CC1")

也就是把col_point向量里面的sign_point点所对应的 "#BCBABE"(灰色编号)变成蓝色("#1B2CC1")

col_point[sign_point] = "#1B2CC1"

此时，存放所有点的向量col_point就变成了，显著的点是蓝色"#1B2CC1"，不显著的点是灰色了，如图：

查看向量col_point
我们用做好的颜色向量重新作图，为显著的点加上颜色

plot(x=log2_foldchange_qc, y=log10_p_value_qc,
     xlim=c(-4,4),ylim=c(0,4),
     col=col_point,pch=16
     )

显著的点加上了蓝色

4.加辅助线——参数abline()

abline(h=-1*log10(0.05),lwd=3,lty=3,col="#4C5B61")
abline(v=log2(2) ,lwd=3,lty=3,col="#4C5B61")
abline(v=log2(1/2) ,lwd=3,lty=3,col="#4C5B61")

终稿图

解释各参数：

• h是加一条水平线（horizontal line）；
• v是加一条竖直线（vertical curve）；
• lwd线条加粗，且选3号样式；
• lty线条做成虚线；
• col给线一个色号"#4C5B61"

大家是否好奇？又不是女生买口红，哪个号对应着哪个颜色都知道。R里我是怎么知道蓝色对应的色号是"#1B2CC1"；灰色对应的色号是"#BCBABE"的
这个我们明天讲~

---

课程分享
生信技能树全球公益巡讲
（https://mp.weixin.qq.com/s/E9ykuIbc-2Ja9HOY0bn_6g）
B站公益74小时生信工程师教学视频合辑
（https://mp.weixin.qq.com/s/IyFK7l_WBAiUgqQi8O7Hxw）
招学徒：
（https://mp.weixin.qq.com/s/KgbilzXnFjbKKunuw7NVfw）