富集分析-7代码

GSVA 条形图代码V1

代码是参考这里：RNA-seq入门实战（八）：GSVA——基因集变异分析 - 简书 (jianshu.com)

### 条状图



sigdiff_10 <- rbind(subset(sigdiff,logFC>0)[1:10,], subset(sigdiff,logFC<0)[1:10,]) #选择上下调前10通路     

dat_plot <- data.frame(id  = row.names(sigdiff_10),
                       p   = sigdiff_10$P.Value,
                       lgfc= sigdiff_10$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 <- 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)

## 设置不同标签数量
# 小于-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()

# 总数量
high1 <- nrow(dat_plot)
high1

# 绘制条形图
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()
  ) +
  geom_text(data = dat_plot[1:low1,],aes(x = id,y = 0.1,label = id),
            hjust = 0,color = 'black') + #黑色标签
  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') # 黑色标签

ggsave("GSVA_barplot_pvalue.pdf",p,width = 20,height  = 20)


library(tidyverse)  # ggplot2 stringer dplyr tidyr readr purrr  tibble forcats
library(ggthemes)
library(ggprism)

GSVA 条形图代码V2

V2版本代码参考：发散条形图/柱形偏差图 - 简书 (jianshu.com)

############################ 条状图V2

library(tidyverse)  # ggplot2 stringer dplyr tidyr readr purrr  tibble forcats
library(ggthemes)
library(ggprism)




sigdiff_20 <- rbind(subset(sigdiff,logFC>0)[1:20,], subset(sigdiff,logFC<0)[1:20,]) #选择上下调前20通路     

dat_plot <- data.frame(id  = row.names(sigdiff_20),
                       p   = sigdiff_20$P.Value,
                       lgfc= sigdiff_20$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设置为负

# 去掉多余文字 "_UP" 和 "_DN"  
dat_plot$id[1:10]
dat_plot$id <- str_replace(dat_plot$id, "_UP" ,"");dat_plot$id[1:10]
dat_plot$id <- str_replace(dat_plot$id, "_DN" ,"");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


# 接着加上对应的分组标签
## 添加标签
# 小于-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()

# 总数量
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)