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

2023-10-07 10:22:15 浏览数 (1)

❝上周推文总结了一下使用自定义基因集对单细胞数据进行打分的函数和R包,所以这周使用真实的单细胞数据来做一下可视化展示。

导入数据

❝这里使用的数据是2022年10月20日推文,细胞注释好的数据,原文链接附上。 ❞

文献数据复现——原发性和转移性肝细胞癌多细胞生态系统的单细胞图谱

代码语言:javascript复制
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
可视化,小提琴图
代码语言:javascript复制
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)

在此选择的pathway通路及基因集(基于文章给出的部分基因)是我自己选用,并没有特别的生物学意义,只是做一下可视化展示。

0 人点赞