有了这个包，猪的GSEA和GSVA分析也不在话下（第一集）

救救孩子吧，GSVA分析都是做人的，有现成的人的数据集，可是其他物种的就惨了，很难下手！
今天我们就说说小鼠，也是常见物种的GSVA分析，结合单细胞的数据！

GSVA的作用不用多说了，大家都熟悉，至少选择要做这个分析，那么他的作用也是清楚的。其实就是对功能富集的量化，然后进行差异分析，寻找感兴趣的通路在样本中的变化。不同于常规的GO、KEGG受差异基因阈值的影响，GSEA受实验分组的影响，GSVA能够对通路量化，看感兴趣通路在多组之间的变化！

首先加载和安装必要的包

library(Seurat)#source("http://bioconductor.org/biocLite.R")#biocLite("GSVA"）library(GSVA)library(tidyverse)library(ggplot2)library(clusterProfiler)library(org.Mm.eg.db)library(dplyr)

导入单细胞转录组数据：

load("D:/SC_re.RData") #加载数据集T.sc <- subset(SC_re, celltype=="T cells")#用subset函数提取我们需要分析的细胞类型T.exp <- as.matrix(Rods.sc@assays$RNA@counts)#提取count矩阵，一定要是as.matrix，如果开始时是dataframe，后期GSVA分析时也要转为matrixmeta <- [e-mail][,c("orig.ident", "sex", "age", "stim", "samples")]#分组信息，为了后续作图

这里有个重要的问题就是，单细胞GSVA到底用什么数据，count or data？之前《单细胞天地》公众号依据测试过这个问题了，count数据和data数据做GSVA分析没有区别，但是不能使用高变基因去做，对结果的影响较大，还是使用原始的count 或者 data数据吧。

==================重点来了==================

之前一直苦于MSigDB数据库只有人的数据集，没有小鼠和其他物种的，网上也有教程如何根据基因同源性进行转化的，但是很麻烦，也不一定成功。还好有一个新的数据包被发现了，简直是福音---msigdbr包，可以做GSEA和GSVA。

#install.packages("msigdbr")library(msigdbr)msigdbr_species() #可以看到，这个包涵盖了20个物种# A tibble: 20 x 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, dome~ 3 Caenorhabditis elegans roundworm 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, rhesu~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 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 5 gs_cat gs_subcat gs_name gene_symbol entrez_gene <chr> <chr> <chr> <chr> <int>1 C3 MIR:MIR_Legacy AAACCAC_MIR140 Abcc4 2392732 C3 MIR:MIR_Legacy AAACCAC_MIR140 Abraxas2 1093593 C3 MIR:MIR_Legacy AAACCAC_MIR140 Actn4 605954 C3 MIR:MIR_Legacy AAACCAC_MIR140 Acvr1 114775 C3 MIR:MIR_Legacy AAACCAC_MIR140 Adam9 11502table(mouse$gs_cat) #查看目录，与MSigDB一样，包含9个数据集###C1 C2 C3 C4 C5 C6 C7 C8 H 20049 533767 795972 92353 1248327 30556 988358 109328 7411

例如，本例中，我们要分析GO，因为mouse文件包含了所有的基因集，所以要查看GO在哪里，然后将需要的文件提出来。

table(mouse$gs_subcat) CGN CGP CM CP 167344 42770 376981 49583 3881 CP:BIOCARTA CP:KEGG CP:PID CP:REACTOME CP:WIKIPATHWAYS 4860 13694 8196 98232 27923 GO:BP GO:CC GO:MF HPO IMMUNESIGDB 660368 100991 105717 381251 944068 MIR:MIR_Legacy MIR:MIRDB TFT:GTRD TFT:TFT_Legacy VAX 34118 372658 235886 153310 44290 mouse_GO_bp = msigdbr(species = "Mus musculus", category = "C5", #GO在C5 subcategory = "GO:BP") %>% dplyr::select(gs_name,gene_symbol)#这里可以选择gene symbol，也可以选择ID，根据自己数据需求来，主要为了方便head(mouse_GO_bp)# A tibble: 6 x 2 gs_name gene_symbol <chr> <chr> 1 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS Aldh1l1 2 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS Aldh1l2 3 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS Mthfd1 4 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS Mthfd1l 5 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS Mthfd2l 6 GOBP_2_OXOGLUTARATE_METABOLIC_PROCESS Aadatmouse_GO_bp_Set = mouse_GO_bp %>% split(x = .$gene_symbol, f = .$gs_name)#后续gsva要求是list，所以将他转化为list

表达矩阵数据有了，通路信息有了，就可以进行GEVA分析了，代码简单就一句！
提示：数据大的话这一步时间较长，出去吃个饭吧（正好午饭时间）！

T_gsva <- gsva(expr = T.exp, gset.idx.list = mouse_GO_bp_Set, kcdf="Poisson", #查看帮助函数选择合适的kcdf方法 parallel.sz = 5) #Setting parallel calculations through a MulticoreParam back-end#with workers=5 and tasks=100.#Estimating GSVA scores for 7410 gene sets.#Estimating ECDFs with Poisson kernels#Estimating ECDFs in parallel#iteration: 100#|=============================================================================| 100%#查看分析结果：变成了GOBP的表达值了head(MG_gsva[1:4, 1:4])#GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS 0.3971866#GOBP_2_OXOGLUTARATE_METABOLIC_PROCESS -0.3306133#GOBP_2FE_2S_CLUSTER_ASSEMBLY -0.2476997#GOBP_3_PHOSPHOADENOSINE_5_PHOSPHOSULFATE_BIOSYNTHETIC_PROCESS 0.1370572

记得将结果保存一下

write.table(T_gsva, 'T_gsva.xls', row.names=T, col.names=NA, sep="\t")

先画个热图展示一下结果吧：

library(pheatmap)library(patchwork)A <- MG_gsva[1:50,] #为了方便展示，我们只展示前50行pheatmap(A, show_rownames=1, show_colnames=0)

结果如下：

之后的差异分析，可视化我们下节再说！

记得点赞加关注，下节不迷路！