在前面的教程中，我们已经详细介绍了METAflux软件分析Bulk RNA_seq的工作流程，感兴趣的小伙伴可以阅读：METAflux | 实操教程2-Bulk RNA-seq数据分析步骤拆解，今天的推文将详细介绍METAflux软件分析scRNA-seq的工作流程。

1.准备工作

METAflux首先将整个肿瘤微环境建模为一个整体，以解释不同组之间的代谢相互作用。在这里，我们不计算单个细胞上通量。相反，我们在簇/细胞类型水平上对通量进行建模。为了估计代谢通量，我们需要提供：

1. 逐个细胞表达矩阵的基因。在使用 METAFlux 之前，应将输入基因表达矩阵标准化（例如，对数转换等）。METAflux 不会对表达式数据执行任何归一化。
2. 每个细胞聚类或按细胞类型分配。
3. 簇/细胞类型分数.这些比例应从实验中获取或使用 CIBERSORTx 从匹配的批量数据中计算 [2]。但是，大多数数据集没有此类信息。仍然可以使用单细胞数据中直接计算的聚类（组）分数。然而，有必要进一步研究来评估这些发现，因为这些比例可能会因组织解离过程中引入的偏差而偏离事实。

通常，每个单细胞数据集都包含多个样本（患者）。为每个受试者分别估计每个簇/细胞类型的比例会更自然（可以把样本拆分出来逐一分析），因为 METAFlux 将对同一 TME 内的相互作用进行建模。然而，根据特定的生物学问题和数据质量，人们仍然可以研究受试者的平均效应（即，使用整合后的所有数据估计每个簇/细胞类型的比例）。

注意，这里的定义是细胞类型和细胞簇，对于簇其实就可以延展到cell state。具体来说就是用富集后的通路对其打分定义state或者结合GRN。

2.主要分析流程

加载必要的数据库

library(METAFlux)
#load medium file.Please read '3.Step by step breakdown:Bulk RNA-seq pipeline section 3.2.3' for more details
data("human_blood")
data("cell_medium")
#load underlying GEM.Please read '3.Step by step breakdown:Bulk RNA-seq pipeline section 3.2.2' for more details
data("human_gem")

使用内置的测试数据集

data("sc_test_example")

#check cell types, here we want to model cell-type level flux
table(sc_test_example$Cell_type)

内置函数来计算 bootstrap 样本的平均表达式。还有许多其他方法可以聚合样本，例如，中位数或几何平均值。因此，我们可以计算自己的 “平均表达概况” 以输入下一步，但布局应采用相同的数据格式。还是我们熟悉的行为基因，列为样本矩阵

#Calculate the mean expression for bootstrap samples
#For the sake of demonstration, we only set the number of bootstraps to 3. 
#In real analysis, the number of bootstraps should be much higher(e.g. 100, 200....)
mean_exp=calculate_avg_exp(myseurat = sc_test_example,myident = 'Cell_type',n_bootstrap=3,seed=1)

head(mean_exp)

打出来检查一下

##          B lymphocytes Epithelial cells Myeloid cells T lymphocytes
## A1BG         0.5903632        0.0000000     0.7048261     0.6375582
## A1BG-AS1     0.1846133        0.0000000     0.2007492     0.1699989
## A1CF         0.0000000        0.0000000     0.0000000     0.0000000
## A2M          0.0000000        0.1032028     1.1101576     0.0000000
## A2M-AS1      0.0000000        0.0000000     0.0000000     0.1172836
## A2ML1        0.0000000        0.0000000     0.0000000     0.0000000
##          B lymphocytes Epithelial cells Myeloid cells T lymphocytes
## A1BG          1.082775       0.00000000     0.6059942     0.5785297
## A1BG-AS1      0.000000       0.00000000     0.1320018     0.0985478
## A1CF          0.000000       0.00000000     0.0000000     0.0000000
## A2M           0.000000       0.06437396     1.0537790     0.0000000
## A2M-AS1       0.000000       0.00000000     0.0000000     0.1172836
## A2ML1         0.000000       0.00000000     0.0000000     0.0000000
##          B lymphocytes Epithelial cells Myeloid cells T lymphocytes
## A1BG          1.207643        0.0000000     0.5241961     0.2085221
## A1BG-AS1      0.000000        0.0000000     0.2405239     0.0000000
## A1CF          0.000000        0.0000000     0.0000000     0.0000000
## A2M           0.000000        0.1818827     0.9224259     0.0000000
## A2M-AS1       0.000000        0.0000000     0.0000000     0.0000000
## A2ML1         0.000000        0.0000000     0.0000000     0.0000000

使用基因 - 蛋白质反应（GPR）从数据中计算单个样本归一化 MRAS。此步骤与bulk RNA-seq流程是相同的

#calculate metabolic reaction scores
scores<-calculate_reaction_score(data=mean_exp)

head(scores,4)

依旧检查一下，保持好习惯

