【原】scRNA|使用scMetabolism完成單細(xì)胞代謝激活分?jǐn)?shù)估計(jì)

生信補(bǔ)給站 2024-04-24 發(fā)布于北京

展開全文

之前介紹過 scRNA分析|使用AddModuleScore 和 AUcell進(jìn)行基因集打分，然后可視化目標(biāo)基因集合的打分，這里介紹scMetabolism包-整合了多個可以完成細(xì)胞代謝相關(guān)通路評估方法的R包。

一載入R包，數(shù)據(jù)

首先根據(jù)官網(wǎng)GitHub - wu-yc/scMetabolism: Quantifying metabolism activity at the single-cell resolution的介紹安裝相關(guān)的R包，需要注意的是VISION要安裝v2.1.0版本。

然后使用之前注釋過的sce.anno.RData數(shù)據(jù) ，為節(jié)省資源，每種細(xì)胞類型隨機(jī)抽取30%的數(shù)據(jù)。

install.packages(c("devtools", "data.table", "wesanderson", "Seurat", "devtools",                    "AUCell", "GSEABase", "GSVA", "ggplot2","rsvd"))#Please note that the version would be v2.1.0devtools::install_github("YosefLab/VISION@v2.1.0") devtools::install_github("wu-yc/scMetabolism")# 加載R包library(scMetabolism)library(tidyverse)library(rsvd)library(Seurat)library(pheatmap)library(ComplexHeatmap)library(ggsci)# 加載數(shù)據(jù)load("sce.anno.RData")sce2@meta.data$CB <- rownames(sce2@meta.data)# 按照細(xì)胞類型抽取一定比例的數(shù)據(jù)sample_CB <- sce2@meta.data %>%   group_by(celltype) %>%   sample_frac(0.3)# 提取數(shù)據(jù)sce3 <- subset(sce2,CB %in% sample_CB$CB) head(sce3，2)

二計(jì)算代謝得分

該包比較簡單，主函數(shù)可以選擇sc.metabolism.Seurat 輸入Seurat的單細(xì)胞對象（推薦），也可以選擇 sc.metabolism 輸入矩陣（作者不太建議）。

Idents(sce3) <- "celltype"countexp.Seurat <- sc.metabolism.Seurat(obj = sce3,  #Seuratde單細(xì)胞object                                      method = "AUCell",                                       imputation = F,                                       ncores = 2,                                       metabolism.type = "KEGG")

其中obj是一個包含 UMI 計(jì)數(shù)矩陣的 Seurat 對象，記得指定Idents 。

method支持VISION、AUCell、ssgsea和gsva四種，默認(rèn)的是VISION 方法。

metabolism.type支持KEGG和REACTOME，分別對應(yīng)不同的代謝相關(guān)通路。

1，查看函數(shù)

可以用過View(sc.metabolism.Seurat) 查看函數(shù)的主體，結(jié)構(gòu)還是比較清楚的，（1）預(yù)設(shè)了KEGG和REACTOME中代謝相關(guān)通路，（2）根據(jù)VISION、AUCell、ssgsea和gsva 四種常見方法計(jì)算代謝通路相關(guān)的得分。

注：gmt可以改為你課題需要的通路，然后放到signatures_KEGG_metab輸出的路徑下。

也可以如 scRNA分析|使用AddModuleScore 和 AUcell進(jìn)行基因集打分使用相關(guān)方法的函數(shù)直接計(jì)算。

signatures_KEGG_metab <- system.file("data", "KEGG_metabolism_nc.gmt",                                       package = "scMetabolism")signatures_KEGG_metab#[1] "C:/Users/XXX/AppData/Local/R/win-library/4.3/scMetabolism/data/KEGG_metabolism_nc.gmt"

2，提取結(jié)果-添加至meta信息

代謝評分的結(jié)果存放在新的assays -- METABOLISM中，可以通過如下方式得到每個基因的代謝通路的活性分?jǐn)?shù)。

