用base plot画火山图
<h1>背景</h1><p>用base plot画美美的火山图</p>
<h1>应用场景</h1>
<p>同时展示多个特征，同时展示P value、ajust P value、分组、size、score等。<br />
例如基因的火山图<br />
<img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/11/20/mi7b6dl7dwo3wc8a8ad.png" alt="" /></p>
<p>出自PMID: 30374049文章</p>
<p>或展示富集分析结果</p>
<p><img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/11/20/mi7b7xjjg40vkad63hs.png" alt="" /></p>
<p>出自PMID: 28957681文章</p>
<h1>环境设置</h1>
<pre><code class="language-{r}">Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor
</code></pre>
<h1>输入文件</h1>
<p>第一列是gene名，后面包含列5列特征值。</p>
<p>可以换成其他值，例如富集分析结果。</p>
<pre><code class="language-{R}">dat <- read.csv("easy_input.csv", header = T, row.names = 1, check.names = F)
head(dat)
dim(dat)
</code></pre>
<h1>把各列数据整理成画图所需的格式</h1>
<pre><code class="language-{r,">### Score列 ###
fc <- dat$score
names(fc) <- rownames(dat)
### -log10P列 ###
p <- dat$`-log10P`
names(p) <- rownames(dat)
### group列 ###
# 给每个pathway的泡泡一种颜色
# 先自定义足够多的颜色
mycol <- c("#B2DF8A","#FB9A99","#33A02C","#E31A1C","#B15928","#6A3D9A","#CAB2D6","#A6CEE3","#1F78B4","#FDBF6F","#999999","#FF7F00")
cols.names <- unique(dat$group)
cols.code <- mycol
names(cols.code) <- cols.names
col <- paste(cols.code,"BB", sep="")
# Highlight your favor
i <- dat$group %in% c("PathwayA","PathwayC","PathwayH","PathwayE")
### size列 ###
sizes <- dat$size
names(sizes) <- rownames(dat)
### pval列 ###
pp <- dat$pval
names(pp) <- rownames(dat)
</code></pre>
<h1>开始画图</h1>
<pre><code class="language-{r}">pdf("base_volcano.pdf", 7, 6)
par(xpd = F, #all plotting is clipped to the plot region
mar = par()$mar + c(0,0,0,6)) #在右侧留出画图例的地方
# base volcano plot
plot(fc, p, log='y',
col=paste(cols.code, "BB", sep=""),
pch=16, #实心圆点
ylab=bquote(~-Log~"P value"), xlab="Enrich score",
cex=ifelse(i, sizes, 1), # 用小泡泡画不感兴趣的pathway
xlim=range(fc * 1.2))
# 添加横线
abline(h=1/0.05, lty=2, lwd=1)
abline(h=1/max(pp), lty=3, lwd=1) #标黑圈和文字的阈值
# 添加竖线
abline(v=-0.5, col="blue", lty=2, lwd=1)
abline(v=0.5, col="red", lty=2, lwd=1)
## 此处用pval列计算adjusted pval来画黑圈和标文字，你还可以另外提供其他信息来画
# 给bonferroni correction pval < 0.001的泡泡标上半透明的文字
w <- which(p.adjust(pp,"bonf") < 0.001) #bonferroni correction
points(fc, p, pch=1, cex=ifelse(i, dat,1))
## Add an alpha value to a colour
add.alpha <- function(col, alpha=1){
if(missing(col))
stop("Please provide a vector of colours.")
apply(sapply(col, col2rgb)/255, 2,
function(x)
rgb(x, x, x, alpha=alpha))
}
cols.alpha <- add.alpha(cols.code$group], alpha=0.6)
text(fc, p, names(fc),
pos=4, #1, 2, 3 and 4, respectively indicate positions below, to the left of, above and to the right of the specified coordinates.
col=cols.alpha)
# 添加size的图例
par(xpd = TRUE) #all plotting is clipped to the figure region
f <- c(0.01,0.05,0.1,0.25)
s <- sqrt(f*50)
legend("topright",
inset=c(-0.2,0), #把图例画到图外
legend=f, pch=16, pt.cex=s, bty='n', col=paste("#88888888"))
# 添加pathway颜色的图例
legend("bottomright",
inset=c(-0.25,0), #把图例画到图外
pch=16, col=cols.code, legend=cols.names, bty="n")
dev.off()
</code></pre>
<p><img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/11/20/mi7bcl9hbqlbuugbx2w.png" alt="" /></p>
<p>文中所用数据可以关注公众号“生信喵实验柴”<br />
<img src="https://tcapi.voiceclouds.cn:8444/uploads/2025/09/07/mf9nu72s90ofyseu9hl.png" alt="" /><br />
发送关键词“20250504”获取</p>
页:
[1]