我们在第一期教程中: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条代谢通路活性的热图。
本期分享到这里就结束了,我们下期再会~~