火山图

2021-8-10

知识来源：R语言实战；训练营

setwd(“D:/R/R-4.0.5/bin/project_visual/visualization”)

rm(list=ls())

准备数据

suppressMessages(library(tidyverse)) 
suppressMessages(library(ggplot2))
suppressMessages(library(ggpubr))
suppressMessages(library(reshape2))
suppressMessages(library(cowplot))
suppressMessages(library(ggrepel))
suppressMessages(library(ggstatsplot))

load("D:/R/R-4.0.5/bin/project_visual/visualization/dif_limma_GSE42872.Rdata")
#处理一下结果，增加gene名和方向
logFC_cutof <- with(DEG_mtx,mean(abs(logFC)) + 2*sd(abs(logFC)))
logFC_cutof <- 0

DEG_mtx$change <- as.factor(ifelse(DEG_mtx$P.Value<0.01 & abs(DEG_mtx$logFC)>logFC_cutof,   ifelse(DEG_mtx$logFC>logFC_cutof,'up','down'),'not'))                         
DEG_mtx$gene <- rownames(DEG_mtx)
head(DEG_mtx)

##             logFC  AveExpr        t      P.Value    adj.P.Val        B change
## DDX3Y    4.256833 8.575889 63.37291 5.671376e-14 1.063610e-09 19.43568     up
## UTY      3.163833 7.632556 43.37164 2.178363e-12 2.042651e-08 17.60272     up
## RPS4Y1   2.869333 7.865222 32.42195 3.548980e-11 2.218586e-07 15.69950     up
## PRUNE2   2.713000 9.497667 29.60698 8.459667e-11 3.768028e-07 15.02565     up
## HLA-DPA1 2.370667 8.736778 29.07918 1.004593e-10 3.768028e-07 14.88817     up
## SETD3    1.473167 9.184111 25.71901 3.241123e-10 1.013067e-06 13.91694     up
##              gene
## DDX3Y       DDX3Y
## UTY           UTY
## RPS4Y1     RPS4Y1
## PRUNE2     PRUNE2
## HLA-DPA1 HLA-DPA1
## SETD3       SETD3

绘图

ggplot(data = DEG_mtx, aes(x = logFC, y = -log10(P.Value),color=change)) +
  #选定数据源和x,y坐标,约定俗成这里用P.value，注意一下
  geom_point(size = 1.5, show.legend = T,alpha = 0.4) +
  #描点，设置大小、颜色分类和图示
  scale_color_manual(values = c("#546de5", "#d2dae2","#ff4757")) + 
  #颜色值
  ylim(0, 12) + 
  #纵坐标的量程
  labs(x = 'Log2(fold change)', y = '-log10(adp)')+ 
  #标签
  geom_hline(yintercept = -log10(0.01), linetype = 'dotdash', color = 'grey20') + 
  #添加水平参考线，p<0.05，dotdash指点虚线
  geom_vline(xintercept = c(-2.5, 2.5), linetype = 'dotdash', color = 'grey20') + 
  #添加垂直参考线，log(fold_change)绝对值大于1
  geom_text_repel(data = subset(DEG_mtx, abs(logFC)>2.5),
                  ###第一步选标签来源数据，利用subset（）选子集
                  size = 2,
                  aes(label = gene))+
  labs(title = "volcano plot")+
  theme(
    legend.title = element_blank()#不显示图例标题
  )+
  theme(plot.title = element_text(size = 16, hjust = 0.5, face = "bold"))+
  #标题的大小、形式(Font face ("plain", "italic", "bold", "bold.italic"))、对齐Horizontal justification (in [0, 1])(0.5居中)
  #theme_set(theme_pubclean())#图片背景
  theme_bw()#修改图片背景

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_text_repel).