1. 火山图？

之前瞥到了小明的跟着Cell学作图：R语言ggplot2作图展示差异表达的基因，但没有细看，根据印象做了下面这个图。

差异表达结果合并

awk 'NR==1{print}FNR>1{print}' *vs0.txt > allvs0.txt

在 R 中读取

library(tidyverse)
df_deg <- read_tsv("allvs0.txt")

作图，其中
不显示padj大于0.5的基因；
df_deg的contrast 列：1vs0 | 2vs0 | 3vs0 | 4vs0 | 5vs0，不同的组的差异表达结果；
df_deg的change列：Up | Down | Stable，可根据log2FoldChange和padj定义；
my_gene 是我想突出的基因构成的向量。

ggplot(filter(df_deg, padj <= 0.05 & !geneids %in% my_genes)) +
  geom_point(aes(contrast, log2FoldChange, color = change, size = -log10(padj), alpha = -log10(padj)),
             # size = 3,
             position = "jitter") +
  geom_point(data = filter(df_deg, (log2FoldChange >= 1 | log2FoldChange <= -1) & padj <= 0.05 & geneids %in% my_genes),
             aes(contrast, log2FoldChange), # size = -log10(padj)
             position = "jitter",
             size = 4,
             shape = 21,
             color = "black",
             fill = "#fee08b") +
  # ylim(-15, 15) +
  geom_hline(yintercept = c(-1, 1), 
             linetype = 'dotdash',
             color = 'grey30') +
  scale_color_manual(values = c(Up = "#e6550d", Down = "#3182bd", Stable = "#636363")) +
  scale_alpha(range = c(0.3, 1)) +
  labs(x = 'Group', y = 'Log2(fold change)') +
  theme_bw()

image.png

2. 常规火山图

只选择一个差异表达结果

library(tidyverse)
library(ggrepel)

df_deg <- read_tsv("1vs0.txt")

ggplot(df_deg) +
  geom_point(aes(log2FoldChange, -log10(padj), color = change),
             show.legend = F) +
  geom_point(data = filter(df_deg, (log2FoldChange >= 1 | log2FoldChange <= -1) & padj <= 0.05 & geneids %in% my_genes),
             aes(log2FoldChange, -log10(padj)),
             # shape = 21,
             # fill = "black",
             color = "#fee08b") +
  geom_text_repel(data = filter(df_deg, (log2FoldChange >= 1 | log2FoldChange <= -1) & padj <= 0.05 & geneids %in% my_genes),
                  aes(log2FoldChange, -log10(padj)),
                  label = "^-^",
                  box.padding = 0.5,
                  max.overlaps = Inf) +
  scale_color_manual(values = c(Up = "#e6550d", Down = "#3182bd", Stable = "#636363")) +
  geom_hline(yintercept = -log10(0.05), 
             linetype = 'dotdash',
             color = 'grey30') +
  geom_vline(xintercept = c(-1, 1), 
             linetype = 'dotdash',
             color = 'grey30') +
  labs(x = 'Log2(fold change)', y = '-log10(p-value)') +
  xlim(-10, 10) +
  ylim(0, 30) +
  theme_bw()