【画图】如何画火山图？
【画图】如何画火山图？
Chris生命科学小站
发布于 2023-02-28 10:36:09
发布于 2023-02-28 10:36:09
代码可运行
运行总次数：0
代码可运行
火山图，可以展现基因总体的表达情况，并且可以凸显某个基因的重要性，是展现效率极高的图形。
如果你不喜欢敲代码可以按照下面的链接在线完成 2分钟！使用小站工具，就能用鼠标点出火山图和GSEA图~回复：我要测试。领取测试数据玩起来吧~
下面的数据与之前预测ACE2基因功能是一套数据，这里举例说明
画图
1. 加载所用到的包
library(limma)
library(annoE)
library(ggplot2)
library(ggrepel)
library(ggsci)
library(stringr)
library(data.table)备注：上面annoE是站长自己写的注释基因的包。
2. 计算差异基因
Genename<-"ACE2"
gene<-geneinfo(Genename)$ID
GTEx_ph<-fread(input='Lung.csv', stringsAsFactors = T,header = T,check.names = F)
lungTMP<-data.frame(GTEx_ph,row.names = 1,check.names = F)
#lungTMP<-read.csv("Lung.csv",header = T,check.names = F,row.names = 1)
ACE2geneexp<-data.frame(t(lungTMP[str_sub(rownames(lungTMP),end=15)==gene,]),check.names = F)
ACE2geneexp$group<-ifelse(ACE2geneexp$ENSG00000130234.10 > median(ACE2geneexp$ENSG00000130234.10),"high","low")
ACE2geneexp$group<-as.factor(ACE2geneexp$group)
group<-ACE2geneexp$group
design<-model.matrix(~0+group)
contrast.matrix <- makeContrasts (contrasts = c("grouphigh-grouplow"), levels = design)
fit1 <- contrasts.fit(lmFit(lungTMP, design), contrast.matrix)
fit2 <- eBayes(fit1)
dif <- topTable(fit2,coef="grouphigh-grouplow",n=nrow(fit2),adjust="BH",sort.by = "P")
difG<- AnnoENSG(data=dif)
dataVol<-difG[,c("SYMBOL","logFC","adj.P.Val")]3. 写个画图函数
volcanoPlot.chris<-function(dif,logFCcut=1,padjcut=0.01,selgene="ACE2",col_pal="ChrisLike1"){
colnames(dif)<-c("SYMBOL","logFC","adj.P.Val")
Vol.data<-dif
Vol.data$change <-as.factor(ifelse(Vol.data[,colnames(Vol.data)=="adj.P.Val"] < padjcut & abs(Vol.data[,colnames(Vol.data)=="logFC"] ) > logFCcut,ifelse(Vol.data[,colnames(Vol.data)=="logFC"] > logFCcut,'UP','DOWN'),'NOT'))
if(col_pal=="Classic"){
mypal_V <- pal_gsea("default", n = 3, alpha = 0.9)(3)}
if(col_pal=="ChrisLike1"){
mypal_V <- c("#4E8098","#CED3DC","#A31621")}
if(col_pal=="ChrisLike2"){
mypal_V <- c("#333333","#E1E1E1","#EB593C")
if(col_pal=="ChrisLike3")
mypal_V <- c("#505050","#FFF4C5","#9DD7D5")
p<-ggplot(data = Vol.data, aes(x = logFC,y = -log10(adj.P.Val), color = change, label = SYMBOL)) + geom_point(alpha = 0.8, size = 1) +
theme_bw(base_size = 18) +
geom_vline(xintercept = logFCcut, color = "grey40", linetype = "longdash", size = 0.5) +
geom_vline(xintercept = -logFCcut, color = "grey40", linetype = "longdash", size = 0.5) +
geom_hline(yintercept = -log10(padjcut), color = "grey40", linetype = "longdash", size = 0.5) +
scale_color_manual(name = "", values = mypal_V, limits = c("DOWN", "NOT", "UP"))
pp<-p+geom_label_repel(data = subset(Vol.data,Vol.data$SYMBOL==selgene), aes(label = SYMBOL), show.legend = F)
pp
}4.出图
volcanoPlot.chris(dif=dataVol)
5.换颜色
volcanoPlot.chris(dif=dataVol,col_pal="Classic")
volcanoPlot.chris(dif=dataVol,col_pal="ChrisLike2")
6.换基因
volcanoPlot.chris(dif=dataVol,selgene = "LINC02512")
volcanoPlot.chris(dif=dataVol,selgene = "ARAP1-AS1")
Tips
1、写个函数方便批量组图哦，可以参照下面的教程批量画图
2、关于配色，大家可以修改函数中mypal_V变量的值。
画图素材：
1、在GTEx上下载其中人肺组织表达谱数据
2、需要annoE包
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2020-09-05，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
腾讯云开发者
