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)
上一篇 下一篇

猜你喜欢

热点阅读