最近微信群有小伙伴说想让出一期scMetabolism单细胞代谢，之前我也接触过这个R包，其实非常简单，而且这个R包不是很完善，也有很大的局限性，例如只能做人的，其他物种的需要同源转化。不过我们还是可以学习一下，代谢我们接触的不多，但是后期其他的代谢R包有机会还是会讲的。我们用人的单细胞数据进行演示，并进行可视化：
# devtools::install_github("YosefLab/VISION")
# devtools::install_github("YosefLab/VISION@v2.1.0")
setwd('D:/KS项目/公众号文章/scMetabolism单细胞代谢分析')
devtools::install_github("wu-yc/scMetabolism")
library(scMetabolism)
library(Seurat)
library(ggplot2)
library(rsvd)
library(pheatmap)
#官网链接：https://github.com/wu-yc/scMetabolism/tree/main
human_data <- readRDS("D:/KS项目/公众号文章/human_data.rds")
human_countexp_Seurat<-sc.metabolism.Seurat(obj = human_data,
method = "AUCell",
imputation =F,
ncores = 2,
metabolism.type = "KEGG")
#可视化1---热图====================================================================
T_metabolism <- subset(human_countexp_Seurat, celltype=='T cell')
df = data.frame(t(T_metabolism@assays[["METABOLISM"]][["score"]]))
names(T_metabolism$orig.ident)
#注意细胞名对应
rownames(df) <- gsub(".", "-", rownames(df), fixed = TRUE)
df = df[names(T_metabolism$orig.ident),]
df$orig.ident <- T_metabolism$orig.ident
avg_df =aggregate(df[,1:ncol(df)-1],list(df$orig.ident),mean)
rownames(avg_df) = avg_df$Group.1
avg_df=avg_df[,-1]
avg_df <- as.data.frame(t(avg_df))
avg_df_sec <- sample_n(avg_df, 20)
pheatmap(avg_df_sec, show_colnames = T,scale='row', cluster_rows = T,
color=colorRampPalette(c('#1A5592','white',"#B83D3D"))(100),
cluster_cols = T)
image.png
也可以用R包自带的气泡图函数可视化，其他的可视化这里就不涉及了。
DotPlot.metabolism(obj = T_metabolism, pathway = rownames(T_metabolism@assays[["METABOLISM"]][["score"]])[1:20],
phenotype = "orig.ident", norm = "y")
image.png
人的数据分析很简单，主要的问题是别的物种，需要同源转化，这里我们演示小鼠的，直接将同源转化和代谢分析包装在一个简单的函数里面，运行就可以得到分析结果了。
#同源转化参考的是之前帖子提高的一个作者写的函数：https://www.jianshu.com/p/6495706bac53
library(Seurat)
library(nichenetr)
library(dplyr)
Mouse.sc.metabolism <- function(obj,
metabolism.type=c("KEGG","REACTOME")){
#将鼠的基因名转化为人的
gene_trans = rownames(obj) %>% convert_mouse_to_human_symbols() %>% as.data.frame()
gene_mouse <- as.data.frame(rownames(obj))
gene_use <- cbind(gene_trans, gene_mouse)
gene_use <- na.omit(gene_use)
colnames(gene_use) <- c('human','mouse')
mouse_data_trans <- subset(obj,features=gene_use$mouse)
#函数参考
#https://www.jianshu.com/p/6495706bac53
RenameGenesSeurat <- function(obj,newnames,gene.use=NULL,de.assay) {
# Replace gene names in different slots of a Seurat object. Run this before integration. Run this before integration.
# It only changes obj@assays$RNA@counts, @data and @scale.data.
print("Run this before integration. It only changes obj@assays$*@counts, @data and @scale.data, @var.features,@reductions$pca@feature.loadings")
lassays <- Assays(obj)
#names(obj@assays)
assay.use <- obj@reductions$pca@assay.used
DefaultAssay(obj) <- de.assay
if (is.null(gene.use)) {
all_genenames <- rownames(obj)
}else{
all_genenames <- gene.use
obj <- subset(obj,features=gene.use)
}
order_name <- function(v1,v2,ref){
v2 <- make.names(v2,unique=T)
df1 <- data.frame(v1,v2)
rownames(df1) <- df1$v1
df1 <- df1[ref,]
return(df1)
}
df1 <- order_name(v1=all_genenames,v2=newnames,ref=rownames(obj))
all_genenames <- df1$v1
newnames <- df1$v2
if ('SCT' %in% lassays) {
if ('SCTModel.list' %in% slotNames(obj@assays$SCT)) {
obj@assays$SCT@SCTModel.list$model1@feature.attributes <- obj@assays$SCT@SCTModel.list$model1@feature.attributes[all_genenames,]
rownames(obj@assays$SCT@SCTModel.list$model1@feature.attributes) <- newnames
}
}
change_assay <- function(a1=de.assay,obj,newnames=NULL,all_genenames=NULL){
RNA <- obj@assays[a1][[1]]
if (nrow(RNA) == length(newnames)) {
if (length(RNA@counts)) RNA@counts@Dimnames[[1]] <- newnames
if (length(RNA@data)) RNA@data@Dimnames[[1]] <- newnames
if (length(RNA@var.features)) {
df1 <- order_name(v1=all_genenames,v2=newnames,ref=RNA@var.features)
all_genenames1 <- df1$v1
newnames1 <- df1$v2
RNA@var.features <- newnames1
}
if (length(RNA@scale.data)){
df1 <- order_name(v1=all_genenames,v2=newnames,ref=rownames(RNA@scale.data))
all_genenames1 <- df1$v1
newnames1 <- df1$v2
rownames(RNA@scale.data) <- newnames1
}
} else {"Unequal gene sets: nrow(RNA) != nrow(newnames)"}
obj@assays[a1][[1]] <- RNA
return(obj)
}
for (a in lassays) {
DefaultAssay(obj) <- a
df1 <- order_name(v1=all_genenames,v2=newnames,ref=rownames(obj))
all_genenames1 <- df1$v1
newnames1 <- df1$v2
obj <- change_assay(obj=obj,a1=a,newnames=newnames1,all_genenames=all_genenames1)
}
hvg <- VariableFeatures(obj,assay=assay.use)
if (length(obj@reductions$pca)){
df1 <- order_name(v1=all_genenames,v2=newnames,ref=hvg)
df1 <- df1[rownames(obj@reductions$pca@feature.loadings),]
all_genenames1 <- df1$v1
newnames1 <- df1$v2
rownames(obj@reductions$pca@feature.loadings) <- newnames1
}
try(obj[[de.assay]]@meta.features <- data.frame(row.names = rownames(obj[[de.assay]])))
return(obj)
}
#转化
mouse_data_trans <- RenameGenesSeurat(mouse_data_trans,
newnames = gene_use$human,
gene.use = gene_use$mouse,
de.assay = 'RNA')
mouse_metabolism <- sc.metabolism.Seurat(obj = mouse_data_trans,
method = "AUCell",
imputation =F,
ncores = 2,
metabolism.type = metabolism.type)
return(mouse_metabolism)
}
mouse_results <- Mouse.sc.metabolism(mouse_data, metabolism.type = 'KEGG')
之后做小鼠的，只需要source这个函数，依据代码就完成了。可视化和前面一样。这个包其实很简单也没有什么好说的了，可能有些小伙伴的mouse数据做的时候还是会出错，需要自己去修改函数了哦。希望这个分享对你有帮助，觉得有用的点个赞、分享一下呗！
网友评论