AddModuleScore 对单细胞基因集打分及可视化

rm(list=ls())
library(Seurat)
library(ggplot2)
library(dplyr)
library(clusterProfiler)
library(org.Hs.eg.db)
library(stringr)
library(msigdbr)

getwd()
#setwd("./3-cell")  
load("sce.all_by_celltype.Rdata")

#取文章中列出的marker gene
paper<-"ALB,SERPINA1,HNF4A,EPCAM,CD3D,CD3E,NKG7,CD68,CD14,CD163,
CD1C,CLEC4C,KIT, IGHG1,JCHAIN,CD79A,VWF,PECAM1,
FCGR2B,ACTA2,COL1A1,COL1A2"
papermarker<-str_to_upper(trimws(strsplit(paper,',')[[1]]))
papermarker

# 将基因转为list 
features <- list(papermarker)
#直接使用文件中基因向量，并转为list形式
homo_KEGG = msigdbr(species = "Homo sapiens",
                     category = "C2",
                     subcategory = "KEGG") %>% dplyr::select(gs_name,gene_symbol)#这里可以选择gene symbol，也可以选择ID
#基因集是list
homo_KEGG_gene = homo_KEGG %>% split(x =.$gene_symbol, f =.$gs_name)

#选择其中一条通路（我这里选择的是氨基酸和核苷酸糖代谢），将其也转为list
features <- list(homo_KEGG_gene$KEGG_AMINO_SUGAR_AND_NUCLEOTIDE_SUGAR_METABOLISM)

#使用Seurat包自带的AddModuleScore函数打分
sce_score<- AddModuleScore(sce.all,
                       features = features,
                       ctrl = 100,
                       name = "features")
head(sce_score@meta.data)
#这里就得到了基因集评分结果，但是注意列名为 features1
colnames(sce_score@meta.data)[22] <- 'Score'

sce_score@meta.data[1:2,1:22]
#                     orig.ident nCount_RNA nFeature_RNA group percent_mito percent_ribo
#HCC01T_AAACCTGAGGGCATGT   HCC01T       5412         1768     T     3.782703     32.07799
#HCC01T_AAACCTGAGTCGCCGT   HCC01T      12726         2813     T     5.395268     31.67196
#                        percent_hb  S.Score G2M.Score Phase old.ident RNA_snn_res.0.01
#HCC01T_AAACCTGAGGGCATGT 0.00000000 -0.0406769 -0.04946085   G1   HCC01T                0
#HCC01T_AAACCTGAGTCGCCGT 0.02163878 -0.0578437 -0.09095275    G1  HCC01T               2
#                        seurat_clusters RNA_snn_res.0.05 RNA_snn_res.0.1 RNA_snn_res.0.2
#HCC01T_AAACCTGAGGGCATGT               0                0               0               0
#HCC01T_AAACCTGAGTCGCCGT               3                1               1               3
 #                       RNA_snn_res.0.3 RNA_snn_res.0.5 RNA_snn_res.0.8 RNA_snn_res.1
#HCC01T_AAACCTGAGGGCATGT               0               0               0             0
#HCC01T_AAACCTGAGTCGCCGT               3               3               3             3
#                        celltype       Score
#HCC01T_AAACCTGAGGGCATGT     T/NK -0.13832027
#HCC01T_AAACCTGAGTCGCCGT  myeloid -0.06321127

library(ggpubr)
ggviolin(sce_score@meta.data, 
         x="celltype", y="Score", 
         width = 0.8,color = "black",#轮廓颜色
         fill="celltype",#填充 
         xlab = F, #不显示x轴的标签 
         add = 'mean_sd', 
         bxp.errorbar=T,#显示误差条         
         bxp.errorbar.width=0.05, #误差条大小          
         size=0.5, #箱型图边线的粗细
         palette = "npg", 
         legend = "right")
ggsave(filename="violin.pdf",width = 8,height = 5)