最近很多刚了解生信的同学问喵学姐：看了一些文献，文献里的各种图怎么看呀，完全看不懂。今天喵学姐就来给大家讲一讲我们平时做的最基础的差异分析——火山图
火山图（Volcano plot）是散点图的一种，它将统计测试中的统计显著性量度（如p-value）和变化幅度相结合，从而能够帮助快速直观地识别那些变化幅度较大且具有统计学意义的数据点（代谢物等）。是一种单变量统计分析方法，常应用于研究基因组、转录组、代谢组、蛋白质组等数据分析。
⌈火山图长啥样？⌋
先给大家看看最最最常见的一种，基础款：
有的人可能觉得，这不就是高级一点的散点图+三根虚线+花里胡哨的颜色吗？
NO NO NO，火山图可不光只有上面那一种，
▲别人家的火山图，是不是很美
▲别人的火山图里竟然还有反比例函数阈值线
▲这一排火山图放在一起是不是很有气势
⌈火山图有啥用？⌋
分析数据用散点图表示看不大明白？
想要表现数据间有统计学意义的明显差异？
我们给散点图换个样子，分区划线做成⌈火山图⌋优势这不就来了？
火山图是生信文章中最常见的表达图，它最直观的作用就是展示组间差异的物质。
比如找到肿瘤细胞跟正常细胞的差异基因，或是研究疾病组跟正常组之间的差异表达基因情况。
绘制火山图通常需要两个关键指标，
①差异倍数Fold Change值
②假设检验P值
⌈火山图怎么看？⌋
用之前讲过的一篇非肿瘤思路文献来举例，：4分＋非肿瘤纯生信，GEO数据集+铁死亡+cytoscape调控网络+miRNA+转录因子
这张图展现了精神分裂组与正常组之间的差异表达基因，
其中每个点代表一个基因。
蓝色的点代表显著下调的基因，红色的点代表显著上调的基因，而灰色的点就表示无显著差异的基因。
那我们怎么去判断它是显著上调还是显著下调呢？
这就涉及到两个坐标轴，
② 横坐标Log2 (Fold Change)
②纵坐标-log10 (P value)
横坐标为log2FC值即差异倍数值。表示两个分组之间的差异倍数，其绝对值越大说明某基因在两组之间的表达差异也越大。该值为正时，表示差异上调；该值为负时，表示差异下调。
纵坐标为-log10（P-value）值，表示某个基因在比较分组之间的表达差异是否足够显著，一般认为p-value<0.05为显著。
根据作者设置的阈值：log2FC的绝对值为0.1，p值<0.05，因此灰色部分的基因不满足条件。
因此，横着看基因与对照组的差异，竖着看有无统计学意义，在欣赏火山图时，我们只需要关注左上和右上两区，然后再研究这两块的基因就好啦。
⌈火山图怎么做？⌋
之前的推文已经为大家详细讲解了整篇文献的思路框架，也教大家学习了GEO数据库以及数据集的下载。
现在呢，我们就要开始复现文章的结果图了，先来做GEO数据差异基因的分析，今天会使用两个方法教大家绘制火山图。
①使用GEO2R绘制火山图
Step 1. 找到数据集
文献中作者用到的数据集是GSE27383，我们直接在GEO数据库搜索这个数据集。
它包含了精神分裂症患者的信息和健康对照者的患者信息。
Step 2. 使用GEO2R进行分析
直接点击GEO2R，
因为需要做两组之间的差异基因，所以要对组别进行定义。
Step 3. 对组别进行定义
首先在Define groups定义精神分裂症组别，命名为SCH，点击enter键，将正常对照组命名为HC。
然后，把相应的样本归属到不同组别中。直接选中相应的精神分裂症的样本，把这43个样本归属到SCH组别中。
接着同样的操作，把正常对造组的29个样本归属到HC。
Step 4. 分析生成火山图
然后点击Analyze分析，就可以生成相应的火山图，底下是所有的差异基因信息。包括基因探针名、校正的P值、P值、T检验、logFC值、基因名，直接下载就可以得到所有的差异基因的信息了。
直接将它复制到excel表格中打开。可以看到一共有45000多个差异基因，但它们并不都参与精神分裂症的发生。
Step 5. 设置阈值，筛选差异基因
第一种方法：手动设置阈值
我们需要根据文献中作者设置的阈值：p<0.05 and ︱logFC︱>0.1，从这两个层面去限定阈值作为精神分裂症与健康对造组之间的显著差异基因。
对刚刚下载的数据进行筛选，首先使用ABS函数对logFC进行绝对值计算。
然后根据作者设置的阈值筛选我们想要的显著性的差异基因。
最终得到4000多个参与精神分裂症的显著差异基因。
第二种方法：直接在GEO2R网站上设置阈值
点击Options这一栏，并在界面右侧进行阈值设置。
0.05为cut-off值， log2FC值为0.1。
点击分析后出来的火山图可以直接下载使用。
Step 6. 差异基因数据下载
底部是经过设置阈值后呈现的所有基因情况，也可以直接下载。
②使用在线分析平台绘制火山图
Step 1. 打开桑格助手分析工具http://vip.sangerbox.com/home.html
对于没有 GEO2R 这个功能的数据集，我们可以直接在桑格助手在线平台去做差异基因分析，只要上传表达谱，然后输入组别信息，就可以得到相应的差异基因情况了。
点击⌈工具列表-limma差异分析工具⌋进入使用界面。
Step 2. 上传表达谱信息以及组别信息
上传处理好的表达谱文件。表达谱加载完后，左边这一列默认为表达谱里的所有样本。
这时我们需要用到临床数据来分别定义精神分裂症组别和健康对照组别。
打开从GEO数据库下载的临床数据集进行组别定义。在右边新建一列，先把健康对照组筛选出来，组别命名为HC，然后对空白处的精神分裂组命名为SCH。
再把组别信息复制到桑格助手里面，这样我们所需要的数据就准备好了。
Step 3. 设置阈值
在右边的实验设计选择实验组vs对照组，差异方法选择limma检验。然后开始分析，进行命名并确认。
此处需要点击选择比较组来显示火山图，PS：因为火山图耗的耗费的内存比较大，所以默认不显示。
左侧为所有的差异基因信息，右边为火山图。也可以选择阈值参数，作者限定的是logFC值，PS：这个工具有一个缺点就是不能定义logFC值，只能定差异倍数，我们一般分析都认为若差异倍数为两倍则显著上调或显著下调。
但是在此数据集中，我们设置差异倍数为2即log2FC为1的时候，基本上没有显著性差异的基因，所以说在文章中，作者就放宽了条件限制，设置为log2FC= 0.1 and p<0.05。
Step 4. 基因差异信息下载
页面底部为基因差异信息，可以展现出就是在不同阈值设置的情况下相应的差异倍数上调或下调的基因数量。
在右边，有未设置阈值/设置阈值后的所有的显著性差异基因数据情况，都可以直接下载使用。
END
只要原始数据准备好，阈值设置好，大家也可以使用其他的在线分析平台去绘制火山图，喵学姐在这里就不赘述了。
火山图绘制的底层逻辑已经告诉大家了，后面大家也可以更近一层，学习用代码把火山图做的更好看。
后面喵学姐再慢慢给大家讲解其他的结果图，新来的同学们可以悄悄点个关注~
最后附上双曲线火山图代码
library(ggplot2)
volcano_plot <- function(df, pvalue_threshold = 0.05, foldchange_threshold = 1) {
xmax <- max(abs(na.omit(df$log2foldchange))) + 0.2
xmin <- min(abs(na.omit(df$log2foldchange)), 0.0001)
x <- seq(xmin, xmax, by = 0.0001)
y <- 1/x + (-log10(pvalue_threshold))
curve_xy <- rbind(data.frame(xpos = x + foldchange_threshold, ypos = y),
data.frame(xpos = -(x + foldchange_threshold), ypos = y))
df$curve_y <- ifelse(df$log2foldchange > 0,
1/(df$log2foldchange - foldchange_threshold) + (-log10(pvalue_threshold)),
1/(-df$log2foldchange - foldchange_threshold) + (-log10(pvalue_threshold)))
df$curve_group <- ifelse(-log10(df$pvalue) > df$curve_y & df$log2foldchange > foldchange_threshold, 'up',
ifelse(-log10(df$pvalue) > df$curve_y & df$log2foldchange < -foldchange_threshold, 'down', 'nosignif'))
df$pvalue <- -log10(df$pvalue)
p <- ggplot(df, aes(x = log2foldchange, y = pvalue, color = curve_group)) +
geom_point(size = 1) +
geom_line(data = curve_xy, aes(x = xpos, y = ypos), lty = 3, col = "black", lwd = 0.6) +
scale_color_manual(values = c('up'='red', 'down'='blue', 'nosignif'='gray')) +
xlim(-xmax, xmax) +
ylim(0, 30) +
labs(x = "log2(FoldChange)", y = "-log10(P-value)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
legend.spacing.x = unit(0.05, 'cm'), plot.title = element_text(hjust = 0.5),
legend.text = element_text(size = 8)) +
guides(color = guide_legend(override.aes = list(size = 2), title = NULL))
return(p)
}
data <- read.table('desktop/sample_dge.txt',header=T,stringsAsFactors=F,sep='\t')
head(data)
gene pvalue log2foldchange
1 ENSG00000000003.15 TSPAN6 6.954955e-04 1.0305811
2 ENSG00000000005.6 TNMD 1.103522e-01 -2.1289526
3 ENSG00000000419.12 DPM1 7.168680e-02 0.5515042
4 ENSG00000000457.14 SCYL3 5.743836e-01 0.1453620
5 ENSG00000000460.17 C1orf112 1.173320e-06 2.1643651
6 ENSG00000000938.13 FGR 1.388476e-13 -4.0345022
p <- volcano_plot(data)
p
本文使用 文章同步助手 同步
