10X单细胞（10X空间转录组）数据分析之代谢推断篇

作者，Evil Genius

最近发生了很多事，电信诈骗40万，心态没有崩溃；被网爆，被骂，心态也没有崩溃；失业，心态也没有崩溃；回家了看着父母的状态，心态也没有崩溃；但是没有了服务器分析数据，心态真的撑不住了，特别焦虑。

答应过大家要做全套视频讲解单细胞空间的分析内容，也因为很多人骂我而暂时搁置了。没有办法，人言可畏，当然，自己也有错。

好了，开始我们今天的内容，单细胞代谢推断，我们从最新的三个软件，compass，scmetabolism以及scFEA入手分析。

Compass---Cell文章引用的代谢推断方法（Metabolic modeling of single Th17 cells reveals regulators of autoimmunity）

细胞代谢可以调节正常细胞功能以及多种疾病的病理生理学。最近，免疫细胞代谢成为一个非常受关注的研究领域，多项研究探索了抗病毒免疫、自身免疫和抗肿瘤反应的炎症代谢调节。由于代谢网络的复杂性，某代谢的扰动可以产生级联效应，从而影响网络中看似相对较远的部分，因此亟待工具将特定反应/酶的观测与代谢及其在疾病中的失调系统地联系起来。

Compass算法展示

已开发并应用在大块组织数据的通量平衡分析(Flux Balance Analysis,FBA)是一个成功的方法，而scRNA-seq的出现为研究单个细胞的代谢异质性提供了基础。Allon Wagner等人在FBA的基础上构建了Compass算法，该算法将scRNA-seq数据与代谢网络相结合从而推断单个细胞的代谢状态。

图1 Compass算法概述

分析实例

如上图所示，使用单细胞转录组（Th17p细胞和Th17n细胞）和代谢模型作为COMPASS的输入，经过流平衡分析（compass代谢分析算法），找到了两群细胞间异常的代谢通路及相关的反应/酶，并进行了分子和功能验证，以及小鼠体内验证。

使用说明

简单运行Compass可以只设置3个参数，包含表达矩阵(--data)、进程个数(--num-processes)、物种(--species，当前支持的物种有人和小鼠。设置物种后，计算中会自动使用与对应物种适配的基因名)，其余参数可以先使用默认值，后续可以根据研究需要进行修改。需要注意的是，Compass计算时间相对较长，每个细胞大约需要30分钟，而另一个用于单细胞代谢预测的工具——scFEA计算速度相对快得多，例如使用8个线程预测100个细胞的测试数据大约只需要1分钟。因此针对较大的数据集，Compass可以通过(--microcluster-size)将细胞划分为cluster，再以cluster的平均值表征该cluster，这里也可以使用其他类似方法，例如MetaCell。此外，根据实际研究问题，需要考虑是否将不同表型的细胞分别进行microcluster-size。

compass --data-mtx expression.mtx genes.tsv sample_names.tsv --num-processes 10 --lambda 0.25 --species homo_sapiens 

运行完上面的命令后，会在当前目录生成reaction.csv文件。其中的值是反应罚分（reaction penalties），高得分表示该反应的可能性较低。

上述获得的是每个反应在每个样品中的反应罚分，我们需要根据Recon2 meta信息（图5）将A列的反应编码转成人们可读的代谢通路名字，将反应罚分进行转化（转成数值大的表示反应活性高，加1，取-log）,然后使用wilcoxon或者其他方法进行差异分析。

目前来看Compass功能是最强大也相对最准确的，其相应的文章于2021年发表在Cell杂志上，但是对人的要求相对较高

最“水”的代谢单细胞分析软件scMetablism

scMetabolism是复旦大学中山医院研发的一款基于单细胞水平量化代谢活性的软件，该软件基于常规单细胞矩阵文件，采用VISION算法对每个细胞进行评分，最终得到细胞在每一条代谢通路中的活性得分。

scMetabolism分析

软件预置了人的代谢基因集，包括85条KEGG通路和82条REACTOME条目。若物种为大、小鼠等哺乳动物，可以通过同源转换后进行分析，其他物种也可选择合适的基因集进行代谢活性分析。

分析使用到的文件包括单细胞转录组数据和对应的代谢基因集，可以使用单细胞转录组的RDS文件和软件自带的代谢基因集。

