这个R包来源于2021年发表在 Cancer Discovery 的一篇文章,
它的github地址为:https://github.com/wu-yc/scMetabolism
在某乎平台某生物公司利用这个包做起了需求,搞起了联系客服。这里我们直接带大家攻破,免费的代码用起来,不过感觉这个包没那么好用,功能并没有那么强大,大家看需使用吧~
此包限制:
细胞数不超过6万5
只支持人类 scRNA-seq 数据
0.安装此包:
install.packages(c("devtools", "data.table", "wesanderson", "Seurat", "devtools", "AUCell", "GSEABase", "GSVA", "ggplot2","rsvd"))
devtools::install_github("YosefLab/VISION")
devtools::install_github("wu-yc/scMetabolism")
1.开始使用加载包和演示数据:
load(file = "pbmc_demo.rda")
library(scMetabolism)
library(ggplot2)
library(rsvd)
演示数据下载地址:https://figshare.com/articles/dataset/scMetabolism_-_pbmc_demo_rda/13670038
2.用 Seurat 对象来做量化单细胞代谢
countexp.Seurat<-sc.metabolism.Seurat(obj = countexp.Seurat, method = "VISION", imputation = F, ncores = 2, metabolism.type = "KEGG")
obj
是一个包含 UMI 计数矩阵的 Seurat 对象。
method
支持VISION
、AUCell
、ssgsea
和gsva
,其中 VISION 是默认方法。
imputation
允许用户选择是否在新陈代谢评分之前估算他们的数据。
ncores
是并行计算的线程数。
metabolism.type
支持KEGG
和REACTOME
,其中 KEGG 包含 85 条代谢途径,而 REACTOME 包含 82 条代谢途径。
运行完以后要提取新陈代谢分数,只需运行metabolism.matrix <- countexp.Seurat@assays
计算前:
计算后:
就是每个细胞在KEGG条代谢通路的一个评分矩阵,行为通路名,列为细胞ID
3.可视化
DimPlot :指定某一通路画评分的 DimPlot
DimPlot.metabolism(obj = countexp.Seurat, pathway = "Glycolysis / Gluconeogenesis", dimention.reduction.type = "umap", dimention.reduction.run = F, size = 1)
countexp.Seurat
是一个包含 UMI 计数矩阵和代谢评分的 Seurat 对象
pathway
是要可视化的感兴趣通路
dimention.reduction.type
支持umap
和tsne
dimention.reduction.run
允许用户选择是否重新运行给定 Seurat 对象的降维。
size
是图中的点大小
该函数返回一个ggplot对象,用户可以自己DIY。
input.pathway <- rownames(countexp.Seurat@assays[["METABOLISM"]][["score"]])[1:30]
DotPlot.metabolism(obj = countexp.Seurat, pathway = input.pathway, phenotype = "ident", norm = "y")
# #Box plot
input.pathway <- rownames(countexp.Seurat@assays[["METABOLISM"]][["score"]])[1:3]
BoxPlot.metabolism(obj = countexp.Seurat, pathway = input.pathway, phenotype = "ident", ncol = 1)
本人使用都报错了,不好用不推荐
sce_Metal_exp = countexp.Seurat
sce_Metal_exp$celltype = sce_Metal_exp$ident
mscore_data = data.frame(t(sce_Metal_exp@assays[["METABOLISM"]][["score"]]),sce_Metal_exp$celltype)
avg_sM=aggregate(mscore_data[,1:ncol(mscore_data)-1],list(mscore_data$sce_Metal_exp.celltype),mean)
rownames(avg_sM) = avg_sM$Group.1
avg_sM=data.frame(t(avg_sM[,-1]))
avg_sM$KEGG = rownames(sce_Metal_exp@assays[["METABOLISM"]][["score"]])
rownames(avg_sM)=avg_sM$KEGG
#取细胞类型的前5
c_k_l = c()
for(c in c(1:ncol(avg_sM))){
c_k=avg_sM[order(avg_sM[,c]),]$KEGG[1:5]
c_k_l=c(c_k_l,c_k)
}
c_k_l= unique(c_k_l)
c_k_d = avg_sM[avg_sM$KEGG %in%c_k_l,]
#画 heatmap. 图
rownames(c_k_d) = c_k_d$KEGG
pheatmap::pheatmap(c_k_d[,-ncol(c_k_d)],show_colnames = T,scale='row')
# 获取的 没种细胞类型 代谢的评价值,保存
write.csv(avg_sM, file='ave_sMetaByCelltype.csv')
简单的使用到这里就完了,发现如果想查看代谢通路评分高的通路里是哪些基因的表达进行贡献的,有点无从下手,还是不知道,又得跑个富集通路分析KEGG才可以呢,如果你使用这个scMetablism R包,你是看上它的啥呢?