如截圖所示細(xì)胞barcode的"-1"變?yōu)榱?.1"，通過str_replace_all簡單處理后添加至meta中，以備后面可能的相關(guān)分析。

#提取score結(jié)果score <- countexp.Seurat@assays$METABOLISM$scorescore[1:4,1:4]#將score中barcode的點(diǎn)轉(zhuǎn)為下劃線score_change <- score %>%   select_all(~str_replace_all(., "\\.", "-"))  #基因ID不規(guī)范會報(bào)錯,下劃線替換-#確定細(xì)胞barcode椅子identical(colnames(score_change) , rownames(countexp.Seurat@meta.data))#[1] TRUEcountexp.Seurat@meta.data <- cbind(countexp.Seurat@meta.data,t(score_change) )
#可以直接使用Seurat的相關(guān)函數(shù)p1 <- FeaturePlot(countexp.Seurat,features = "Glycolysis / Gluconeogenesis")p2 <- VlnPlot(countexp.Seurat,features = "Glycolysis / Gluconeogenesis")p1 + p2

三可視化

可以使用scMetabolism自帶的函數(shù)完成一些可視化展示。

1，umap展示某條通路的代謝得分

DimPlot.metabolism(obj = countexp.Seurat,                    pathway = "Glycolysis / Gluconeogenesis",                    dimention.reduction.type = "umap",                    dimention.reduction.run = F, size = 1)

2，指定通路-細(xì)胞類型點(diǎn)圖

可以選擇直接指定目標(biāo)通路或者展示前幾個，注意將phenotype 參數(shù)改為需要展示的列。

#直接指定input.pathway<-c("Glycolysis / Gluconeogenesis",                  "Oxidative phosphorylation",                  "Citrate cycle (TCA cycle)")#展示前10個input.pathway <- rownames(countexp.Seurat@assays$METABOLISM$score)[1:10]
DotPlot.metabolism(obj = countexp.Seurat,                    pathway = input.pathway,                    phenotype = "celltype",  #更改phenotype 參數(shù)                   norm = "y")

3，指定通路-箱線圖

可以使用ggsci 包修改一下顏色

BoxPlot.metabolism(obj = countexp.Seurat,                    pathway = input.pathway[1:4],                    phenotype = "celltype",                    ncol = 2) +  scale_fill_nejm()

4，自定義熱圖

首先計(jì)算每種細(xì)胞類型的相關(guān)代謝通路得分的均值，然后可以使用pheatmap 直接繪制熱圖，或者參照scRNA｜ComplexHeatmap自定義單細(xì)胞轉(zhuǎn)錄組celltype-level 熱圖可視化繪制復(fù)雜熱圖

#可以計(jì)算celltype均值，然后繪制df <- countexp.Seurat@meta.data#19列開始是代謝通路的得分，按照celltype計(jì)算均值avg_df = aggregate(df[,19:ncol(df)],                   list(df$celltype),                   mean)#熱圖需要轉(zhuǎn)為矩陣avg_df <- avg_df %>%   select(1:20) %>% #展示前20個  column_to_rownames("Group.1") avg_df[1:4,1:4]

也可以手動選擇想展示的代謝通路。

（1）直接pheatmap繪制

pheatmap(t(avg_df),          show_colnames = T,         scale='row',          cluster_rows = T,         color=colorRampPalette(c('#1A5592','white',"#B83D3D"))(100),         cluster_cols = T)

（2）組合復(fù)雜熱圖

作為復(fù)雜熱圖的一個組件。為使圖形更好看，我們先手動對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化。

exp <- apply(avg_df, 2, scale)rownames(exp) <- rownames(avg_df)# 組件h_state <- Heatmap(t(exp),                   column_title = "state_gsva",                   col = colorRampPalette(c('#1A5592','white',"#B83D3D"))(100),                   name= "gsva ",                   show_row_names = TRUE,                   show_column_names = TRUE)
h_state