scMetabolism：一个简单的单细胞代谢分析R包

最近微信群有小伙伴说想让出一期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数据做的时候还是会出错，需要自己去修改函数了哦。希望这个分享对你有帮助，觉得有用的点个赞、分享一下呗！