生信分析学习笔记：（1）火山图

由于一直没认真学习，经常性三天打鱼两天晒网的，导致自己落后于他人太多【大哭】，所以要立个flag，争取早日追上大家。
今天第一天，先从最简单的 火山图 开始学起。

现在看一下效果图

背景知识介绍：

1、每个点代表一个检测到的基因。
2、横轴和纵轴用于固定点在空间的位置。
一般横轴是Log2(foldchange)，点越偏离中心，表示差异倍数越大。
纵轴是-Log 10 (adjusted P-value)，点越靠图的顶部表示差异越显著。
3、点的大小和颜色也可以表示更多的属性，如下图中点的颜色标记其对应的基因是上调, 下调还是无差异。
4、大小也可用于展示基因表达的平均丰度，一般我们关注表达水平较高且差异较大的基因用于后续的分析和验证

1、什么是fold change?
翻译成中文是差异倍数，简单来说就是基因在一组样品中的表达值的均值除以其在另一组样品中的表达值的均值。所以火山图只适合展示两组样品之间的比较。

2、 为什么要做Log 2转换？
两个数相除获得的结果 （fold change）要么大于1，要么小于1，要么等于1。这是一句正确的废话吧？那么对应于基因差异呢？简单说，大于1表示上调（可以描述为上调多少倍），小于1表示下调（可以描述为下调为原来的多少分之多少）。大于1可以到多大呢？多大都有可能。小于1可以到多小呢？最小到0。用原始的fold change描述上调方便，描述下调不方便。绘制到图中时，上调占的空间多，下调占的空间少，展示起来不方便。所以一般会做Log2转换。默认我们都会用两倍差异 （fold change == 2 | 0.5）做为一个筛选标准。Log2转换的优势就体现出来了，上调的基因转换后Log2 (fold change)都大于等于1，下调的基因转换后Log2 (fold change)都小于等于-1。无论是展示还是描述是不是都更方便了。

3、P-value都比较熟悉，统计检验获得的是否统计差异显著的一个衡量值，约定成俗的P-value<0.05为统计检验显著的常规标准。

4、 什么是adjusted P-value?
这里面就涉及到一个统计学问题了。做差异基因检测时，要对成千上万的基因分别做差异统计检验。统计学家认为做这么多次的检验，本身就会引入假阳性结果，需要做一个多重假设检验校正。
这个校正怎么做呢？最简单粗暴的方法是每一次统计检验获得的P-value都乘以总的统计检验的次数获得adjusted P-value (这就是Bonferroni correction）。
但这样操作太严苛了，很容易降低统计检出力，找不到有差异的基因。后续又有统计学家提出相对不这么严苛的计算方法，如holm, hochberg, hommel, BH, BY, fdr等。BH是我们比较常用的一个校正方法，获得的值是假阳性率 FDR （false discovery rate）。
FDR筛选时就可以不用遵循0.05这个标准了。我们可以设置FDR<0.05表示我们容许数据中存在至多5%假阳性率；FDR<0.1表示我们对假阳性率的容忍度至多是10%。当然如果说我们设置FDR<0.5，即数据中最多可能有一半是假阳性就说不过去了。

5、 同样为什么做 -Log 10转换呢？
因为FDR值是0-1之间，数值越小越是统计显著，也越是我们关注的。-Log 10 > (adjusted P-value)转换后正好是反了多来，数值越大越显著，而且以10为底很容易换算回去。

数据结构

代码展示

###需要的程序包-------------------------------------
library(openxlsx) ##读取数据
library(ggplot2) ##绘图
library(ggrepel)


setwd("E:\\volcano") #设置工作路径
data<-read.xlsx("volcanoIN.xlsx")
data$threshold = factor(ifelse(data$Pvalue < 0.05 & abs(data$logFC) >= 1, 
                        ifelse(data$logFC>= 1 ,'Up','Down'),'NoSignifi'),
                        levels=c('Up','Down','NoSignifi'))

p2 <- ggplot(data, aes(x=logFC, y=-log10(Pvalue), color=threshold))+
  geom_point()+
  scale_color_manual(values=c("#DC143C","#00008B","#808080"))+#确定点的颜色
  geom_text_repel(
    data = data[data$Pvalue<0.05 & abs(data$logFC)>