代谢活性分析

首先，使用Vision算法进行代谢打分计算得分。结果为打分均值统计表，average_KEGG_score表格的第一列为匹配到代谢基因集的代谢条目名称，每一列为不同分组的得分均值。

类似于平时做的基于已有的功能数据集GSVA或者说ssGESA分析，其分析出的代谢结果可参考的价值并不大。

发表于Genome Research的代谢分析软件---scFEA

主要内容

细胞间的代谢异质性，和代谢物之间的相互作用被认为是疾病治疗耐药性的重要影响因素。然而，由于缺乏成熟的高通量单细胞代谢组学技术，我们尚未对组织内代谢异质性和协同机制建立系统的认识。为了缓解这一技术短板，作者开发了一种新的计算方法，即单细胞通量估计分析(scFEA)，以从scRNA-seq数据推断细胞内代谢物的通量研究。

scFEA分析代谢主要是基于两个假设：（1）代谢模块的通量变化可以建模为催化酶转录组水平变化的非线性函数；(2) 所有中间底物的总通量不平衡应在所有单细胞中最小化。

scFEA同样是基于已有的一些代谢相关的功能通路，不过作者将代谢图谱中的网络图案被简化为代谢模块，其中反应和代谢物用黑色和蓝色的矩形表示，模块和代谢物用绿色和粉色表示。

%%bash
cd /Users/chang/Documents/work/flux/scFEA
python src/scFEA.py --data_dir data --input_dir input /
                    --test_file Melissa_full.csv /
                    --moduleGene_file module_gene_m168.csv /
                    --stoichiometry_matrix cmMat_c70_m168.csv

from IPython.display import Image
Image(filename='./scFEA/output/loss_20210923-112543.png')

%%R
library(tidyverse)
library(rstatix)
library(ggpubr)
library(reshape)
library(ggridges)
library(ggplot2)

# load predicted flux
data_c <- read.csv('./scFEA/output/mouse_flux.csv')
data_c0 <- as.matrix(data_c[,-1])
rownames(data_c0) <- as.character(data_c[,1])
data_c0 <- t(data_c0)
ppp_all <-c()

# load cell label of mouse example data
load('./scFEA/input/mouse_example_cell_ident.RData')
yyy <- cell_id[colnames(data_c0)]
for(ii in 1:nrow(data_c0)){
    xxx <- data_c0[ii,]
    final_df <- cbind(paste('X', 1:length(xxx), sep=''), xxx, yyy)
    final_df <- as.data.frame(final_df)
    final_df[,2] <- as.numeric(final_df[,2])
    colnames(final_df) <- c('var', 'flux', 'cellType')
    pp <- sd(final_df$flux)/abs(mean(final_df$flux))
    ppp_all <- c(ppp_all, pp)
}
tg_ids <- which(ppp_all > 1e-10)

# load mouse module info
load('./scFEA/input/mouse_module_info.RData')
if(length(tg_ids) > 0){
    jj = 1 # check module 2
    xxx <- data_c0[tg_ids[jj], ]
    #final_df <- cbind(paste('X', 1:length(xxx), sep=''), xxx, yyy)
    #final_df <- as.data.frame(final_df)
    #final_df[,2] <- as.numeric(final_df[,2])
    #colnames(final_df) <- c('var', 'flux', 'cellType')
    final_df <- data.frame(var = paste('X', 1:length(xxx), sep=''),
                          flux = xxx,
                          cellType = yyy)
        
    title <- mouse_module_info[rownames(data_c0)[tg_ids[jj]], 'M_name']
    aa <- ggplot(final_df, aes(x = flux, y = cellType, fill = cellType)) +
                  geom_density_ridges() +
                  theme_ridges() +
                  theme(legend.position = 'none') +
                  ggtitle(title) +
                  theme(plot.title = element_text(hjust = 0.5))
    plot(aa)
    
}

总的来说，scFEA预测的细胞通量组可用于一系列下游分析，包括识别具有共同代谢变化的代谢模块或细胞群，对酶对整个代谢通量影响的敏感性评估，以及对细胞组织和细胞间代谢通信的推断。

好了，今天就到这里了，生活很好，有你更好