本文首发于“生信补给站”公众号 https://mp.weixin.qq.com/s/8XoLkXLbDEDpyof1c8QrPA
单细胞数据完成差异分析后，可以根据结果进行后续的GO ,KEGG,GSEA富集分析，推荐使用clusterProfiler-R包，可参考 R|clusterProfiler-富集分析clusterProfiler|GSEA富集分析及可视化 。
此外还可以进行GSVA分析，基因集变异分析即GSVA（Gene set variation analysis）, 是一种非参数、无监督的分析方法，可以分析 不同的目标基因集 在不同样本中的富集程度。
输入文件比较简单，只需要 A：表达量矩阵 和 B：目标基因集 即可分析。
一 载入R包 数据
1, 获取表达矩阵
如果想计算celltype的GSVA结果，可以使用 AverageExpression 函数计算 不同celltype之间的表达量均值矩阵；
如果计算每个细胞的GSVA结果，直接提取表达量矩阵即可；
library(Seurat) #source("http://bioconductor.org/biocLite.R") #biocLite("GSVA"） library(GSVA) library(tidyverse) library(org.Hs.eg.db) #加载单细胞数据 load("sce.anno.RData") #计算各celltype的表达量均值 Idents(sce2) <- "Anno" expr <- AverageExpression(sce2, assays = "RNA", slot = "data")[[1]] expr <- expr[rowSums(expr)>0,] #过滤细胞表达量全为零的基因 expr <- as.matrix(expr) head(expr)
Epi Fibroblast un T Endo Myeloid OR4F5 2.762724e-04 0.0000000000 0.00000000 0.000000000 0.0000000 0.000000000 AL627309.1 4.399548e-02 0.0080237642 0.01036203 0.007629199 0.0000000 0.047063397 AL627309.5 0.000000e+00 0.0008855168 0.00000000 0.000000000 0.0000000 0.000000000 AP006222.2 1.068307e-01 0.2828329948 0.12706163 0.050051575 0.1942506 0.244211331 LINC01409 1.299390e-04 0.0005598456 0.00000000 0.000000000 0.0000000 0.000000000 FAM87B 7.746596e-05 0.0026794766 0.02404758 0.000000000 0.0000000 0.001711278
单细胞表达量多为0，可以先过滤一下在所有细胞中均无表达的基因。
2，获取目标基因集
根据自己的需要选择MSigDB数据库中的基因集
2.1 手动下载
进入http://www.gsea-msigdb.org/gsea/msigdb/index.jsp后选择需要下载的基因集，然后使用R读取下载好的gmt格式的文件。
2.2 msigdbr包
直接使用msigdbr包内置好的基因集，含有多个物种 以及 多个基因集，通过参数选择物种以及数据集，较为方便。推荐！
library(msigdbr) msigdbr_species() #列出有的物种 #选择基因集合 human_KEGG = msigdbr(species = "Homo sapiens", #物种 category = "C2", subcategory = "KEGG") %>% dplyr::select(gs_name,gene_symbol)#这里可以选择gene symbol或者ID human_KEGG_Set = human_KEGG %>% split(x = .$gene_symbol, f = .$gs_name)#list
# A tibble: 20 × 2 species_name species_common_name <chr> <chr> 1 Anolis carolinensis Carolina anole, green anole 2 Bos taurus bovine, cattle, cow, dairy cow, domestic cattle, domestic cow, ox, oxen 3 Caenorhabditis elegans NA 4 Canis lupus familiaris dog, dogs 5 Danio rerio leopard danio, zebra danio, zebra fish, zebrafish 6 Drosophila melanogaster fruit fly 7 Equus caballus domestic horse, equine, horse 8 Felis catus cat, cats, domestic cat 9 Gallus gallus bantam, chicken, chickens, Gallus domesticus 10 Homo sapiens human 11 Macaca mulatta rhesus macaque, rhesus macaques, Rhesus monkey, rhesus monkeys 12 Monodelphis domestica gray short-tailed opossum 13 Mus musculus house mouse, mouse 14 Ornithorhynchus anatinus duck-billed platypus, duckbill platypus, platypus 15 Pan troglodytes chimpanzee 16 Rattus norvegicus brown rat, Norway rat, rat, rats 17 Saccharomyces cerevisiae baker's yeast, brewer's yeast, S. cerevisiae 18 Schizosaccharomyces pombe 972h- NA 19 Sus scrofa pig, pigs, swine, wild boar 20 Xenopus tropicalis tropical clawed frog, western clawed frog
A：如果你的研究是其中的物种就可以无缝做GSEA 和 GSVA了。
B：如果研究的物种不在其中，也可以自定义基因集，注意转为对应的形式。human_KEGG_Set 为基因集合的列表形式。
二 GSVA分析
1, GSVA分析
数据准备好后，加载GSVA包，一个gsva函数就可以得到GSVA的结果了。
library(GSVA) gsva.kegg <- gsva(expr, gset.idx.list = human_KEGG_Set, kcdf="Gaussian", method = "gsva", parallel.sz=1) head(gsva.kegg)
Epi Myeloid Fibroblast T Endo un KEGG_ABC_TRANSPORTERS -0.30909556 -0.38967397 -0.4080869 -0.36390157 -0.23292000 -0.2954922 KEGG_ACUTE_MYELOID_LEUKEMIA -0.51620884 0.03252670 -0.4041246 0.15837659 0.27282853 -0.3265242 KEGG_ADHERENS_JUNCTION -0.25225184 -0.24621834 -0.1800216 -0.35910806 0.12149392 -0.2252523 KEGG_ADIPOCYTOKINE_SIGNALING_PATHWAY -0.49448870 -0.13413340 -0.4633362 -0.24019340 -0.12668157 -0.3430883 KEGG_ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM 0.03912373 -0.25381877 -0.3541902 -0.03168169 -0.38696763 0.3222462 KEGG_ALDOSTERONE_REGULATED_SODIUM_REABSORPTION -0.28822490 0.07870691 -0.2787129 0.05187873 -0.03042265 0.1007009
行为目标基因集，列为celltype ，数值为gsva分数。
这里需要注意，如果输入矩阵为log转化后的连续表达矩阵指则设置kcdf参数为"Gaussian"，如果是counts矩阵则设置kcdf为"Poisson"。
2, 绘制热图
以结果的前20个绘制示例热图，可以自选择重点的通路
library(pheatmap) pheatmap(gsva.kegg[1:20,], show_colnames = T, scale = "row",angle_col = "45", cluster_row = T,cluster_col = T, color = colorRampPalette(c("navy", "white", "firebrick3"))(50))es = T, scale = "row",angle_col = "45", cluster_row = T,cluster_col = T, color = colorRampPalette(c("navy", "white", "firebrick3"))(50))
3, limma差异分析
和正常的转录组差异分析一样，构建分组信息 以及 比较矩阵，然后使用limma进行差异分析。
此处分析 注释出来的5种celltype 和 注释为unknown之间的差异GSVA结果。
library(limma) # 构建分组文件 group_list <- data.frame(celltype = colnames(gsva.kegg), group = c(rep("con", 5), rep("case", 1))) design <- model.matrix(~ 0 + factor(group_list$group)) colnames(design) <- levels(factor(group_list$group)) rownames(design) <- colnames(gsva.kegg) # 构建差异比较矩阵 contrast.matrix <- makeContrasts(con-case, levels = design) # 差异分析，case vs. con fit <- lmFit(gsva.kegg, design) fit2 <- contrasts.fit(fit, contrast.matrix) fit2 <- eBayes(fit2) diff <- topTable(fit2, coef = 1, n = Inf, adjust.method = "BH", sort.by = "P") head(diff)
logFC AveExpr t P.Value adj.P.Val B KEGG_RIBOSOME -0.7314315 -0.240532302 -2.289899 0.03492937 0.9966699 -4.382757 KEGG_ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM -0.5197531 -0.110881401 -2.100093 0.05077169 0.9966699 -4.418969 KEGG_MATURITY_ONSET_DIABETES_OF_THE_YOUNG -0.5827932 -0.078191940 -1.986597 0.06314026 0.9966699 -4.440295 KEGG_PHENYLALANINE_METABOLISM -0.4828199 -0.125947781 -1.888280 0.07597847 0.9966699 -4.458485 KEGG_ARRHYTHMOGENIC_RIGHT_VENTRICULAR_CARDIOMYOPATHY_ARVC -0.4559099 -0.080461313 -1.826270 0.08522283 0.9966699 -4.469792 KEGG_FOLATE_BIOSYNTHESIS -0.4870146 -0.004848985 -1.811745 0.08752649 0.9966699 -4.472419
4, GSVA差异结果可视化
文章中常见的为双向的柱形图，需要先进行一些设置：
（1）根据logFC 和 adj.P.Val 参数设置，确定上调 ，下调 。为展示结果，以下参数较离谱 ，无参考价值。
实际项目多设置为logFC > 0 & adj.P.Val < 0.05；
#设置分组 diff$group <- ifelse( diff$logFC > 0 & diff$P.Value < 0.3 ,"up" , ifelse(diff$logFC < 0 & diff$P.Value < 0.3 ,"down","noSig") )
（2）设置label的位置
本示例中过滤掉了不显著的通路，在filter行首添加#注释掉即不进行过滤
diff2 <- diff %>% mutate(hjust2 = ifelse(t>0,1,0)) %>% mutate(nudge_y = ifelse(t>0,-0.1,0.1)) %>% filter(group != "noSig") %>% #注释掉该行即可 arrange(t) %>% rownames_to_column("ID")
（3）通过factor设置展示顺序
diff2$ID <- factor(diff2$ID, levels = diff2$ID) limt = max(abs(diff2$t))
（4）ggplot2 可视化展示
ggplot(diff2, aes(ID, t,fill=group)) + geom_bar(stat = 'identity',alpha = 0.7) + scale_fill_manual(breaks=c("down","up"), #设置颜色 values = c("#008020","#08519C"))+ geom_text(data = diff2, aes(label = diff2$ID, #添加通路标签 y = diff2$nudge_y), nudge_x =0,nudge_y =0,hjust =diff2$hjust, size = 3)+ #设置字体大小 labs(x = "KEGG pathways", #设置标题 和 坐标 y=paste0("t value of GSVA score\n","celltype-unknown"), title = "GSVA")+ scale_y_continuous(limits=c(-limt,limt))+ coord_flip() + theme_bw() + #去除背景色 theme(panel.grid =element_blank(), #主题微调 panel.border = element_rect(size = 0.6), plot.title = element_text(hjust = 0.5,size = 18), axis.text.y = element_blank(), axis.title = element_text(hjust = 0.5,size = 18), axis.line = element_blank(), axis.ticks.y = element_blank(), legend.position = "none" #去掉legend )
三 GSVA 样本分组
1, 表达量文件
如果是按照样本分组的话就无需计算每个celltype的表达量均值，直接使用每个细胞的表达量；
expr2 <- as.matrix(sub@assays$RNA@data) gsva.kegg2 <- gsva(expr2, gset.idx.list = keggSet, kcdf="Gaussian",method = "gsva", parallel.sz=1)
2, 分组文件
分组文件因为是每个barcode的粒度，在metadata中构建分组列信息
#之前定义过分组信息 sce2@meta.data$group <- ifelse( grepl("MET",sce2@meta.data$sample ) ,"MET" ,"PT" ) group_list2 <- sce2@meta.data[,c("group")] #标准limma分析 design <- model.matrix(~0+factor(group_list2)) colnames(design)=levels(factor(group_list2)) rownames(design)=colnames(expr) contrast.matrix<-makeContrasts(contrasts = "MET-PT",levels = design) fit <- lmFit(gsva.kegg2,design) fit2 <- contrasts.fit(fit, contrast.matrix) fit2 <- eBayes(fit2) diff2 <- topTable(fit2, coef = 1, n = Inf, adjust.method = "BH", sort.by = "P")
好了，单细胞GSVA分析完成 ，处理数据以及代码都不复杂。
