试一试这个R包做单细胞分析-scMetablism
这个R包来源于2021年发表在 Cancer Discovery 的一篇文章，
- Yingcheng Wu, Shuaixi Yang, Jiaqiang Ma, Zechuan Chen, Guohe Song, Dongning Rao, Yifei Cheng, Siyuan Huang, Yifei Liu, Shan Jiang, Jinxia Liu, Xiaowu Huang, Xiaoying Wang, Shuangjian Qiu, Jianmin Xu, Ruibin Xi, Fan Bai, Jian Zhou, Jia Fan, Xiaoming Zhang, and Qiang Gao. Spatiotemporal Immune Landscape of Colorectal Cancer Liver Metastasis at Single-Cell Level. Cancer Discovery. 2021.
它的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)
4. 直接使用表达矩阵以及网页在线版量化单细胞代谢（不推荐）
本人使用都报错了，不好用不推荐
5. 不同细胞类型的代谢通路评分均值
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包，你是看上它的啥呢？