单细胞基因集打分方法——AUCell

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)

###本地安装AUCell包###
#github包：https://github.com/aertslab/AUCell
devtools::install("../AUCell-master/")
library(AUCell)
#cells_AUC <- AUCell_run(sce.all@assays$RNA@data, geneSets)
cells_rankings <- AUCell_buildRankings(sce.all@assays$RNA@data,splitByBlocks=TRUE) 
cells_rankings
cells_AUC <- AUCell_calcAUC(homo_KEGG_gene, cells_rankings, 
                            aucMaxRank=nrow(cells_rankings)*0.1)
cells_AUC
#提取P53 Pathway
homo_KEGG_gene$KEGG_P53_SIGNALING_PATHWAY
geneSet <- "KEGG_P53_SIGNALING_PATHWAY"
AUCell_auc <- as.numeric(getAUC(cells_AUC)[geneSet, ])
AUCell_auc
#添加至metadata中
sce.all$AUCell <- AUCell_auc
head(sce.all@meta.data)

cells_rankings <- AUCell_buildRankings(sce.all@assays$RNA@data,splitByBlocks=TRUE) 
#Error in AUCell_buildRankings(sce.all@assays$RNA@data, splitByBlocks = TRUE) : 
 # could not find function "AUCell_buildRankings"

cells_rankings
# Ranking for 24590 genes (rows) and 69997 cells (columns).
# 
# Top-left corner of the ranking:
#   cells
# genes           HCC01T_AAACCTGAGGGCATGT HCC01T_AAACCTGAGTCGCCGT
# RP11-34P13.7                    23497                   19215
# FO538757.2                       8795                   20627
# AP006222.2                       1598                    2164
# RP4-669L17.10                   18808                    6330
# RP5-857K21.4                     3508                   17640
# RP11-206L10.9                   19009                   18844
# cells
# genes           HCC01T_AAACCTGGTCACACGC HCC01T_AAACCTGTCCAGTATG
# RP11-34P13.7                    17115                   22521
# FO538757.2                      10350                   11762
# AP006222.2                      16968                   12078
# RP4-669L17.10                   16280                   16136
# RP5-857K21.4                     4398                    8877
# RP11-206L10.9                   22727                   17715
# cells
# genes           HCC01T_AAACGGGTCGCTTGTC
# RP11-34P13.7                     2166
# FO538757.2                      12879
# AP006222.2                      16149
# RP4-669L17.10                   16618
# RP5-857K21.4                    20486
# RP11-206L10.9                   21591

cells_AUC <- AUCell_calcAUC(homo_KEGG_gene, cells_rankings, 
                            aucMaxRank=nrow(cells_rankings)*0.1)