irGSEA：基于秩次的單細(xì)胞基因集富集分析整合框架（修訂版） - 簡書 (jianshu.com)
GitHub - chuiqin/irGSEA: The integration of single cell rank-based gene set enrichment analysis
計算時間，這個很關(guān)鍵
image.png
數(shù)據(jù)評分
library(Seurat)
obj<- readRDS(file = 'after_annotation.rds')
library(irGSEA)
library(RcppML)
obj<-JoinLayers(obj)
Idents(obj)<-'customclassif'
head(Idents(obj))
obj <- CreateSeuratObject(counts = CreateAssay5Object(GetAssayData(obj, assay = "RNA", slot = "counts")),
meta.data = obj[[]])
obj <- NormalizeData(obj)
obj <- irGSEA.score(object = obj,assay = "RNA",
slot = "data", seeds = 123, ncores = 1,
min.cells = 3, min.feature = 0,
custom = F, geneset = NULL, msigdb = T,
species = "Mus musculus", category = "C2",
subcategory = "CP:KEGG", geneid = "symbol",
method = c("AUCell","UCell","singscore","AddModuleScore"),
aucell.MaxRank = NULL, ucell.MaxRank = NULL,
kcdf = 'Gaussian')
VlnPlot(obj, features = "nCount_AddModuleScore", split.by = "group",
group.by = 'customclassif')
obj$celltype.stim <- paste(obj$customclassif, obj$group, sep = "_")
Idents(obj) <- "celltype.stim"
# 計算差異基因集，并進(jìn)行RRA
# 如果報錯还蹲，考慮加句代碼options(future.globals.maxSize = 100000 * 1024^5)
result.dge <- irGSEA.integrate(object = obj,
group.by = "celltype.stim",
method = c("AUCell","UCell","singscore","AddModuleScore"))
result.dge <- irGSEA.integrate(object = obj,
group.by = NULL,
metadata = obj$celltype.stim,col.name = "ident",
method = c("AUCell","UCell","singscore","AddModuleScore"))
可視化
??irGSEA.heatmap()
irGSEA.heatmap.plot <- irGSEA.heatmap(object = result.dge,
method = "RRA",
top = 50,
show.geneset = NULL,
heatmap.width = 30,
heatmap.heigh = 25)
irGSEA.heatmap.plot
ggsave(plot = irGSEA.heatmap.plot,filename = 'irGSEA.heatmap.plot.pdf',width = 15,height = 12)
# 查看RRA識別的在多種打分方法中都普遍認(rèn)可的差異基因集
geneset.show <- result.dge$RRA %>%
dplyr::filter(pvalue <= 0.05) %>%
dplyr::pull(Name) %>% unique(.)
irGSEA.heatmap.plot1 <- irGSEA.heatmap(object = result.dge,
method = "RRA",
show.geneset = geneset.show)
irGSEA.heatmap.plot1
irGSEA.bubble.plot <- irGSEA.bubble(object = result.dge,
method = "RRA",
show.geneset = NULL)
irGSEA.bubble.plot
ggsave(plot = irGSEA.bubble.plot,filename = 'irGSEA.bubble.plot.pdf',width = 15,height = 12)
irGSEA.upset.plot <- irGSEA.upset(object = result.dge, upset.width = 30,
upset.height = 15,set.size = 1,
method = "RRA")
irGSEA.upset.plot
ggsave(plot = irGSEA.upset.plot,filename = 'irGSEA.upset.plot.pdf',width = 15,height = 12)
irGSEA.barplot.plot <- irGSEA.barplot(object = result.dge,
method = c("AUCell","UCell","singscore",
"AddModuleScore","RRA"))
irGSEA.barplot.plot
ggsave(plot = irGSEA.barplot.plot,filename = 'irGSEA.barplot.plot.pdf',width = 15,height = 12)
aaa<-as.data.frame(geneset.show)
KEGG-ECM-RECEPTOR-INTERACTION
halfvlnplot <- irGSEA.halfvlnplot(object = obj,
method = "AddModuleScore",
show.geneset = "KEGG-ECM-RECEPTOR-INTERACTION")
halfvlnplot
vlnplot <- irGSEA.vlnplot(object = obj,
method = c("AUCell", "UCell", "singscore", "AddModuleScore"),
show.geneset = "KEGG-ECM-RECEPTOR-INTERACTION")
vlnplot
ridgeplot <- irGSEA.ridgeplot(object = obj,
method = "UCell",
show.geneset = "KEGG-ECM-RECEPTOR-INTERACTION")
ridgeplot
#> Picking joint bandwidth of 0.00533
densityheatmap <- irGSEA.densityheatmap(object = obj,
method = "AUCell",
show.geneset = "KEGG-ECM-RECEPTOR-INTERACTION")
densityheatmap
hub.result <- irGSEA.hub(object = obj, assay = "RNA", slot = "data",
method = c("AUCell","UCell","singscore", "AddModuleScore"
),
show.geneset = c("KEGG-ECM-RECEPTOR-INTERACTION",
"KEGG-GLYCEROPHOSPHOLIPID-METABOLISM"),
ncores = 4, type = "rank", maxRank = 2000, top = 5,
correlation.color = c("#0073c2","white","#efc000"),
method.color = NULL)
head(hub.result$hub_result)
head(hub.result$hub_plot$`KEGG-ECM-RECEPTOR-INTERACTION`)
head(hub.result$hub_plot$`KEGG-GLYCEROPHOSPHOLIPID-METABOLISM`)
image.png
R包的安裝
# install packages from CRAN
cran.packages <- c("aplot", "BiocManager", "circlize", "cowplot","data.table",
"devtools", "doParallel", "doRNG", "dplyr", "ggfun", "gghalves",
"ggplot2", "ggplotify", "ggridges", "ggsci", "irlba",
"magrittr", "Matrix", "msigdbr", "pagoda2", "plyr", "pointr",
"purrr", "RcppML", "readr", "reshape2", "reticulate",
"rlang", "RMTstat", "RobustRankAggreg", "roxygen2",
"Seurat", "SeuratObject", "stringr", "tibble", "tidyr",
"tidyselect", "tidytree", "VAM")
for (i in cran.packages) {
if (!requireNamespace(i, quietly = TRUE)) {
install.packages(i, ask = F, update = F)
}
}
# install packages from Bioconductor
bioconductor.packages <- c("AUCell", "BiocParallel", "ComplexHeatmap",
"decoupleR", "fgsea", "ggtree", "GSEABase",
"GSVA", "Nebulosa", "scde", "singscore",
"SummarizedExperiment", "UCell",
"viper","sparseMatrixStats")
for (i in bioconductor.packages) {
if (!requireNamespace(i, quietly = TRUE)) {
BiocManager::install(i, ask = F, update = F)
}
}
# install packages from Github
if (!requireNamespace("irGSEA", quietly = TRUE)) {
devtools::install_github("chuiqin/irGSEA", force =T)
}
#### install packages from Github
# VISION
if (!requireNamespace("VISION", quietly = TRUE)) {
devtools::install_github("YosefLab/VISION", force =T)
}
# mdt need ranger
if (!requireNamespace("ranger", quietly = TRUE)) {
devtools::install_github("imbs-hl/ranger", force =T)
}
# gficf need RcppML (version > 0.3.7) package
if (!utils::packageVersion("RcppML") > "0.3.7") {
message("The version of RcppML should greater than 0.3.7 and install RcppML package from Github")
devtools::install_github("zdebruine/RcppML", force =T)
}
# please first `library(RcppML)` if you want to perform gficf
if (!requireNamespace("gficf", quietly = TRUE)) {
devtools::install_github("gambalab/gficf", force =T)
}
# GSVApy and ssGSEApy need SeuratDisk package
if (!requireNamespace("SeuratDisk", quietly = TRUE)) {
devtools::install_github("mojaveazure/seurat-disk", force =T)
}
# sargent
if (!requireNamespace("sargent", quietly = TRUE)) {
devtools::install_github("Sanofi-Public/PMCB-Sargent", force =T)
}
# pagoda2 need scde package
if (!requireNamespace("scde", quietly = TRUE)) {
devtools::install_github("hms-dbmi/scde", force =T)
}
# if error1 (functio 'sexp_as_cholmod_sparse' not provided by package 'Matrix')
# or error2 (functio 'as_cholmod_sparse' not provided by package 'Matrix') occurs
# when you perform pagoda2, please check the version of irlba and Matrix
# It's ok when I test as follow：
# R 4.2.2 irlba(v 2.3.5.1) Matrix(1.5-3)
# R 4.3.1 irlba(v 2.3.5.1) Matrix(1.6-1.1)
# R 4.3.2 irlba(v 2.3.5.1) Matrix(1.6-3)
#### create conda env
# If error (Unable to find conda binary. Is Anaconda installed) occurs,
# please perform `reticulate::install_miniconda()`
if (! "irGSEA" %in% reticulate::conda_list()$name) {
reticulate::conda_create("irGSEA")
}
# if python package exist
python.package <- reticulate::py_list_packages(envname = "irGSEA")$package
require.package <- c("anndata", "scanpy", "argparse", "gseapy", "decoupler")
for (i in seq_along(require.package)) {
if (i %in% python.package) {
reticulate::conda_install(envname = "irGSEA", packages = i, pip = T)
}
}
部分用戶可以通過鏡像加速安裝
options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
options("repos" = c(CRAN="http://mirrors.cloud.tencent.com/CRAN/"))
# install packages from CRAN
cran.packages <- c("aplot", "BiocManager", "circlize", "cowplot", "data.table",
"devtools", "doParallel", "doRNG", "dplyr", "ggfun", "gghalves",
"ggplot2", "ggplotify", "ggridges", "ggsci", "irlba",
"magrittr", "Matrix", "msigdbr", "pagoda2", "plyr", "pointr",
"purrr", "RcppML", "readr", "reshape2", "reticulate",
"rlang", "RMTstat", "RobustRankAggreg", "roxygen2",
"Seurat", "SeuratObject", "stringr", "tibble", "tidyr",
"tidyselect", "tidytree", "VAM")
for (i in cran.packages) {
if (!requireNamespace(i, quietly = TRUE)) {
install.packages(i, ask = F, update = F)
}
}
# install packages from Bioconductor
bioconductor.packages <- c("AUCell", "BiocParallel", "ComplexHeatmap",
"decoupleR", "fgsea", "ggtree", "GSEABase",
"GSVA", "Nebulosa", "scde", "singscore",
"SummarizedExperiment", "UCell", "viper")
for (i in bioconductor.packages) {
if (!requireNamespace(i, quietly = TRUE)) {
BiocManager::install(i, ask = F, update = F)
}
}
