双阈值的火山图

标准的火山图：

双阈值的火山图：
耶比较常见咯，比如这篇文章里https://pubmed.ncbi.nlm.nih.gov/32278127/

其实也没啥。就是设置了两组logFC和Pvalue的阈值，把中间的基因用浅一点的颜色表示出来。我试了同一颜色设置不同透明度，结果失败咯。因为点太密集，设置了透明度，颜色耶会叠加，不好看。

1.示例数据

随便拿个芯片数据过来，做完他的差异分析，差异分析结果表格。

if(!require(tinyarray))devtools::install_github("xjsun1221/tinyarray")
library(tinyarray)
library(stringr)
geo = geo_download("GSE4107")
Group = factor(ifelse(str_detect(geo$pd$title,"control"),"control","treat"))
ids = idmap(geo$gpl)
dcp = get_deg_all(log2(geo$exp+1),Group,ids)

## [1] "743 down genes,1908 up genes"

这个10列的表格就是limma差异分析结果，后面添加了几列。

head(dcp$deg)

##      logFC   AveExpr         t      P.Value    adj.P.Val        B    probe_id
## 1 6.861133  7.838697 15.686305 9.489229e-14 2.091036e-09 19.58011   202768_at
## 2 5.219628 11.381501 15.545324 1.147345e-13 2.091036e-09 19.43798   209189_at
## 3 2.453844 13.346921 11.538846 5.019678e-11 3.920727e-07 14.57411 201041_s_at
## 4 4.538726  7.466782 10.617313 2.533831e-10 1.731715e-06 13.19062   223316_at
## 5 3.403667  9.523239 10.411289 3.687011e-10 2.239859e-06 12.86560 201693_s_at
## 6 5.024796  6.872653  9.961547 8.508265e-10 4.651894e-06 12.13533   220276_at
##   symbol change ENTREZID
## 1   FOSB     up     2354
## 2    FOS     up     2353
## 3  DUSP1     up     1843
## 4  CCDC3     up    83643
## 5   EGR1     up     1958
## 6  RERGL     up    79785

2.常规火山图

如果要常规的火山图，已经有啦。我在写这个tinyarray包的时候，顺便把画图的代码也加进去了。仅作简化代码用哦。

dcp$plots

3.双阈值的火山图

每句代码都简单，实际应用起来略复杂，搞起来。

library(ggplot2)
logFC_t1 = 1
P.Value_t1 = 0.05
logFC_t2 = 2
P.Value_t2 = 0.01
dat = dcp$deg
library(dplyr)
k1 = with(dat,logFC > logFC_t2 & P.Value<P.Value_t2);table(k1)

## k1
## FALSE  TRUE 
## 19467   283

k2 = with(dat,logFC < -logFC_t2 & P.Value<P.Value_t2);table(k2)

## k2
## FALSE  TRUE 
## 19692    58

k3 = with(dat,logFC > logFC_t1 & P.Value < P.Value_t1 );table(k3)

## k3
## FALSE  TRUE 
## 17842  1908

k4 = with(dat,logFC < -logFC_t1 & P.Value <P.Value_t1 );table(k4)

## k4
## FALSE  TRUE 
## 19007   743

# 设置不同颜色和大小
my_color = case_when(k1~"red",
                     k2~"blue",
                     k3~"#f08080",
                     k4~"#4169e1",
                     TRUE~"#696969")
my_size = case_when(k1|k2~2.5,
                    k3|k4~2,
                    TRUE~1.5)

p = ggplot(data = dat, 
       aes(x = logFC, 
           y = -log10(P.Value))) +
  geom_point(alpha=0.5, size=my_size, 
             color=my_color) +
  geom_vline(xintercept = c(-logFC_t1,logFC_t1,-logFC_t2,logFC_t2),lty= 4,lwd=0.8,alpha = c(0.5,0.5,1,1)) +
  geom_hline(yintercept = c(-log10(P.Value_t1),-log10(P.Value_t2)),lty= 4,lwd=0.8,alpha = c(0.5,1)) +
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_x_continuous(limits = c(-7, 7), expand = c(0,0))
p

# 加标签
for_label <- dat %>% 
  filter(symbol %in% c("FOSB", "FOS", "DUSP1"))

p +
  geom_point(size = 3, shape = 1, data = for_label) +
  ggrepel::geom_label_repel(
    aes(label = symbol),
    data = for_label,
    color="black"
  )

只要顺利得到了差异分析结果表格，上面的图应该是很轻松搞定咯。如果你仔细看了代码，没太整明白case_when的话，分享下这个函数的核心思想，应该会有茅塞顿开的感觉。可以看到我写的k1 k2 k3 k4 是两组阈值得到差异基因的条件，为什么不是分更多种情况讨论？这是因为case_when函数的条件之间有优先级顺序哦，不满足第一个条件，才往下看第二个条件，以此类推。就像ggplot2的代码也有先后顺序一样！hadly大佬永远的神 ~