单细胞METAFlux分析学习和整理

既往的代谢分析是基于通量平衡分析(Flux balance analysis, FBA)，旨在通过稳态条件下的代谢通量优化来预测细胞在不同环境或遗传背景下的代谢行为。通量平衡分析假设代谢网络处于稳态，通过线性规划优化目标函数（例如生长率或产物生成率），以模拟细胞在给定条件下的代谢通量分布。

METAFlux是基于METAbolic Flux balance analysis (通量平衡分析的框架之上)，以营养感知的方式表征整个代谢回路，并输出非退化性的通量，一共有13082个代谢物可以计算。

步骤流程

1.导入

rm(list = ls())
library(stringr)
library(Seurat)  开发者建议使用SeuratV4版本数据
library(METAFlux)
library(scater)
library(scran)
library(dplyr)
library(patchwork)
library(qs)
library(BiocParallel)
register(MulticoreParam(workers = 8, progressbar = TRUE))

sc_dataset <- qread("sc_dataset.qs") # 单细胞
load("consensus.Rdata") # bulk数据(TCGA-HNSC)
#load("scRNA.Rdata")

2.bulk数据分析

# bulkRNA数据
# 行是基因名(必须是gene symbol)，列是样本名
# 表达矩阵必须是Normalized之后的比如log转换
# 表达值不能为负值
# Calculate Metabolic Reaction Activity Score(mras) for cell line data
scores <- calculate_reaction_score(exprSet)
# Calculate flux for cell line data
# 设置medium的原因是其反应了TME中的营养物质情况，其可被细胞所吸收
# 如果没有medium的先验知识就根据数据是细胞系还是人类RNAseq数据来选择介质文件：cell_medium/human_blood
# 如果了解medium的成为，就需要修改内置文件中的信息，可见github流程(开发者提供了Nutrient look up files)
flux <- compute_flux(mras=scores,
                     medium = human_blood) 

# 可选: flux scores cubic root normalization
# cbrt <- function(x) {
#     sign(x) * abs(x)^(1/3)
# }
# 
# flux=cbrt(flux)

3.可视化

# 柱状图
# 提取特定代谢物的ID
# 代谢通量的sign代表方向，正值表达化合物的释放，负值表示化合物的吸收
# 数值大小可认为是释放/吸收的强度
data("nutrient_lookup_files")
glucose <- data.frame(glucose = flux[grep("HMR_9034",human_gem$ID),])
head(glucose)
#                       glucose
# TCGA-CR-7374-01A -0.004602662
# TCGA-CV-A45V-01A -0.004763098
# TCGA-CV-7102-01A -0.004122545
# TCGA-MT-A67D-01A -0.004472756
# TCGA-P3-A6T4-01A -0.005428906
# TCGA-CV-7255-01A -0.003568878

library(ggplot2)
ggplot(glucose,aes(y=-glucose,x="sample"))+
  geom_boxplot()+
  ggtitle("Glucose uptake level")+
  xlab("")+
  ylab("Glucose uptake scores")

# 其他反应
# HMR_4363: 2-phospho-D-glycerate[c] <=> H2O[c] + PEP[c]
HMR_4363 <- data.frame(hmr4363=flux[grep("HMR_4363",human_gem$ID),])

# 热图可视化
#compute pathway level activity for all samples
pathway <- unique(unlist(human_gem$SUBSYSTEM))
pathway_score<-list()
for (i in pathway){
path=i
activity_score<-c()
for (d in 1:ncol(flux)){
activity_score[d]<-mean(abs(flux[which(unlist(human_gem$SUBSYSTEM)==i),d]))
} 
pathway_score[[i]]<-activity_score
}

all_pathway_score <- as.data.frame(do.call(rbind,pathway_score))

#heatmap 
pdf("metabolic_res.pdf",width = 20,height = 30)
mapal <- colorRampPalette(RColorBrewer::brewer.pal(11,"RdBu"))(256)
pheatmap::pheatmap(all_pathway_score,
                   cluster_cols = F,
                   color = rev(mapal),
                   scale = "row")
dev.off()

