分享是一种态度
基因集富集分析是很常见的分析内容，可视化展示的方式也多样。本文提供两组分组间的差异表达基因集的功能富集结果的一些相对美观的可视化方式。
1 读取Seurat对象
生成差异表达基因
代码语言：javascript复制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 基因名转换
获取上调和下调基因
代码语言：javascript复制## 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富集分析
代码语言：javascript复制## 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富集分析结果
代码语言：javascript复制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))
