R: 相关数据的火山图
2022-09-12 本文已影响0人
胡童远
输入
相关分析
p = c()
cor = c()
for(i in 1:ncol(df_tmp))
{
res = cor.test(df_tmp[,i], index$Matsuda_change, method = c("pearson"))
p = c(p, res$p.value)
cor = c(cor, as.numeric(res$estimate))
}
out = data.frame(species=colnames(df_tmp),p_raw=p,cor_value=cor)
out$p_log10 = -log10(out$p_raw)
write.table(out,file="03_associ/out_correlation_matsuda.txt",
row.names=F, sep="\t", quote=F)
分组
sign=c()
label=c()
for(i in 1:nrow(out))
{
if(out$cor_value[i]<0 & out$p_raw[i]<0.05)
{sign=c(sign,"Negative");label=c(label,out$species[i])}
else if(out$cor_value[i]>0 & out$p_raw[i]<0.05)
{sign=c(sign,"Positive");label=c(label,out$species[i])}
else
{sign=c(sign, "Not significant");label=c(label,"")}
}
out$sign=sign
out$label=label
作图
p=
ggplot(out, aes(x=cor_value,y=p_log10,color=sign)) +
geom_point(size=3) +
labs(x="Correlation value", y="-log10(p value)") +
theme_classic() +
scale_color_manual(
values = c("Positive"="red","Negative"="blue",
"Not significant"="gray")) +
theme(axis.title = element_text(size = 15),
axis.text = element_text(size = 16),
axis.text.x = element_text(hjust = 1),
axis.line = element_line(size = 1),
axis.ticks = element_line(size = 1),
title = element_text(size = 12)) +
geom_text_repel(aes(label = out$label), size = 3)
ggsave(p,file="03_associ/out_correlation_matsuda.pdf",width=8)