【原】R可視化之美化功能富集條形圖

健明 2021-11-02

展開全文

基因集富集分析是很常見的分析內容，可視化展示的方式也多樣。本文提供兩組分組間的差異表達基因集的功能富集結果的一些相對美觀的可視化方式。

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))