作者，Evil Genius

最近發(fā)生了很多事，電信詐騙40萬，心態(tài)沒有崩潰；被網(wǎng)爆，被罵，心態(tài)也沒有崩潰；失業(yè)，心態(tài)也沒有崩潰；回家了看著父母的狀態(tài)，心態(tài)也沒有崩潰；但是沒有了服務(wù)器分析數(shù)據(jù)，心態(tài)真的撐不住了，特別焦慮。

答應(yīng)過大家要做全套視頻講解單細胞空間的分析內(nèi)容，也因為很多人罵我而暫時擱置了。沒有辦法，人言可畏，當然，自己也有錯。

好了，開始我們今天的內(nèi)容，單細胞代謝推斷，我們從最新的三個軟件，compass，scmetabolism以及scFEA入手分析。

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

細胞代謝可以調(diào)節(jié)正常細胞功能以及多種疾病的病理生理學。最近，免疫細胞代謝成為一個非常受關(guān)注的研究領(lǐng)域，多項研究探索了抗病毒免疫、自身免疫和抗腫瘤反應(yīng)的炎癥代謝調(diào)節(jié)。由于代謝網(wǎng)絡(luò)的復(fù)雜性，某代謝的擾動可以產(chǎn)生級聯(lián)效應(yīng)，從而影響網(wǎng)絡(luò)中看似相對較遠的部分，因此亟待工具將特定反應(yīng)/酶的觀測與代謝及其在疾病中的失調(diào)系統(tǒng)地聯(lián)系起來。

Compass算法展示

已開發(fā)并應(yīng)用在大塊組織數(shù)據(jù)的通量平衡分析(Flux Balance Analysis,FBA)是一個成功的方法，而scRNA-seq的出現(xiàn)為研究單個細胞的代謝異質(zhì)性提供了基礎(chǔ)。Allon Wagner等人在FBA的基礎(chǔ)上構(gòu)建了Compass算法，該算法將scRNA-seq數(shù)據(jù)與代謝網(wǎng)絡(luò)相結(jié)合從而推斷單個細胞的代謝狀態(tài)。

圖1 Compass算法概述

分析實例

使用說明

簡單運行Compass可以只設(shè)置3個參數(shù)，包含表達矩陣(--data)、進程個數(shù)(--num-processes)、物種(--species，當前支持的物種有人和小鼠。設(shè)置物種后，計算中會自動使用與對應(yīng)物種適配的基因名)，其余參數(shù)可以先使用默認值，后續(xù)可以根據(jù)研究需要進行修改。需要注意的是，Compass計算時間相對較長，每個細胞大約需要30分鐘，而另一個用于單細胞代謝預(yù)測的工具——scFEA計算速度相對快得多，例如使用8個線程預(yù)測100個細胞的測試數(shù)據(jù)大約只需要1分鐘。因此針對較大的數(shù)據(jù)集，Compass可以通過(--microcluster-size)將細胞劃分為cluster，再以cluster的平均值表征該cluster，這里也可以使用其他類似方法，例如MetaCell。此外，根據(jù)實際研究問題，需要考慮是否將不同表型的細胞分別進行microcluster-size。

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

運行完上面的命令后，會在當前目錄生成reaction.csv文件。其中的值是反應(yīng)罰分（reaction penalties），高得分表示該反應(yīng)的可能性較低。

上述獲得的是每個反應(yīng)在每個樣品中的反應(yīng)罰分，我們需要根據(jù)Recon2 meta信息（圖5）將A列的反應(yīng)編碼轉(zhuǎn)成人們可讀的代謝通路名字，將反應(yīng)罰分進行轉(zhuǎn)化（轉(zhuǎn)成數(shù)值大的表示反應(yīng)活性高，加1，取-log）,然后使用wilcoxon或者其他方法進行差異分析。

目前來看Compass功能是最強大也相對最準確的，其相應(yīng)的文章于2021年發(fā)表在Cell雜志上，但是對人的要求相對較高

最“水”的代謝單細胞分析軟件scMetablism

scMetabolism是復(fù)旦大學中山醫(yī)院研發(fā)的一款基于單細胞水平量化代謝活性的軟件，該軟件基于常規(guī)單細胞矩陣文件，采用VISION算法對每個細胞進行評分，最終得到細胞在每一條代謝通路中的活性得分。

scMetabolism分析

軟件預(yù)置了人的代謝基因集，包括85條KEGG通路和82條REACTOME條目。若物種為大、小鼠等哺乳動物，可以通過同源轉(zhuǎn)換后進行分析，其他物種也可選擇合適的基因集進行代謝活性分析。

分析使用到的文件包括單細胞轉(zhuǎn)錄組數(shù)據(jù)和對應(yīng)的代謝基因集，可以使用單細胞轉(zhuǎn)錄組的RDS文件和軟件自帶的代謝基因集。

代謝活性分析

首先，使用Vision算法進行代謝打分計算得分。結(jié)果為打分均值統(tǒng)計表，average_KEGG_score表格的第一列為匹配到代謝基因集的代謝條目名稱，每一列為不同分組的得分均值。

發(fā)表于Genome Research的代謝分析軟件---scFEA

主要內(nèi)容

細胞間的代謝異質(zhì)性，和代謝物之間的相互作用被認為是疾病治療耐藥性的重要影響因素。然而，由于缺乏成熟的高通量單細胞代謝組學技術(shù)，我們尚未對組織內(nèi)代謝異質(zhì)性和協(xié)同機制建立系統(tǒng)的認識。為了緩解這一技術(shù)短板，作者開發(fā)了一種新的計算方法，即單細胞通量估計分析(scFEA)，以從scRNA-seq數(shù)據(jù)推斷細胞內(nèi)代謝物的通量研究。

scFEA分析代謝主要是基于兩個假設(shè)：（1）代謝模塊的通量變化可以建模為催化酶轉(zhuǎn)錄組水平變化的非線性函數(shù)；(2) 所有中間底物的總通量不平衡應(yīng)在所有單細胞中最小化。

scFEA同樣是基于已有的一些代謝相關(guān)的功能通路，不過作者將代謝圖譜中的網(wǎng)絡(luò)圖案被簡化為代謝模塊，其中反應(yīng)和代謝物用黑色和藍色的矩形表示，模塊和代謝物用綠色和粉色表示。

%%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預(yù)測的細胞通量組可用于一系列下游分析，包括識別具有共同代謝變化的代謝模塊或細胞群，對酶對整個代謝通量影響的敏感性評估，以及對細胞組織和細胞間代謝通信的推斷。

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