我们在第一期教程中：METAflux | 实操教程1-简介和快速上手 ，简单介绍了Bulk RNA-seq和scRNA_seq的大致流程，没有对其进行详细解析。本期推文小编将详细介绍Bulk RNA-seq的工作流程。
1.Load the library
library(METAFlux)
2.Load data
2.1 Load gene expresssion data
METAFlux 需要基因表达数据作为输入。
基因表达矩阵应为基因 x 样本矩阵，其中行名称为人类基因名称（基因符号），列名称应为样本名称。请注意，METAFlux 不支持其他基因 ID。 在使用 METAFlux 之前，应将输入基因表达矩阵标准化（例如，对数转换等）。METAflux 不会对表达式数据执行任何归一化。 基因表达数据不能为负值。
2.2 Load the METAFlux underlying GEM information
在开始之前，我们可以看一下Human-GEM 文件。我们使用 Human-GEM（由 13082 个代谢反应和 8378 个代谢物组成）作为我们的基础代谢模型。对于每个反应，我们都有一个唯一的反应 ID 和 SUBSYSTEM（即该反应属于哪个途径）。
library(METAFlux)data("human_gem")
对于每个反应，我们还从该文件中获取其他重要信息，例如：反应是否可逆，发生在哪个隔室中，以及涉及哪些代谢物和基因。1对1关系在下面的卡通中用颜色编码。
2.3 Load the METAFlux medium
加载代谢物数据
data("cell_medium")# for cell line model
data("human_blood")# for human derived samples
2.3.1 what's is "medium" and why do we need it?
在METAFlux流程过程中，"medium"代表的是在肿瘤微环境中（TME）的营养物质（代谢产物）,这些营养物质可以被细胞吸收
option 1. if users do not have prior knowledge about their medium composition
如果对代谢物没有先验知识，这里提供2种通用代谢物。METAFlux中通用有 2 个 medium代谢物文件：cell_medium 和 human_blood medium。Cell_medium 可用于细胞系模型，人血培养基可用于患者来源的样品。细胞系培养基含有 44 种代谢物，其相应的交换反应可在“reaction_name”列下找到。人血培养基在人血中含有 64 种代谢物。
option 2. if users have prior knowledge about their medium composition
要删除营养物质，只需删除该营养物质所在的整行，但请不要更改列名称或删除列；要添加额外的营养物质，请参阅营养物质查找文件。
data("nutrient_lookup_files")
“营养查找文件”包含 1648 个交换反应（注意：交换反应是人工反应！将它们包括在这里，因为这是细胞如何摄取或释放代谢物的数学表示）。例如，如果我们需要添加代谢物 “萘”，可以在搜索框中搜索，然后就会找到相应的营养反应和方程式。接下来，需要将新的萘和HMR_7110行添加到当前的 medium 文件中。此更新的 medium 文件稍后将用于建模。
3.Calculate MARS (Metabolic Reaction Activity Score）
我们可以用 GPR（基因-蛋白质-反应）从基因表达中计算单个样品归一化后的代谢反应活性MRAS。
scores<-calculate_reaction_score(bulk_test_example)
head(scores)
4.Calculate flux
我们计算13082个反应的代谢通量。
#flux<-compute_flux(mras=scores,medium = cell_medium)#if data are cell line samples
flux<-compute_flux(mras=scores,medium = human_blood)#if data are human derived samples
5.Inspecting and interpreting the flux data
通量的迹象具有生物学意义。这里 通量 的 “sign” 表示方向。例如，对于营养物质吸收/释放情况（营养物质查找文件中的 1648 个交换反应），正值表示化合物的释放，负值表示化合物的吸收。对于其他反应，正值表示净通量是向前的，负的通量表示净通量是反向的。绝对值表示磁通量的大小。 需要注意的一点是，由于我们正在最小化模型中总通量的总和，因此我们将得到一个简洁的通量数据输出，这意味着许多反应将接近 0 通量。（例如，有些反应应该只向前进行，但预测的通量有一个非常小的负号，可以认为接近 0 通量） 提取目标数据。如果对代谢物的摄取或释放非常感兴趣，可以搜索 “nutrient lookup file” 来获得代谢物交换反应 ID。例如，如果我们想知道葡萄糖摄取，我们需要搜索 'glucose' 来获得葡萄糖摄取反应，即 HMR_9034。这些可以被视为葡萄糖代谢物摄取率。接下来，我们可以提取数据（提取感兴趣的代谢物的摄取率）：
data("nutrient_lookup_files")
glucose<-data.frame(glucose=flux[grep("HMR_9034",human_gem$ID),])
glucose
## glucose
## Sample1 -0.1576912
## Sample2 -0.1896533
## Sample3 -0.1522470
## Sample4 -0.1692082
## Sample5 -0.1902032
## 这里就看到了具体样本的具体代谢物活性
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),])
我们还可以计算142条通路的通路水平活性。对于给定的通路，活性水平可以计算为与该通路相关的反应的平均通量，通路的活性水平是由相关代谢反应的通路反映。
#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
mapal <- colorRampPalette(RColorBrewer::brewer.pal(11,"RdBu"))(256)
pheatmap::pheatmap(all_pathway_score,cluster_cols = F,color = rev(mapal),scale = "row")
如图所示展示就是研究对象的142条代谢通路活性的热图。
本期分享到这里就结束了，我们下期再会~~
