在前面的教程中,我们已经详细介绍了METAflux软件分析Bulk RNA_seq的工作流程,感兴趣的小伙伴可以阅读:METAflux | 实操教程2-Bulk RNA-seq数据分析步骤拆解,今天的推文将详细介绍METAflux软件分析scRNA-seq的工作流程。
1.准备工作
METAflux首先将整个肿瘤微环境建模为一个整体,以解释不同组之间的代谢相互作用。在这里,我们不计算单个细胞上通量。相反,我们在簇/细胞类型水平上对通量进行建模。为了估计代谢通量,我们需要提供:
逐个细胞表达矩阵的基因。在使用 METAFlux 之前,应将输入基因表达矩阵标准化(例如,对数转换等)。METAflux 不会对表达式数据执行任何归一化。 每个细胞聚类或按细胞类型分配。 簇/细胞类型分数.这些比例应从实验中获取或使用 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 群落)。