当前位置：首页 > news >正文

单细胞scMetabolism代谢相关通路分析学习和整理

news 2025/5/18 19:18:29

scMetabolism一个单细胞水平的代谢相关通路分析工具。

内置了KEGG_metabolism_nc和REACTOME_metabolism两个库的代谢通路信息。

分析方法可选择VISION、AUCell、ssgsea和gsva这四种，默认是VISION。

其他没有什么需要特殊介绍的。

分析流程

1.导入

rm(list = ls())
# V5_path = "/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/seurat5/"
# .libPaths(V5_path)
# .libPaths()
library(stringr)
library(Seurat)
library(rsvd)
library(qs)
library(scMetabolism)
library(BiocParallel)
register(MulticoreParam(workers = 4, progressbar = TRUE))scRNA <- qread("sc_dataset.qs") # V4 data
sub_scRNA <- subset(scRNA,subset = celltype %in% c("epithelial/cancer cells","myeloid cells"))
load("scRNA.Rdata") # V5 data

2.V5版本(需修改代码)

sc.metabolism.SeuratV5 <- function (obj, method = "VISION", imputation = F, ncores = 2, metabolism.type = "KEGG") 
{countexp <- GetAssayData(obj, layer='counts')countexp <- data.frame(as.matrix(countexp))signatures_KEGG_metab <- system.file("data", "KEGG_metabolism_nc.gmt", package = "scMetabolism")signatures_REACTOME_metab <- system.file("data", "REACTOME_metabolism.gmt", package = "scMetabolism")if (metabolism.type == "KEGG") {gmtFile <- signatures_KEGG_metabcat("Your choice is: KEGG\n")}if (metabolism.type == "REACTOME") {gmtFile <- signatures_REACTOME_metabcat("Your choice is: REACTOME\n")}if (imputation == F) {countexp2 <- countexp}if (imputation == T) {cat("Start imputation...\n")cat("Citation: George C. Linderman, Jun Zhao, Yuval Kluger. Zero-preserving imputation of scRNA-seq data using low-rank approximation. bioRxiv. doi: https://doi.org/10.1101/397588 \n")result.completed <- alra(as.matrix(countexp))countexp2 <- result.completed[[3]]row.names(countexp2) <- row.names(countexp)}cat("Start quantify the metabolism activity...\n")if (method == "VISION") {library(VISION)n.umi <- colSums(countexp2)scaled_counts <- t(t(countexp2)/n.umi) * median(n.umi)vis <- Vision(scaled_counts, signatures = gmtFile)options(mc.cores = ncores)vis <- analyze(vis)signature_exp <- data.frame(t(vis@SigScores))}if (method == "AUCell") {library(AUCell)library(GSEABase)cells_rankings <- AUCell_buildRankings(as.matrix(countexp2), nCores = ncores, plotStats = F)geneSets <- getGmt(gmtFile)cells_AUC <- AUCell_calcAUC(geneSets, cells_rankings)signature_exp <- data.frame(getAUC(cells_AUC))}if (method == "ssGSEA") {library(GSVA)library(GSEABase)geneSets <- getGmt(gmtFile)gsva_es <- gsva(as.matrix(countexp2), geneSets, method = c("ssgsea"), kcdf = c("Poisson"), parallel.sz = ncores)signature_exp <- data.frame(gsva_es)}if (method == "ssGSEA") {library(GSVA)library(GSEABase)geneSets <- getGmt(gmtFile)gsva_es <- gsva(as.matrix(countexp2), geneSets, method = c("gsva"), kcdf = c("Poisson"), parallel.sz = ncores)signature_exp <- data.frame(gsva_es)}cat("\nPlease Cite: \nYingcheng Wu, Qiang Gao, et al. Cancer Discovery. 2021. \nhttps://pubmed.ncbi.nlm.nih.gov/34417225/   \n\n")obj@assays$METABOLISM$score <- signature_expobj
}Idents(scRNA) <- "celltype"
res <-sc.metabolism.SeuratV5(obj = scRNA,method = "VISION", # VISION、AUCell、ssgsea和gsvaimputation =F, ncores = 2, metabolism.type = "KEGG") # KEGG和REACTOME

3.V5可视化(需修改代码)

# 需要修改DimPlot.metabolism中的UMAP大小写
DimPlot.metabolismV5 <- function (obj, pathway, dimention.reduction.type = "umap", dimention.reduction.run = T, size = 1) 
{cat("\nPlease Cite: \nYingcheng Wu, Qiang Gao, et al. Cancer Discovery. 2021. \nhttps://pubmed.ncbi.nlm.nih.gov/34417225/   \n\n")if (dimention.reduction.type == "umap") {if (dimention.reduction.run == T) obj <- Seurat::RunUMAP(obj, reduction = "pca", dims = 1:40)umap.loc <- obj@reductions$umap@cell.embeddingsrow.names(umap.loc) <- colnames(obj)signature_exp <- obj@assays$METABOLISM$scoreinput.pathway <- pathwaysignature_ggplot <- data.frame(umap.loc, t(signature_exp[input.pathway, ]))library(wesanderson)pal <- wes_palette("Zissou1", 100, type = "continuous")library(ggplot2)plot <- ggplot(data = signature_ggplot, aes(x = umap_1, y = umap_2, color = signature_ggplot[, 3])) + geom_point(size = size) + scale_fill_gradientn(colours = pal) + scale_color_gradientn(colours = pal) + labs(color = input.pathway) + xlab("UMAP 1") + ylab("UMAP 2") + theme(aspect.ratio = 1) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"))}if (dimention.reduction.type == "tsne") {if (dimention.reduction.run == T) obj <- Seurat::RunTSNE(obj, reduction = "pca", dims = 1:40)tsne.loc <- obj@reductions$tsne@cell.embeddingsrow.names(tsne.loc) <- colnames(obj)signature_exp <- obj@assays$METABOLISM$scoreinput.pathway <- pathwaysignature_ggplot <- data.frame(tsne.loc, t(signature_exp[input.pathway, ]))pal <- wes_palette("Zissou1", 100, type = "continuous")library(ggplot2)plot <- ggplot(data = signature_ggplot, aes(x = tSNE_1, y = tSNE_2, color = signature_ggplot[, 3])) + geom_point(size = size) + scale_fill_gradientn(colours = pal) + scale_color_gradientn(colours = pal) + labs(color = input.pathway) + xlab("tSNE 1") + ylab("tSNE 2") + theme(aspect.ratio = 1) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"))}plot
}# check一下有哪些pathway
pathways <- res@assays$METABOLISM$score
head(rownames(pathways))
# [1] "Glycolysis / Gluconeogenesis"            
# [2] "Citrate cycle (TCA cycle)"               
# [3] "Pentose phosphate pathway"               
# [4] "Pentose and glucuronate interconversions"
# [5] "Fructose and mannose metabolism"         
# [6] "Galactose metabolism"# Dimplot
DimPlot.metabolismV5(obj = res, pathway = "Glycolysis / Gluconeogenesis", dimention.reduction.type = "umap",  dimention.reduction.run = F, size = 1)# DotPlot
input.pathway<-c("Glycolysis / Gluconeogenesis","Oxidative phosphorylation","Citrate cycle (TCA cycle)")
DotPlot.metabolism(obj = res,pathway = input.pathway,phenotype = "celltype", #这个参数需按需修改norm = "y")# BoxPlot
# 由于开发者默认吧obj定义为countexp.Seurat，所以还需要重新命名一下
countexp.Seurat <- res
BoxPlot.metabolism(obj = countexp.Seurat,pathway = input.pathway, phenotype = "celltype", #这个参数需按需修改ncol = 1)

4.V4可直接用

res <-sc.metabolism.Seurat(obj = sub_scRNA,method = "AUCell", # VISION、AUCell、ssgsea和gsvaimputation =F, ncores = 2, metabolism.type = "KEGG") # KEGG和REACTOME

5.V4可视化

# check一下有哪些pathway
pathways <- res@assays$METABOLISM$score
head(rownames(pathways))
# [1] "Terpenoid backbone biosynthesis"                
# [2] "Caffeine metabolism"                            
# [3] "Neomycin, kanamycin and gentamicin biosynthesis"
# [4] "Metabolism of xenobiotics by cytochrome P450"   
# [5] "Drug metabolism - cytochrome P450"              
# [6] "Drug metabolism - other enzymes" # Dimplot绘图
DimPlot.metabolism(obj = res, pathway = "Glycolysis / Gluconeogenesis", dimention.reduction.type = "umap",  dimention.reduction.run = F, size = 1)# DotPlot
input.pathway<-c("Folate biosynthesis","Nicotinate and nicotinamide metabolism","Pyrimidine metabolism")
DotPlot.metabolism(obj = res,pathway = input.pathway,phenotype = "celltype", #这个参数需按需修改norm = "y")#BoxPlot
# 由于开发者默认吧obj定义为countexp.Seurat，所以还需要重新命名一下
countexp.Seurat <- res
BoxPlot.metabolism(obj = countexp.Seurat,pathway = input.pathway, phenotype = "celltype", #这个参数需按需修改ncol = 1)