发散条形图/柱形偏差图

该图来源于文章：Phenotype molding of stromal cells in the lung tumor microenvironment | Nature Medicine

• 发散条形图/柱形偏差图展示出了通路上下调信息与t值，最有特点的是将上下调通路放于坐标两侧用不同颜色标注，并将小于阈值的柱形和标签设为灰色，很有意思。尝试参考以下文章绘制这个不知道到底叫发散条形图还是柱形偏差图的图：
  跟着 Nat Med. 学作图 | GSVA+limma差异通路分析+发散条形图 - 简书 (jianshu.com)
  R 绘制柱形偏差图 (qq.com)

由于一般展示通路的显著性都会用Pvalue，t值其实较少使用，因此尝试用Pvalue绘制发散条形图/柱形偏差图，既适用于GSVA结果展示，也适用于差异基因富集分析的通路结果展示。以下数据结果承接RNA-seq入门实战（八）：GSVA——基因集变异分析 中的KEGG的gsva差异分析结果进行绘图。

• 首先对上下调通路进行分组上色绘图放于坐标两侧

#### 发散条形图绘制 ####
library(tidyverse)  # ggplot2 stringer dplyr tidyr readr purrr  tibble forcats
library(ggthemes)
library(ggprism)

degs <- gsva_kegg_degs  #载入gsva的差异分析结果
Diff <- rbind(subset(degs,logFC>0)[1:20,], subset(degs,logFC<0)[1:20,]) #选择上下调前20通路     
dat_plot <- data.frame(id  = row.names(Diff),
                       p   = Diff$P.Value,
                       lgfc= Diff$logFC)
dat_plot$group <- ifelse(dat_plot$lgfc>0 ,1,-1)    # 将上调设为组1，下调设为组-1
dat_plot$lg_p <- -log10(dat_plot$p)*dat_plot$group # 将上调-log10p设置为正，下调-log10p设置为负

# 去掉多余文字
dat_plot$id[1:10]
dat_plot$id <- str_replace(dat_plot$id, "KEGG_","");dat_plot$id[1:10]

# 根据阈值分类
p_cutoff=0.001
dat_plot$threshold <- factor(ifelse(abs(dat_plot$p) <= p_cutoff,
                                   ifelse(dat_plot$lgfc >0 ,'Up','Down'),'Not'),
                            levels=c('Up','Down','Not'))
table(dat_plot$threshold)

# 根据p从小到大排序
dat_plot <- dat_plot %>% arrange(lg_p)
# id变成因子类型
dat_plot$id <- factor(dat_plot$id,levels = dat_plot$id)
# 绘制条形图
p <- ggplot(data = dat_plot,aes(x = id, y = lg_p, 
                                fill = threshold)) +
  geom_col()+
  coord_flip() + #坐标轴旋转
  scale_fill_manual(values = c('Up'= '#36638a','Not'='#cccccc','Down'='#7bcd7b')) +
  geom_hline(yintercept = c(-log10(p_cutoff),log10(p_cutoff)),color = 'white',size = 0.5,lty='dashed') +
  xlab('') + 
  ylab('-log10(P.Value) of GSVA score') + 了
  guides(fill="none")+ # 不显示图例
  theme_prism(border = T) +
  theme(
    plot.margin=unit(c(2,2,2,2),'lines'),#图片四周上右下左间距
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )
p

p

• 接着加上对应的分组标签

## 添加标签
# 小于-cutoff的数量
low1 <- dat_plot %>% filter(lg_p < log10(p_cutoff)) %>% nrow(); low1
# 小于0总数量
low0 <- dat_plot %>% filter(lg_p < 0) %>% nrow(); low0 
# 小于cutoff总数量
high0 <- dat_plot %>% filter(lg_p < -log10(p_cutoff)) %>% nrow(); high0 
# 总数量
high1 <- nrow(dat_plot); high1 

# 依次从下到上添加标签
p1 <- p + geom_text(data = dat_plot[1:low1,],aes(x = id,y = 0.1,label = id),
                   hjust = 0,color = 'black') + # 小于-cutoff的为黑色标签
  geom_text(data = dat_plot[(low1 +1):low0,],aes(x = id,y = 0.1,label = id),
            hjust = 0,color = 'grey') + # 灰色标签
  geom_text(data = dat_plot[(low0 + 1):high0,],aes(x = id,y = -0.1,label = id),
            hjust = 1,color = 'grey') + # 灰色标签
  geom_text(data = dat_plot[(high0 +1):high1,],aes(x = id,y = -0.1,label = id),
            hjust = 1,color = 'black') # 大于cutoff的为黑色标签
p1
ggsave("GSVA_barplot_pvalue.pdf",p1,width = 15,height  = 15)

大功告成：

p1