有了这个包，猪的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 <- T.sc.sc@meta.data[,c("orig.ident", "sex", "age", "stim", "samples")]#分组信息，为了后续作图

#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

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

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)