R可视化之美化功能富集条形图

基因集富集分析是很常见的分析内容，可视化展示的方式也多样。本文提供两组分组间的差异表达基因集的功能富集结果的一些相对美观的可视化方式。

1 读取Seurat对象

生成差异表达基因

library(Seurat)
library(dplyr)
sce <- readRDS('F:/R_Language/R_Practice/scRNA_Seq_column/src/scRNA-seq_advance/Data/sce.rds')

## 1 设置分组
sce.all = sce
table(Idents(sce.all))
## 
##  Naive CD4 T   CD14+ Mono Memory CD4 T            B        CD8 T FCGR3A+ Mono 
##          711          480          472          344          279          162 
##           NK           DC     Platelet 
##          144           32           14
two_groups <- c("Naive CD4 T", "Memory CD4 T")
sce.sub = sce.all[, Idents(sce.all) %in% two_groups] # 挑选细胞

## 2 生成差异表达基因
sce.markers = FindMarkers(object = sce.sub, ident.1 = two_groups[1], ident.2 = two_groups[2], test.use='MAST' )  ## MAST在单细胞领域较为常用
sce.markers <- sce.markers %>% tibble::rownames_to_column('gene')
head(sce.markers)
##      gene        p_val avg_log2FC pct.1 pct.2    p_val_adj
## 1  S100A4 2.616570e-71 -1.3010105 0.686 0.951 3.588363e-67
## 2     B2M 2.251490e-40 -0.3812789 1.000 1.000 3.087693e-36
## 3 S100A11 7.402574e-35 -1.0582548 0.271 0.633 1.015189e-30
## 4   ANXA1 7.945638e-35 -1.0066778 0.482 0.805 1.089665e-30
## 5    IL32 1.321968e-31 -0.7828632 0.768 0.949 1.812947e-27
## 6  MALAT1 1.848406e-31  0.4423242 1.000 1.000 2.534905e-27


2 基因名转换

获取上调和下调基因

## 3-1 Load packages
library(ggpubr)
library(clusterProfiler)
options(connectionObserver = NULL) #加载org.Hs.eg.db失败时的解决方法
options(stringsAsFactors = F)
suppressMessages(library('org.Hs.eg.db'))
keytypes(org.Hs.eg.db)
##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
##  [6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"    
## [11] "GO"           "GOALL"        "IPI"          "MAP"          "OMIM"        
## [16] "ONTOLOGY"     "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
## [21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"       "UNIGENE"     
## [26] "UNIPROT"

## 3-2 Specify species name
go_organism <- "org.Hs.eg.db" # org.Hs.eg.db/org.Mm.eg.db
kegg_organism <- 'hsa' # hsa/mmu

ids = bitr(sce.markers$gene, fromType="SYMBOL", toType="ENTREZID", OrgDb=go_organism)
sce.markers = merge(sce.markers, ids, by.x='gene', by.y='SYMBOL')
sce.markers$group <- factor(ifelse(sce.markers$avg_log2FC < 0, -1, 1), levels = c(-1, 1))
gcSample = split(sce.markers$ENTREZID, sce.markers$group)
formula_res <- compareCluster(gcSample, fun="enrichKEGG", organism = kegg_organism, pvalueCutoff=0.05)
dotplot(formula_res) 


## 获取上下调基因
up_genes <- subset(sce.markers, group==1)$ENTREZID
down_genes <- subset(sce.markers, group==-1)$ENTREZID


3 KEGG富集分析

## KEGG Annotation
##########################
# FunctionName函数名首字母大写, 不用点分隔 (所含单词首字母大写)，函数命名应为动词或动词性短语。eg: CalculateAvgClicks
EnrichGeneKegg <- function(gene_set){
  kegg_anno <- enrichKEGG(gene_set, organism = kegg_organism, keyType = 'kegg', pvalueCutoff = 0.05,pAdjustMethod = 'BH',
                        minGSSize = 10,maxGSSize = 500, qvalueCutoff = 0.2,use_internal_data = FALSE)
  kegg_anno <- as.data.frame(kegg_anno)
  return(kegg_anno)
}

up_kegg <- EnrichGeneKegg(gene_set=up_genes) %>% dplyr::mutate(group=rep(1, n()))
down_kegg <- EnrichGeneKegg(gene_set=down_genes) %>% dplyr::mutate(group=rep(-1, n()))
dat <- rbind(down_kegg, up_kegg)

