微信公眾號(hào)：生物信息學(xué)習(xí)
火山圖（volcano plot）常用于顯著差異基因表達(dá)的展示，包含顯著和差異兩個(gè)重要信息艾船。
那么如何看懂一張火山圖所包含的信息呢？首先需要知道屿岂，火山圖的橫坐標(biāo)通常用log2（fold change）表示鲸匿，差異越大的基因分布在兩端，縱坐標(biāo)用-log10（pvalue）表示带欢，T檢驗(yàn)顯著性P值的負(fù)對(duì)數(shù)烤惊。由于P值越小表示越顯著吁朦，所以我們進(jìn)行-log10（P value）轉(zhuǎn)化后，轉(zhuǎn)化值越大表示差異越顯著扰才。通常差異倍數(shù)越大的基因T檢驗(yàn)越顯著，所以左上角和右上角的值往往是我們關(guān)注的纺讲。
隨機(jī)生成log2FoldChange和padj的數(shù)值
log2FoldChange<-sample(rnorm(10000,mean=0,sd=4),1000,replace=T)
padj<-sample(runif(10000,0,0.08),1000,replace=T)
#log2FoldChange和padj合并到一個(gè)數(shù)據(jù)框
mydata <- data.frame(log2FoldChange,padj)
##對(duì)滿足不同條件的數(shù)據(jù)給不同的標(biāo)記，放入Condition列逢渔，顏色放入color列
mydata$Condition=ifelse(mydata$log2FoldChange>=1 & mydata$padj<=0.05,"up",ifelse(mydata$log2FoldChange<=-1 & mydata$padj<=0.05,"down","normal"))
mydata$color=ifelse(mydata$log2FoldChange>=1 & mydata$padj<=0.05,"red",ifelse(mydata$log2FoldChange<=-1 & mydata$padj<=0.05,"green","gray"))
接下來(lái)就開(kāi)始畫(huà)圖乡括。
1.直接用plot畫(huà)圖
plot(mydata$log2FoldChange,-log10(mydata$padj),col=mydata$color,pch=20,cex=.3,xlab="log2 fold change",ylab="-log10 p-value",font.lab=2,font.axis=2,cex.axis=.9)
##用abline給幾條虛線
abline(v=1,lty=3)
abline(v=-1,lty=3)
abline(h=-log10(0.05),lty=3)
結(jié)果圖如下所示：
image
2.用ggplot2畫(huà)圖
library(ggplot2)
p <-ggplot(data=mydata, aes(x=log2FoldChange, y=-log10(padj), colour=Condition)) + geom_point(alpha=0.8, size=1) + xlab("log2 fold change") + ylab("-log10 padj")+geom_hline(yintercept=-log10(0.05),linetype=4)+geom_vline(xintercept=c(-1,1),linetype=4)+scale_color_manual(values=c('up'='red','down'='green','normal'='gray'))
p
結(jié)果如下圖所示：
image
關(guān)注微信公眾號(hào)：生物信息學(xué)習(xí)粟判，并私信：volcano亿昏，即可獲得包含了如何保存圖片的代碼档礁。