##          B lymphocytes Epithelial cells Myeloid cells T lymphocytes
## HMR_3905  0.0001110376     0.0016033661  6.089889e-04  0.0011867417
## HMR_3907  0.0018399207     0.0011775454  1.781516e-03  0.0005612051
## HMR_4097  0.0000000000     0.0001369900  7.021776e-05  0.0000000000
## HMR_4099  0.0005489067     0.0009259649  3.490022e-04  0.0001699136
##          B lymphocytes.1 Epithelial cells.1 Myeloid cells.1 T lymphocytes.1
## HMR_3905    0.0001052063       0.0017716880    0.0007462208    0.0009607519
## HMR_3907    0.0019568699       0.0012993935    0.0020043234    0.0006483234
## HMR_4097    0.0000000000       0.0001750853    0.0001416029    0.0000000000
## HMR_4099    0.0003529752       0.0010112861    0.0004172000    0.0003275986
##          B lymphocytes.2 Epithelial cells.2 Myeloid cells.2 T lymphocytes.2
## HMR_3905    0.0001941826       0.0019331392    0.0007114677    0.0006899241
## HMR_3907    0.0025411176       0.0011940906    0.0018005957    0.0003950606
## HMR_4097    0.0000000000       0.0003202937    0.0001778084    0.0000000000
## HMR_4099    0.0030136719       0.0008887043    0.0006402115    0.0001799770

3.计算代谢流

计算代谢通量。请注意，_分数参数输入_的顺序应与 scores 结果中的顺序一致。例如，细胞类型 1（B 淋巴细胞）是第一列评分，细胞类型 2（上皮细胞）是第二列评分，依此类推。更重要的是，分数需要求和为 1。

#calculate the fractions of celltype/clusters
round(table(sc_test_example$Cell_type)/nrow(sc_test_example@meta.data),1)

##    B lymphocytes Epithelial cells    Myeloid cells    T lymphocytes 
##              0.1              0.3              0.3              0.3

#calculate flux,here we used human blood as our medium, but please change to cell line medium if you sample is from cell line.
#num_cell:number of cell types/clusters=4
flux=compute_sc_flux(num_cell = 4,fraction =c(0.1,0.3,0.3,0.3),fluxscore=scores,medium = human_blood)

4.数据整合和可视化

• 预测通量数据的总维度将为(num_cell∗13082+1648）∗ number_of_bootstrap.行名称用 celltype 1、celltype2 等注释，以指示该反应所指的细胞类型/簇。带有 “external_medium” 的行名是整个 TME 与外部环境交换（摄取或释放）代谢物的反应。这些反应是指整个肿瘤群落，而不是一种特定的细胞类型。带有 “internal_medium ” 的行名是特定细胞类型/簇与 TME 交换（摄取或释放）代谢物的反应。

dim(flux)
## [1] 53976     3
#(4*13082+1648)*3=53976*3

• 通量的sign具有生物学意义。这在 bulk pipeline 中已经提到过，这里的逻辑是相同的。对于营养吸收/释放情况（由营养查找文件中的 1648 个交换反应表示的情况），正值表示化合物的释放，负值表示化合物的吸收。对于其他反应，正值表示净通量是向前的，负的通量表示净通量是反向的。绝对值表示磁通量的大小
• 细胞类型特异性通量是平均通量。细胞类型特异性通量是每种细胞类型的平均通量或单位通量。更具体地说，如果反应 A 的细胞类型 1 通量等于 0.2，则仅意味着细胞类型 1 中平均细胞的反应 A 通量为 0.2。external_medium 通量不是平均代谢通量，而是加权总代谢通量。external_medium 和 internal_medium 之间的通量可以用以下这样交换反应方程解释：

Proportioncelltype1⋅Intenal−MediumCelltype1exchange−reactioni+Proportioncelltype2⋅Intenal−MediumCelltype2exchange−reactioni+……..=Entenal−Mediumexchange−reactioni

• 我们可以提取我们感兴趣的反应。同样，如果对代谢物的吸收或释放非常感兴趣，可以搜索“营养查找文件”以获得代谢物交换反应 ID（请阅读“批量分步分解第 2.3.1 节”了解更多详情）。例如，如果我们想知道葡萄糖摄取，我们需要搜索 'glucose' 来获得葡萄糖摄取反应，即 HMR_9034。接下来，我们可以通过以下方式提取数据：

data("nutrient_lookup_files")
glucose<-data.frame(glucose=flux[grep("HMR_9034",rownames(flux)),])
glucose

##                                       glucose.V1   glucose.V2   glucose.V3
## celltype 1 internal_medium HMR_9034 -0.010137361 -0.016898285 -0.006488233
## celltype 2 internal_medium HMR_9034 -0.001911976 -0.002269979 -0.003828839
## celltype 3 internal_medium HMR_9034 -0.003646470 -0.006834260 -0.005491176
## celltype 4 internal_medium HMR_9034 -0.001676353 -0.001327100 -0.002836066
## external_medium HMR_9034            -0.003184176 -0.004819230 -0.004295648

#if needed we can apply cubic root normalization to normalize the scores
cbrt <- function(x) {
    sign(x) * abs(x)^(1/3)
}

# glucose=cbrt(glucose)

综上所述，看到 Glucose.V1、Glucose.v2、Glucose.V3 是 3 个 bootstrap 样本。人们可以查看所有 bootstraps 的分布，以比较 b、上皮、髓系和 T 细胞的营养吸收概况。

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))

• “celltype 1 internal_medium HMR_9034”是 1 型细胞（指示例中的 B 细胞）的葡萄糖摄取通量。
• “细胞类型 2 internal_medium HMR_9034”是 2 型细胞（指示例中的上皮细胞）的葡萄糖摄取通量。
• “celltype 3 internal_medium HMR_9034”是 2 型细胞（指示例中的髓系细胞）的葡萄糖摄取通量。
• “4 型 internal_medium HMR_9034”是 2 型细胞（指示例中的 T 细胞）的葡萄糖摄取通量。
• “external_medium HMR_9034” 是整个肿瘤群落的葡萄糖摄取通量（在示例中指组合的 TME 群落）。

关注公众号，下回更新不迷路