今天我们继续学习METAflux这个工具，在前面的推文中我们一起详细的学习了METAflux的文献和代码实操部分，感兴趣的小伙伴可以先阅读这部分内容：
计算反应活性
scores<-calculate_reaction_score(bulk_test_example)
相关源码及其解读
#' Calculate metabolic reaction scores (MRAS) for 13082 reactions
#'
#' @param data gene expression data.1.The gene expression matrix should be gene by sample matrix where row names are human gene names (gene symbols),
#' and column names should be sample names. Please note that METAFlux does not support other gene IDs.
#' 2.The input gene expression matrix should be normalized (e.g., log-transformed, etc.) before using METAFlux.
#' METAflux will not perform any normalization on expression data.
#' 3.Gene expression data cannot have negative values.
#'
#' @return
#' @export
#'
#' @examples calculate_reaction_score(bulk_test_example)
calculate_reaction_score <- function(data) {
if (sum(data < 0) > 0)
stop("Expression data needs to be all positive")
# make sure features are present
features <- rownames(data)
if (sum(features %in% rownames(gene_num)) == 0)
stop("Requested gene names cannot be found. Rownames of input data should be human gene names.Please check the rownames of input data.") #change
message(paste0(round(sum(features %in% rownames(gene_num))/3625 * 100,
3), "% metabolic related genes were found......"))
#计算匹配的基因名占总基因名的百分比
gene_num <- METAFlux:::gene_num#获取 gene_num 数据
Hgem <- METAFlux:::Hgem#代谢反应的信息
iso <- METAFlux:::iso#同位素标记代谢物的信息
multi_comp <- METAFlux:::multi_comp#多组分复合体的信息
simple_comp <- METAFlux:::simple_comp#简单复合体的信息
message("Computing metabolic reaction activity scores......")
core <- do.call(rbind, lapply(1:length(iso), calculate_iso_score, data = data,
list = iso, gene_num = gene_num))
core2 <- do.call(rbind, lapply(1:length(simple_comp), calculate_simple_comeplex_score,
data, list = simple_comp, gene_num = gene_num))
core3 <- do.call(rbind, lapply(1:length(multi_comp), calculate_multi_comp,
data = data, gene_num = gene_num))
message("Preparing for score matrix......")
rownames(core) <- names(iso)
rownames(core2) <- names(simple_comp)
rownames(core3) <- names(multi_comp)
big_score_matrix <- rbind(core, core2, core3)
big_score_matrix <- apply(big_score_matrix, 2, stdize, na.rm = T)
# handle NAN
big_score_matrix[is.na(big_score_matrix)] <- 0
empty_helper <- as.data.frame(Hgem$Reaction)
colnames(empty_helper) <- "reaction"
Final_df <- merge(empty_helper, big_score_matrix, all.x = T, by.x = 1,
by.y = 0)
Final_df[is.na(Final_df)] <- 1
rownames(Final_df) <- Final_df$reaction
Final_df$reaction <- NULL
Final_df <- Final_df[Hgem$Reaction, , drop = F]
#以上是对计算结果处理，首先合并成一个大的得分矩阵；、
#将得分矩阵中的NaN值替换为0；将得分矩阵和反应信息合并；
#将合并后矩阵中的NaN值替换为1；删除不再需要的反应
if (all.equal(rownames(Final_df), Hgem$Reaction)) {
message("Metabolic reaction activity scores successfully calculated \n")
} else {
message("Calculation not reliable Check input data format \n")
}
Final_df[which(Hgem$LB == 0 & Hgem$UB == 0), ] <- 0 #this is biomass reaction, do not use
return(Final_df)
}
#确保最终得分矩阵的行名与 Hgem 数据集中的反应名称匹配
#检查最终得分矩阵的行名是否与 Hgem 数据集中的反应名称完全匹配，并给出相应的信息
#将生物量反应（下界和上界都为0的反应）的得分设置为0，得到最终的代谢反应活性得分矩阵
计算代谢反应通量
#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
相关源码及其解读
#' Final optimization step for flux calculation
#'
#' @param mras metabolic reaction activity scores
#' @param medium input medium file which indicates the nutrients available in the medium.
#' We provide 2 general mediums if you have no prior knowledge about your medium: cell line medium and human blood medium if prior knowledge of medium is not available.
#' Please see tutorial for more details.
#'
#' @return Calculated fluxes
#' @export
#' @import osqp
#'
#' @examples
compute_flux <- function(mras, medium) {
message("Setting up for optimization.....")
Hgem <- METAFlux:::Hgem
P <- diag(1, ncol(Hgem$S), ncol(Hgem$S))
#创建一个与代谢网络模型的 S 矩阵列数相同大小的单位对角矩阵 P
q <- rep(0, ncol(Hgem$S))#创建一个全零向量 q，其长度与 S 矩阵的列数相同
q[which(Hgem$Obj == 1)] <- -10000
#在目标函数中，对于目标反应设置一个很大的负值，以确保在优化过程中这些反应的通量最小化
A <- as.matrix(rbind(Hgem$S, P))#将 S 矩阵和单位矩阵 P 合并，形成线性规划问题的 A 矩阵
resu <- list()#初始化两个列表，用于存储每次优化的结果和通量向量
flux_vector <- list()
message("Computing bulk RNA-seq flux.....")
Seq <- seq(1, ncol(mras))
pb <- txtProgressBar(0, length(Seq), style = 3)
#循环遍历 mras 矩阵的每一列，对每个代谢反应活性评分进行优化计算
for (i in Seq) {
setTxtProgressBar(pb, i)
origlb <- Hgem$LB#获取原始的下界
origlb[Hgem$rev == 1] <- -mras[, i][Hgem$rev == 1]#对于可逆反应，设置下界为 mras 的负值
origlb[Hgem$rev == 0] <- 0#对于不可逆反应，设置下界为 0
origlb <- origlb[, 1]#确保下界是一维向量
origub <- mras[, i]#设置上界（upper bound）为当前迭代的 mras 值
origlb[Hgem$Reaction %in% Hgem$Reaction[which(Hgem$pathway == "Exchange/demand reactions")]] <- 0
#对于交换/需求反应，设置下界为 0
origlb[Hgem$Reaction %in% medium$reaction_name] <- -1
l <- c(rep(0, nrow(Hgem$S)), origlb)#创建线性规划问题的下界向量 l
u <- c(rep(0, nrow(Hgem$S)), origub)#创建线性规划问题的上界向量 u
settings <- osqpSettings(max_iter = 1000000L, eps_abs = 1e-04,
eps_rel = 1e-04, adaptive_rho_interval = 50, verbose = FALSE)
#设置 OSQP 优化器的参数
model <- osqp(P, q, A, l, u, settings)#创建并解决 OSQP 优化模型
# Solve problem
res <- model$Solve()#解决优化问题并获取结果
resu[[i]] <- res
flux_vector[[i]] <- res$x#将优化结果和通量向量存储到列表中
}
close(pb)
flux_vector <- do.call(cbind, flux_vector)#将通量向量的列表合并成一个矩阵
colnames(flux_vector) <- colnames(mras)#设置通量向量矩阵的列名为 mras 矩阵的列名
rownames(flux_vector) <- Hgem$Reaction#设置通量向量矩阵的行名为反应名称
return(flux_vector)
}
建模理论知识
在代谢通量计算中，S矩阵（也称为斯托基矩阵或反应系数矩阵）和A矩阵（约束矩阵）是两个关键的数学对象，它们在构建和解决代谢网络模型中起着至关重要的作用。S矩阵提供了代谢网络的结构信息，而A矩阵则包含了确保解决方案符合生物学和实验约束的条件。通过结合这两个矩阵，可以构建和求解代谢网络模型，从而预测和理解细胞在不同条件下的代谢行为。
S矩阵（斯托基矩阵）
S矩阵是一个二维矩阵，它描述了代谢网络中所有代谢物（如酶、辅酶等）与所有代谢反应之间的关系。每一行代表一个代谢物，每一列代表一个特定的代谢反应。矩阵中的每个元素表示该代谢物在相应反应中的系数，即该代谢物在反应中被消耗或产生的数量。如果一个代谢物在某个反应中被消耗，那么相应的元素是负数；如果它是产物，则元素是正数；如果该代谢物不参与某个反应，则元素是零。S矩阵的主要作用是：
表示代谢网络的结构。 用于计算代谢物的平衡，即输入和输出的代谢物的净变化。 作为线性规划问题中的约束矩阵，用于求解代谢通量。
A矩阵（约束矩阵）
A矩阵是一个更大的矩阵，它不仅包括了S矩阵的信息，还包含了额外的约束条件，如代谢反应的方向性（可逆或不可逆）、代谢物的边界条件（如细胞边界或细胞器边界）等。A矩阵通常用于线性规划问题中，以确保求解过程中遵守这些约束条件。A矩阵的主要作用是：
将S矩阵的约束与其他类型的约束（如上界和下界）结合起来。 用于定义线性规划问题的可行解空间。 确保计算出的代谢通量符合生物学意义和实验条件。
线性规划问题
在代谢通量计算中，通常将问题表述为一个线性规划问题，目标是最大化或最小化某个特定的代谢目标（如生长速率、能量产生等）。S矩阵和A矩阵在这个过程中扮演着核心角色：
目标函数：通常由S矩阵中的特定列（代表目标代谢物）构成，通过调整这些代谢物的通量来优化目标函数。 约束条件：由A矩阵定义，确保所有代谢反应的通量都符合代谢网络的结构和细胞的生理条件。
通过解决这个线性规划问题，可以得到每个代谢反应在给定条件下的最优通量，即代谢网络在特定条件下的工作状态。
注意：由于最小化模型中总通量的总和，因此我们将得到一个简洁的通量数据输出，这意味着许多反应将接近 0 通量。（例如，有些反应应该只向前进行，但预测的通量有一个非常小的负号，可以认为接近 0 通量）
计算表达量矩阵
mean_exp=calculate_avg_exp(myseurat = sc_test_example,myident = 'Cell_type',n_bootstrap=3,seed=1)
相关源码及其解读
generate_boots
#' Generate bootstrap index data
#'
#' @param celltype colnames of single cell data.Colnames should be labeled as cell type or cluster.
#' @param n number of bootstrap
#' @import dplyr
#'
#' @return
#'
#'
#' @examples
generate_boots <- function(celltype, n) {
dt <- data.frame(cluster = celltype, id = 1:length(celltype))
index <- do.call(cbind, sapply(1:n, function(x) {
splits <- dt %>%
group_by(cluster) %>%
sample_n(dplyr::n(), replace = TRUE) %>%
ungroup() %>%
dplyr::select("id")
}))
return(index)
}
这个generate_boots函数的目的是为单细胞数据分析生成bootstrap索引。函数接受两个参数:celltype和n。celltype是一个向量，包含每个细胞的类型或簇标签；n是一个整数，表示要生成的bootstrap样本的数量。函数的输出是一个矩阵，其中每一列代表一个bootstrap样本的索引
get_ave_exp
#
#' Calculate mean expression for one bootstrap
#'
#' @param i index
#' @param myseurat single cell Seurat object.METAFlux will calculate on "data" slot
#' @param samples generated bootstrap index data
#' @param myident Seurat idents.This will be a character string indicating the grouping of the seurat object
#'
#' @return
#' @import Seurat
#'
#' @examples
get_ave_exp <- function(i, myseurat, samples,myident) {
meta.data=myseurat@meta.data[samples[,i],]
#从Seurat对象的元数据中提取出第 i 列的bootstrap样本索引对应的样本数据
sample <-myseurat@assays$RNA@counts[,samples[,i]]
#从Seurat对象的RNA矩阵中提取出第 i 列的bootstrap样本索引对应的表达数据
SeuratObject<-suppressWarnings(
CreateSeuratObject(count=sample,meta.data = meta.data))
#使用提取的样本数据和表达量数据创建一个新的Seurat对象
SeuratObject<-NormalizeData(SeuratObject,verbose = FALSE)
#新创建的Seurat对象进行数据标准化
ave<-AverageExpression(SeuratObject,group.by = myident,return.seurat = T)[["RNA"]]@data
#计算新Seurat对象的平均表达量，按照 myident 指定的分组进行，然后从结果中提取RNA检测的平均表达量数据
return(ave)
}
这个函数get_ave_exp用于计算每一次bootstrap样本的平均表达量
calculate_avg_exp
#' Calculate bootstrap mean expression for single cell data
#'
#' @param n_bootstrap number of bootstrap
#' @param seed random seed
#' @param myseurat Seurat object. METAFlux will calculate on "data" slot
#' @param myident Seurat idents for grouping.This will be a character string indicating the grouping of the seurat object
#'
#' @return mean expression data
#' @import Seurat
#' @export
#'
#' @examples
calculate_avg_exp <- function(myseurat,myident,n_bootstrap,seed) {
set.seed(seed)#保证计算的可重复性
samples=generate_boots(myseurat@meta.data[,myident],n_bootstrap)
#函数生成bootstrap索引数据。这里使用Seurat对象的样本数据和指定的分组idents来生成索引
exp <- lapply(1:n_bootstrap,get_ave_exp,myseurat,samples,myident)
#对每个bootstrap样本调用 get_ave_exp 函数计算平均表达量
exp <- do.call(cbind, exp)
#将exp中的所有平均表达量矩阵按列合并，形成一个矩阵，其中每一列代表一个bootstrap样本的平均表达量
return(exp)
}
统计理论知识
在统计分析中，bootstrap是一种重采样方法，用于估计一个统计量的分布（例如均值、中位数、方差等）。通过从原始数据集中多次有放回地随机抽取样本，可以创建多个“bootstrap样本”，然后对每个样本计算感兴趣的统计量。这个过程可以帮助评估统计量的不确定性和置信区间。这里是示例练习，所以bootstrap设置为3；真实分析过程时bootstrap应该设置为几百bootstrap样本的定义是指从原始单细胞数据集中有放回地随机抽取的子集。这个过程模拟了多次实验或测量，目的是为了估计统计量的分布，如平均表达量、方差等。
本期的内容就到这里啦，我们下期再见~
