|
单细胞数据完成差异分析后，可以根据结果进行后续的GO ,KEGG,GSEA富集分析，推荐使用clusterProfiler-R包，可参考 R|clusterProfiler-富集分析 clusterProfiler|GSEA富集分析及可视化 。 此外还可以进行GSVA分析，基因集变异分析即GSVA（Gene set variation analysis）, 是一种非参数、无监督的分析方法，可以分析 不同的目标基因集 在不同样本中的富集程度。 输入文件比较简单，只需要 A：表达量矩阵 和 B：目标基因集 即可分析。 1, 获取表达矩阵如果想计算celltype的GSVA结果，可以使用 AverageExpression 函数计算 不同celltype之间的表达量均值矩阵； 如果计算每个细胞的GSVA结果，直接提取表达量矩阵即可； library(Seurat) #source("http:///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 MyeloidOR4F5 2.762724e-04 0.0000000000 0.00000000 0.000000000 0.0000000 0.000000000AL627309.1 4.399548e-02 0.0080237642 0.01036203 0.007629199 0.0000000 0.047063397AL627309.5 0.000000e+00 0.0008855168 0.00000000 0.000000000 0.0000000 0.000000000AP006222.2 1.068307e-01 0.2828329948 0.12706163 0.050051575 0.1942506 0.244211331LINC01409 1.299390e-04 0.0005598456 0.00000000 0.000000000 0.0000000 0.000000000FAM87B 7.746596e-05 0.0026794766 0.02404758 0.000000000 0.0000000 0.001711278
单细胞表达量多为0，可以先过滤一下在所有细胞中均无表达的基因。
2，获取目标基因集根据自己的需要选择MSigDB数据库中的基因集 2.1 手动下载 进入http://www./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或者IDhuman_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 为基因集合的列表形式。 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 unKEGG_ABC_TRANSPORTERS -0.30909556 -0.38967397 -0.4080869 -0.36390157 -0.23292000 -0.2954922KEGG_ACUTE_MYELOID_LEUKEMIA -0.51620884 0.03252670 -0.4041246 0.15837659 0.27282853 -0.3265242KEGG_ADHERENS_JUNCTION -0.25225184 -0.24621834 -0.1800216 -0.35910806 0.12149392 -0.2252523KEGG_ADIPOCYTOKINE_SIGNALING_PATHWAY -0.49448870 -0.13413340 -0.4633362 -0.24019340 -0.12668157 -0.3430883KEGG_ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM 0.03912373 -0.25381877 -0.3541902 -0.03168169 -0.38696763 0.3222462KEGG_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))
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. confit <- 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 BKEGG_RIBOSOME -0.7314315 -0.240532302 -2.289899 0.03492937 0.9966699 -4.382757KEGG_ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM -0.5197531 -0.110881401 -2.100093 0.05077169 0.9966699 -4.418969KEGG_MATURITY_ONSET_DIABETES_OF_THE_YOUNG -0.5827932 -0.078191940 -1.986597 0.06314026 0.9966699 -4.440295KEGG_PHENYLALANINE_METABOLISM -0.4828199 -0.125947781 -1.888280 0.07597847 0.9966699 -4.458485KEGG_ARRHYTHMOGENIC_RIGHT_VENTRICULAR_CARDIOMYOPATHY_ARVC -0.4559099 -0.080461313 -1.826270 0.08522283 0.9966699 -4.469792KEGG_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 )
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分析完成 ，处理数据以及代码都不复杂。 精心整理（含图PLUS版）|R语言生信分析，可视化（R统计，ggplot2绘图，生信图形可视化汇总）
|