可以提取特定代谢物进行分析，这里可以提取数据进行样本分组就可以对比干预与否下代谢物的摄取情况也可以对多个代谢物和样本进行热图绘制

4.单细胞数据分析-单样本

# 计算cluster/细胞类型水平上对通量建模，而非cell水平
# sample code for 4 clusters in multiple patients
# 纳入分析的数据需要进行标准化处理
obj.list <- SplitObject(sc_input,
                        split.by = "orig.ident")
sc <- obj.list[[1]]
mean_exp = calculate_avg_exp(myseurat = sc,
                             myident = 'celltype',
                             n_bootstrap= 1, # 自行分析的时候把该值写大来
                             seed=1)
scores <- calculate_reaction_score(data = mean_exp)

# 计算细胞类型/样本的分数
#g=round(table(sc$Cell_type)/nrow(sc@meta.data),3)
g=table(sc$celltype)/nrow(sc@meta.data)
print(g)
#          B/plasma cells         dendritic cells       endothelial cells 
#              0.01680672              0.02941176              0.49579832 
# epithelial/cancer cells             fibroblasts              mast cells 
#              0.11764706              0.17226891              0.02521008 
#           myeloid cells              T/NK cells 
#              0.07563025              0.06722689 
flux=compute_sc_flux(num_cell = 8,
                     fraction =c(g[1],g[2],g[3],g[4],
                                 g[5],g[6],g[7],g[8]),
                     fluxscore=scores,
                     medium = human_blood)
qsave(flux,paste0("object",i,".qs"))

5.可视化

library(ggplot2)
glucose$celltype=rownames(glucose)
long_glucose=reshape2::melt(glucose,id.vars='celltype')
ggplot(long_glucose,aes(y=-value,x=celltype))+
  geom_boxplot()+
  ggtitle("Glucose uptake level for different cell types")+
  xlab("")+
  ylab("Glucose uptake scores")+ 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
ggsave("total.pdf",width = 30,height = 4)

粗看之下图片很奇怪，明明是柱状图为什么没有出现柱子的情况~ 这是因为n_bootstrap只设置了1，自行分析的时候需要重新设置并且把X轴的名称更改一下。

6.单细胞数据分析-多样本

这个函数运行的非常慢，本来想多样本试一试，但是发现实在是太慢了！但是流程还是差不多的，只要拆解之后就能得到不同样本下的不同细胞系的代谢物结果。

# 多样本的单细胞数据
# 计算cluster/细胞类型水平上对通量建模，而非cell水平
# sample code for 4 clusters in multiple patients
# 纳入分析的数据需要进行标准化处理
DefaultAssay(sc_dataset) <- "RNA"
obj.list <- SplitObject(sc_dataset,
                        split.by = "orig.ident")
for (i in c(1:length(obj.list))){
sc <- obj.list[[i]]
mean_exp = calculate_avg_exp(myseurat = sc,
                             myident = 'celltype',
                             n_bootstrap= 1, # 自行分析的时候把该值写大来
                             seed=1)
scores <- calculate_reaction_score(data = mean_exp)

# 计算细胞类型/样本的分数
#g=round(table(sc$Cell_type)/nrow(sc@meta.data),3)
g=table(sc$celltype)/nrow(sc@meta.data)
print(g)
flux=compute_sc_flux(num_cell = 8,
                     fraction =c(g[1],g[2],g[3],g[4],
                                 g[5],g[6],g[7],g[8]),
                     fluxscore=scores,
                     medium = human_blood)
qsave(flux,paste0("object",i,".qs"))
}

总之这个工具可以作为代谢分析中的辅助分析手段，帮助研究者探索或者验证自己的结果。

参考资料：

1. Characterizing cancer metabolism from bulk and single-cell RNA-seq data using METAFlux. Nat Commun. 2023 Aug 12;14(1):4883.
2. METAFlux: https://github.com/KChen-lab/METAFlux
3. GenePulse: https://mp.weixin.qq.com/s/E49dgvDI75dpw2wKxrdVuQ
4. 生信宝库：https://mp.weixin.qq.com/s/EwGDgGiEB47gH7EVDraZrQ

注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