肿瘤单细胞代谢预测实战|全流程python分析

编者按

在重构R算法这个栏目中，我们将介绍omicverse开源重构生态最新的成果以及完整的教程，欢迎关注。

本第三期教程里，由首都医科大学的林子恒给大家带来单细胞代谢分析：这是omicverse中 Metabolism 模块的内容。我们重构了scmetabolism和scFEA。

单细胞转录组测的是RNA，不是代谢物。但很多时候，代谢状态会在转录层面留下非常强的信号。比如糖酵解、氧化磷酸化、戊糖磷酸途径这些通路，如果一类细胞真的在代谢层面被激活，那么相关基因的表达也会跟着发生系统性变化。

本文将复现Xiao, Dai & Locasale 2019在Nature Communications发表的肿瘤微环境单细胞代谢图谱分析。

本文AI率20%，阅读大约需要5min。

OmicVerse团队

1. 引言

肿瘤不是一团均匀的细胞。

一个真实的肿瘤组织里，有恶性细胞、成纤维细胞、内皮细胞、T细胞、B细胞、巨噬细胞、树突细胞。它们待在同一个微环境里，但代谢方式完全不一样。

如果你用bulk RNA-seq去看，得到的是所有细胞混在一起的平均值。这个平均值当然有用，但它会把很多真正有意思的东西抹平。

比如我们真正想问的是：

• 哪些细胞的糖酵解最强？
• 恶性细胞是不是整体代谢更活跃？
• 氧化磷酸化是不是只在某一群肿瘤细胞里升高？
• 缺氧和糖酵解是不是在单细胞层面相关？
• 不同免疫/基质细胞是否有各自的代谢偏好？

这些问题都需要单细胞分辨率。

2019年，Xiao, Dai和Locasale在Nature Communications发表了文章，系统分析了肿瘤微环境的单细胞代谢图谱。他们的一个主要结论是：在头颈鳞癌数据中，恶性细胞上调的代谢通路数量最多，肿瘤内部还存在明显的代谢异质性。

2. Metabolism模块能做什么

OmicVerse里统一通过这个接口做代谢分析：

ov.single.Metabolism(adata, method='scmetabolism')

目前这个模块里有三类思路：

在本期教程，我们主要使用 scmetabolism 来计算85条KEGG代谢通路活性，最后再用 scFEA 做一个代谢通量的补充。

scmetabolism是复现的R包，我们在确保计算得分一致性的前提下，对计算速度也进行了较大的优化，最终得到的是一个“细胞 × 代谢通路”的矩阵，我们可以用这个矩阵进行差异分析、UMAP可视化、热图和细胞类型比较等内容。

3. 数据读入

首先载入OmicVerse并设置绘图风格：

import omicverse as ov

ov.plot_set()

这里使用的是Puram et al. 2017的头颈部鳞状细胞癌数据，也就是GSE103322。这个数据包含19位患者的5,578个单细胞，也是Xiao 2019那篇文章分析的数据。

adata = ov.datasets.metabolism_hnsc()
adata

AnnData object with n_obs × n_vars = 5578 × 23686
    obs: 'patient', 'malignant', 'celltype', 'lymph_node', 'maxima_enzyme'
    uns: 'dataset', 'expression_units'

我们看一下细胞类型组成：

adata.obs['celltype'].value_counts()

celltype
Malignant      2215
Fibroblast     1440
T cell         1237
Endothelial     260
B cell          138
Mast            120
Macrophage       98
Dendritic        51
myocyte          19
Name: count, dtype: int64

在这个数据中，恶性细胞数量足够多，免疫和基质细胞也比较完整，可以直接比较肿瘤细胞和微环境细胞之间的代谢差异。

4. 构建UMAP

由于这个数据已经是标准化后的表达矩阵，所以我们直接选高变基因，然后做scale、PCA、邻居图和UMAP标准的单细胞分析流程。

ov.pp.highly_variable_genes(adata, n_top_genes=2000)
adata.raw = adata

hvg = adata[:, adata.var.highly_variable].copy()
ov.pp.scale(hvg)
ov.pp.pca(hvg, layer='scaled', n_pcs=30)

ov.pp.neighbors(
    hvg,
    n_neighbors=15,
    n_pcs=30,
    use_rep='scaled|original|X_pca',
)
ov.pp.umap(hvg)
adata.obsm['X_umap'] = hvg.obsm['X_umap']

首先看一下UMAP上细胞类型和恶性状态的分布：

ov.pl.embedding(
    adata,
    basis='X_umap',
    color=['celltype', 'malignant'],
    frameon='small',
    wspace=0.5,
)

可以看到恶性细胞和其他细胞在UMAP上基本可以分开。后面我们把代谢通路分数投到这个UMAP上，就能看出哪些区域的代谢活性更强。

5. 计算代谢通路活性

接下来运行 scmetabolism。

met = ov.single.Metabolism(adata, method='scmetabolism')
met.run(score_method='AUCell', metabolism_type='KEGG')

这里用AUCell计算所有KEGG代谢通路在每个细胞里的活性。结果会写入：

adata.obsm['X_metabolism']

看一下矩阵大小：

adata.obsm['X_metabolism'].shape

(5578, 85)

也就是说，我们现在得到了一个5,578个细胞 × 85条KEGG代谢通路的活性矩阵。

这一步其实就是后面所有代谢分析的基础。

6. 恶性细胞上调最多代谢通路

Xiao 2019的一个主要结论是：恶性细胞上调的代谢通路数量最多。

我们可以用 ov.single.differential_metabolism 比较恶性细胞和非恶性细胞的代谢通路活性，这个算法调用的也是omicverse重构生态的成果py-limma包。

deg = ov.single.differential_metabolism(
    adata,
    groupby='malignant',
    group1='Malignant',
)

up = deg.query('padj < 0.05 and log2fc > 0')
print(f'{len(up)} of {len(deg)} KEGG pathways up-regulated in malignant cells')

53 of 85 KEGG pathways up-regulated in malignant cells

在这个数据中，我们发现在85条KEGG代谢通路里，有53条在恶性细胞中显著上调。

我们看一下排名靠前的通路：

deg.head(8)

我觉得这个结果还是比较符合肿瘤代谢的背景知识。糖酵解、戊糖磷酸途径、氧化磷酸化、嘧啶代谢、嘌呤代谢都排在前面。

接着我们把不同细胞类型分别拿出来比较，统计每类细胞显著上调的代谢通路数量：

import pandas as pd

n_up = {
    ct: len(
        ov.single.differential_metabolism(
            adata,
            groupby='celltype',
            group1=ct,
        ).query('padj < 0.05 and log2fc > 0')
    )
for ct in ['Malignant', 'Fibroblast', 'Macrophage', 'B cell', 'T cell']
}

pd.Series(n_up).sort_values(ascending=False)

Malignant     53
Macrophage    38
Fibroblast    27
B cell        17
T cell        15
dtype: int64

可视化：