# Rename group information
sample_names <- c(paste0(two_groups[1]," up-regulated"), paste0(two_groups[2]," up-regulated"))
dat <- dat %>% 
  dplyr::mutate(group_type = factor(ifelse(group == 1, sample_names[1], sample_names[2]), levels = sample_names)) %>% 
  dplyr::mutate(Gene_Number = Count * group)

# 为便于图形展示，提取部分子集
dat <- dat %>% dplyr::group_by(group_type) %>% dplyr::do(head(., n = 5)) 


4 可视化KEGG富集分析结果

library(ggplot2)
library(cowplot)
## 设置统一的主题
th <- theme(axis.text.y = element_blank(), 
            axis.ticks.y = element_blank(),
            axis.title.y = element_blank(),
            panel.background = element_rect(fill = 'white'), 
            panel.border = element_rect(color = 'black', fill = NA),
            panel.grid.major.y = element_blank(),
            panel.grid.minor.y = element_blank(),
            panel.grid.major.x = element_line(colour = "grey80", linetype = "dashed"),
            panel.grid.minor.x = element_blank())

## 图一：常规分组条形图
ggplot(dat, aes(x = reorder(Description, Gene_Number), y = Gene_Number)) +
  geom_bar(aes(fill = group_type), stat = "identity") + 
  labs(x = '', y = 'Gene Number', fill = '') +
  coord_flip() + theme(legend.position = 'bottom') + guides(fill = guide_legend(ncol = 1))


## 图二：当功能名写在条形里
ggplot(dat, aes(x = reorder(Description, Gene_Number), y = Gene_Number)) +
  geom_bar(aes(fill = group_type), stat = "identity") + 
  labs(y = 'Gene Number', fill = '') +
  coord_flip() + 
  geom_text(aes(y = 0, label = Description, hjust = as.numeric(Gene_Number < 0))) +  # label text based on value
  th + theme(legend.position = c(0.25, 0.9)) + guides(fill = guide_legend(ncol = 1))


## 图三：当功能名写在条形外
ggplot(dat, aes(x = reorder(Description, Gene_Number), y = Gene_Number)) +
  geom_bar(aes(fill = group_type), stat = "identity") + 
  labs(y = 'Gene Number', fill = '') +
  coord_flip() + 
  geom_text(aes(y = 0, label = Description, hjust = as.numeric(Gene_Number > 0))) +  # label text based on value
  th + theme(legend.position = c(0.25, 0.9)) + guides(fill = guide_legend(ncol = 1))


## 图四：用条形颜色的深浅和长度同时表达数据
nbreaks = 10
minimum <- floor(min(dat$Gene_Number)/nbreaks)*nbreaks
maximum <- ceiling(max(dat$Gene_Number)/nbreaks)*nbreaks

ggplot(dat, aes(x = reorder(Description, Gene_Number), y = Gene_Number)) +
  geom_bar(aes(fill = Gene_Number), stat = "identity") + 
  labs(y = 'Gene Number', fill = '') +
  coord_flip() +  
  geom_text(aes(y = 0, label = Description, hjust = as.numeric(Gene_Number > 0))) +  # label text based on value
  th + 
  scale_y_continuous(limits = c(minimum, maximum), breaks = seq(minimum, maximum, nbreaks)) +
  scale_fill_gradient2(low = 'darkblue', high = 'red', mid = 'white',
                       limits = c(minimum, maximum), breaks = seq(minimum, maximum, nbreaks)) #修改图例名字以及图中颜色


## 图五：借用ggpubr包绘制条形图
library(ggpubr)
dat2 <- dat %>% 
  dplyr::group_by(group_type) %>%
  dplyr::arrange(Gene_Number)

ggbarplot(dat2, x = "Description", y = "Gene_Number",
          fill = "group_type",           # change fill color by mpg_level
          color = "white",            # Set bar border colors to white
          palette = "jco",            # jco journal color palett. see ?ggpar
          # sort.val = "desc",          # Sort the value in descending order
          sort.by.groups = F,     # Don't sort inside each group
          # x.text.angle = 90,          # Rotate vertically x axis texts
          ylab = "Gene Number",
          legend.title = "",
          rotate = TRUE,
          ggtheme = theme_classic()) + 
            theme(legend.position = 'bottom') + guides(fill = guide_legend(ncol = 1))