之前我们发过GSVA分析（有了这个包，猪的GSEA和GSVA分析也不在话下（第一集），【后续来了】有了这个包，猪的GSEA和GSVA分析也不在话下（第二集）），接着单细胞系列，重新说一下GSVA分析。
首先是数据集的问题，通常只做人和小鼠的，想要做其他物种的，苦于没有数据集，不过这里说的这个包，即使猪的GSVA分析也都可以做。
这里我们以小鼠为示例，也是常见物种的GSVA分析，结合单细胞的数据！GSVA其实就是对功能富集的量化，然后进行差异分析，寻找感兴趣的通路在样本中的变化。不同于常规的GO、KEGG受差异基因阈值的影响，GSEA受实验分组的影响，GSVA能够对通路量化，看感兴趣通路在多组之间的变化！
首先加载和安装必要的包并加载单细胞数据。
library(Seurat)#source("/biocLite.R")#biocLite("GSVA"）library(GSVA)library(tidyverse)library(ggplot2)library(clusterProfiler)library(org.Mm.eg.db)library(dplyr)immuneT <- subset(immune, celltype=="T cells")#提取我们需要分析的细胞类型immuneT <- as.matrix(immuneT@assays$RNA@counts)#提取count矩阵meta <- immuneT@meta.data[,c("orig.ident", "sex", "age", "stim", "samples")]#分组信息，为了后续作图
之前一直苦于MSigDB数据库只有人的数据集，没有小鼠和其他物种的，网上也有教程如何根据基因同源性进行转化的，但是很麻烦，也不一定成功。还好有一个新的数据包被发现了，简直是福音---msigdbr包，可以做GSEA和GSVA。
#install.packages("msigdbr")library(msigdbr)msigdbr_species() #可以看到，这个包涵盖了20个物种# A tibble: 20 x 2species_name species_common_name<chr> <chr> 1 Anolis carolinensisCarolina anole, green anole 2 Bos taurus bovine, cattle, cow, dairy cow, domestic cattle, dome~3 Caenorhabditis elegans roundworm4 Canis lupus familiaris dog, dogs5 Danio rerio leopard danio, zebra danio, zebra fish, zebrafish6 Drosophila melanogasterfruit fly7 Equus caballusdomestic horse, equine, horse8 Felis catus cat, cats, domestic cat 9 Gallus gallusbantam, chicken, chickens, Gallus domesticus10 Homo sapiens human 11 Macaca mulattarhesus macaque, rhesus macaques, Rhesus monkey, rhesu~12 Monodelphis domestica gray short-tailed opossum 13 Mus musculus house mouse, mouse14 Ornithorhynchus anatinusduck-billed platypus, duckbill platypus, platypus15 Pan troglodytes chimpanzee 16 Rattus norvegicus brown rat, Norway rat, rat, rats 17 Saccharomyces cerevisiaebaker's yeast, brewer's yeast, S. cerevisiae18 Schizosaccharomyces pombe 9~ NA 19 Sus scrofa pig, pigs, swine, wild boar 20 Xenopus tropicalis tropical clawed frog, western clawed frog 查看下小鼠的基因集，是否与MSigDB数据库一样mouse <- msigdbr(species = "Mus musculus")mouse[1:5,1:5]# A tibble: 5 x 5gs_cat gs_subcatgs_name gene_symbol entrez_gene<chr> <chr><chr><chr> <int>1 C3MIR:MIR_Legacy AAACCAC_MIR140 Abcc4 2392732 C3MIR:MIR_Legacy AAACCAC_MIR140 Abraxas2 1093593 C3MIR:MIR_Legacy AAACCAC_MIR140 Actn4 605954 C3MIR:MIR_Legacy AAACCAC_MIR140 Acvr1 114775 C3MIR:MIR_Legacy AAACCAC_MIR140 Adam9 11502table(mouse$gs_cat) #查看目录，与MSigDB一样，包含9个数据集###C1C2C3C4C5C6C7C8 H 9 533767 795972 92353 1248327 30556 988358 109328 7411
本例中，我们要分析GO，因为mouse文件包含了所有的基因集，所以要查看GO在哪里，然后将需要的文件提出来。
table(mouse$gs_subcat)CGN CGP CM CP 167344 42770376981 49583 3881 CP:BIOCARTA CP:KEGGCP:PIDCP:REACTOME CP:WIKIPATHWAYS 4860 13694 8196 98232 27923 GO:BP GO:CC GO:MF HPOIMMUNESIGDB 660368100991105717381251944068 MIR:MIR_Legacy MIR:MIRDB TFT:GTRD TFT:TFT_Legacy VAX 34118372658235886153310 44290 mouse_GO_bp = msigdbr(species = "Mus musculus",category = "C5", #GO在C5subcategory = "GO:BP") %>% dplyr::select(gs_name,gene_symbol)#这里可以选择gene symbol，也可以选择ID，根据自己数据需求来，主要为了方便mouse_GO_bp_Set = mouse_GO_bp %>% split(x = .$gene_symbol, f = .$gs_name)#后续gsva要求是list，所以将他转化为list
表达矩阵数据有了，通路信息有了，就可以进行GEVA分析了，代码简单就一句！保存结果！
T_gsva <- gsva(expr = immuneT, gset.idx.list = mouse_GO_bp_Set,kcdf="Poisson", #查看帮助函数选择合适的kcdf方法 parallel.sz = 5)write.table(T_gsva, 'T_gsva.xls', row.names=T, col.names=NA, sep="\t")
接着差异分析可以用limma包，类似于转录组芯片数据分析流程。
group <- c(rep("control", 50), rep("test", 71)) %>% as.factor()#设置分组，对照在前desigN <- model.matrix(~ 0 + group) #构建比较矩阵colnames(desigN) <- levels(group)fit = lmFit(test_control, desigN)fit2 <- eBayes(fit)diff=topTable(fit2,adjust='fdr',coef=2,number=Inf)#校准按照需求修改 ？topTablewrite.csv(diff, file = "Diff.csv")
最后对差异的感兴趣的通路进行可视化：
up <- c("GOBP_EGG_ACTIVATION","GOBP_TENDON_DEVELOPMENT","GOBP_SOMITE_SPECIFICATION","GOBP_THREONINE_CATABOLIC_PROCESS","GOBP_REGULATION_OF_GLUTAMATE_RECEPTOR_CLUSTERING","GOBP_NEGATIVE_CHEMOTAXIS","GOBP_NEGATIVE_REGULATION_OF_FAT_CELL_PROLIFERATION","GOBP_REGULATION_OF_T_HELPER_17_CELL_LINEAGE_COMMITMENT","GOBP_REGULATION_OF_ANTIMICROBIAL_HUMORAL_RESPONSE")down <- c("GOBP_DETERMINATION_OF_PANCREATIC_LEFT_RIGHT_ASYMMETRY","GOBP_MITOTIC_DNA_REPLICATION","GOBP_EOSINOPHIL_CHEMOTAXIS","GOBP_NEUTROPHIL_MEDIATED_CYTOTOXICITY","GOBP_POTASSIUM_ION_EXPORT_ACROSS_PLASMA_MEMBRANE","GOBP_REGULATION_OF_LEUKOCYTE_MEDIATED_CYTOTOXICITY","GOBP_REGULATION_OF_SEQUESTERING_OF_ZINC_ION","GOBP_ENDOTHELIN_RECEPTOR_SIGNALING_PATHWAY","GOBP_PRE_REPLICATIVE_COMPLEX_ASSEMBLY_INVOLVED_IN_CELL_CYCLE_DNA_REPLICATION","GOBP_ESTABLISHMENT_OF_PLANAR_POLARITY_OF_EMBRYONIC_EPITHELIUM")TEST <- c(up,down)diff$ID <- rownames(diff) Q <- diff[TEST,]group1 <- c(rep("treat", 9), rep("control", 10)) df <- data.frame(ID = Q$ID, score = Q$t,group=group1 )# 按照t score排序sortdf <- df[order(df$score),]sortdf$ID <- factor(sortdf$ID, levels = sortdf$ID)#增加通路ID那一列ggplot(sortdf, aes(ID, score,fill=group)) + geom_bar(stat = 'identity',alpha = 0.7) + coord_flip() + theme_bw() + #去除背景色theme(panel.grid =element_blank())+theme(panel.border = element_rect(size = 0.6))+labs(x = "",y="t value of GSVA score")+scale_fill_manual(values = c("#008020","#08519C"))#设置颜色
整个流程就结束了，希望对你们的研究能有启发，GO、KEGG做多了，可以换着做一下GSVA分析！
