空间转录组数据注释分析：Cell2location反卷积(Nature Biotechnology IF: 33.1)
cell2location 于2022年发表于Nature Biotechnology [IF: 33.1]，是一种基于分层贝叶斯统计模型的空间转录组反卷积分析工具，旨在用于精确解析组织中不同细胞类型的空间分布及其丰度。它需要使用单细胞数据作为参考，将单细胞水平的细胞类型注释信息“映射”到空间位点中，揭示复杂组织中细胞类型、状态与空间位置的关联，为理解组织微环境提供重要线索。
这个方法需要性能比较高的计算资源配置，如果有GPU显卡会大大的缩短运行建模的时间，如果你没有，可以考虑生信技能树的共享服务器，只要800一年，内存2T，配置20G显卡，购买详情可看：满足你生信分析计算需求的低价解决方案
软件原理图如下
本次学习来自官网的教程：
.html
0.环境配置
conda 安装：
代码语言：javascript代码运行次数：0运行复制conda create -n cell2location -y python=3.9
conda activate cell2location
pip install -i jupyterlab
# 永久设置镜像
pip config set global.index-url /
# /home/zhangj/.config/pip/pip.conf
pip install scanpy
pip install cell2location[tutorials]
pip install pandas==2.1.1
加载：
代码语言：javascript代码运行次数：0运行复制import os
# this line should go before importing cell2location
os.environ["THEANO_FLAGS"] = 'device=cuda,floatX=float32,force_device=True'
# Loading packages
import sys
import scanpy as sc
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import cell2location
from matplotlib import rcParams
rcParams['pdf.fonttype'] = 42 # enables correct plotting of text for PDFs
1.加载Visium空间数据
1.1 数据介绍
参考的单细胞数据来自人的次级淋巴器官包括34种细胞类型，空间数据来自10X Genomics的人类淋巴结Visium dataset，为10X Space Ranger软件的标准输出文件，来自 Kleshchevnikov et al (section 4, Fig 4)。
Visium空间数据下载地址为10X Genomics 官网: .0.0/V1_Human_Lymph_Node。格式为10X Space Ranger软件的输出，去网站下载到本地后导入：
代码语言：javascript代码运行次数：0运行复制adata_vis = sc.read_visium("../data/V1_Human_Lymph_Node/1.0.0/", count_file='filtered_feature_bc_matrix.h5', load_images=True)
adata_vis
adata_vis.var_names
# rename genes to ENSEMBL ID for correct matching between single cell and spatial data
adata_vis.obs['sample'] = list(adata_vis.uns['spatial'].keys())[0]
adata_vis.var['SYMBOL'] = adata_vis.var_names
adata_vis.var.set_index('gene_ids', drop=True, inplace=True)
adata_vis
adata_vis.var_names
注意，数据 data/V1_Human_Lymph_Node/1.0.0/ 目录结构为：
1.2 去除线粒体基因
线粒体编码的基因(基因名称以前缀 MT-或 MT-开头)与空间mapping无关，因为它们的表达代表单细胞和细胞核数据中的技术artifacts，而不是线粒体的生物丰度。然而，这些基因在每个spot组成15-40% 的 mRNA 。因此，为了避免绘制artifacts，强烈建议去除线粒体基因。
代码语言：javascript代码运行次数：0运行复制# find mitochondria-encoded (MT) genes
adata_vis.var['MT_gene'] = [gene.startswith('MT-') for gene in adata_vis.var['SYMBOL']]
# remove MT genes for spatial mapping (keeping their counts in the object)
adata_vis.obsm['MT'] = adata_vis[:, adata_vis.var['MT_gene'].values].X.toarray()
adata_vis = adata_vis[:, ~adata_vis.var['MT_gene'].values]
2.加载单细胞参考集数据
由于患者供体年龄的原因，已发表的淋巴结scRNA-seq数据集通常缺乏生发中心相关免疫细胞群的充分代表。因此，作者将综合淋巴结、脾脏和扁桃体的scRNA-seq数据集纳入单细胞参考，以确保捕获了空间转录组数据集中可能存在的免疫细胞状态的全部多样性。
2.1 导入单细胞数据
下载地址：.h5ad
代码语言：javascript代码运行次数：0运行复制# Read scrna data
adata_ref = sc.read('../data/V1_Human_Lymph_Node/1.0.0/sc.h5ad')
# 将数据的genes转为ENSEMBL ID，与空间数据对应同样的id
adata_ref.var['SYMBOL'] = adata_ref.var.index
# rename 'GeneID-2' as necessary for your data
adata_ref.var.set_index('GeneID-2', drop=True, inplace=True)
adata_ref.var
2.2 very permissive genes选择
在估计参考细胞类型特征之前，建议进行very permissive genes选择。与标准的高可变基因选择相比，这种方法保留了罕见标记基因，同时删除了大多数无信息的基因。默认参数为：cell_count_cutoff=5, cell_percentage_cutoff2=0.03, nonz_mean_cutoff=1.12
在实战项目过程中，可以适当增加这些阈值以便过滤更多基因，推荐 cell_count_cutoff=5，cell_percentage_cutoff2和nonz_mean_cutoff可以适当增加到基因个数为 8k-16k 个基因。
from cell2location.utils.filtering import filter_genes
selected = filter_genes(adata_ref, cell_count_cutoff=5, cell_percentage_cutoff2=0.03, nonz_mean_cutoff=1.12)
# filter the object
adata_ref = adata_ref[:, selected].copy()
过滤掉的基因示意图：
- 每一个点表示一个基因
- x轴：检测到该基因的细胞中的平均 RNA count值
- y轴：表达该基因的细胞数量
橙色矩形突出显示了基于x轴和y轴的组合而排除的基因。
3.参考细胞类型signatures评估
根据scRNA-seq数据估计特征，考虑到批此效应，使用负二项回归模型进行分析。
3.1 首先，为回归模型准备一个数据对象:
代码语言：javascript代码运行次数：0运行复制# prepare anndata for the regression model
cell2location.models.RegressionModel.setup_anndata(adata=adata_ref,
# 10X reaction / sample / batch
batch_key='Sample',
# cell type, covariate used for constructing signatures
labels_key='Subset',
# multiplicative technical effects (platform, 3' vs 5', donor effect)
categorical_covariate_keys=['Method']
)
# create the regression model
from cell2location.models import RegressionModel
mod = RegressionModel(adata_ref)
# view anndata_setup as a sanity check
mod.view_anndata_setup()
显示的内容如下：
代码语言：javascript代码运行次数：0运行复制Anndata setup with scvi-tools version 1.1.6.post2.
Setup via `RegressionModel.setup_anndata` with arguments:
{
│ 'layer': None,
│ 'batch_key': 'Sample',
│ 'labels_key': 'Subset',
│ 'categorical_covariate_keys': ['Method'],
│ 'continuous_covariate_keys': None
}
Summary Statistics
┏━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━┓
┃ Summary Stat Key ┃ Value ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━┩
│ n_batch │ 23 │
│ n_cells │ 73260 │
│ n_extra_categorical_covs │ 1 │
│ n_extra_continuous_covs │ 0 │
│ n_labels │ 34 │
│ n_vars │ 10237 │
└──────────────────────────┴───────┘
Data Registry
┏━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Registry Key ┃ scvi-tools Location ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│ X │ adata.X │
│ batch │ adata.obs['_scvi_batch'] │
│ extra_categorical_covs │ adata.obsm['_scvi_extra_categorical_covs'] │
│ ind_x │ adata.obs['_indices'] │
│ labels │ adata.obs['_scvi_labels'] │
└────────────────────────┴────────────────────────────────────────────┘
batch State Registry
┏━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓
┃ Source Location ┃ Categories ┃ scvi-tools Encoding ┃
┡━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩
│ adata.obs['Sample'] │ 4861STDY7135913 │ 0 │
│ │ 4861STDY7135914 │ 1 │
│ │ 4861STDY7208412 │ 2 │
│ │ 4861STDY7208413 │ 3 │
│ │ 4861STDY7462253 │ 4 │
│ │ 4861STDY7462254 │ 5 │
│ │ 4861STDY7462255 │ 6 │
│ │ 4861STDY7462256 │ 7 │
│ │ 4861STDY7528597 │ 8 │
│ │ 4861STDY7528598 │ 9 │
│ │ 4861STDY7528599 │ 10 │
│ │ 4861STDY7528600 │ 11 │
│ │ BCP002_Total │ 12 │
│ │ BCP003_Total │ 13 │
│ │ BCP004_Total │ 14 │
│ │ BCP005_Total │ 15 │
│ │ BCP006_Total │ 16 │
│ │ BCP008_Total │ 17 │
│ │ BCP009_Total │ 18 │
│ │ Human_colon_16S7255677 │ 19 │
│ │ Human_colon_16S7255678 │ 20 │
│ │ Human_colon_16S8000484 │ 21 │
│ │ Pan_T7935494 │ 22 │
└─────────────────────┴────────────────────────┴─────────────────────┘
labels State Registry
┏━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓
┃ Source Location ┃ Categories ┃ scvi-tools Encoding ┃
┡━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩
│ adata.obs['Subset'] │ B_Cycling │ 0 │
│ │ B_GC_DZ │ 1 │
│ │ B_GC_LZ │ 2 │
│ │ B_GC_prePB │ 3 │
│ │ B_IFN │ 4 │
│ │ B_activated │ 5 │
│ │ B_mem │ 6 │
│ │ B_naive │ 7 │
│ │ B_plasma │ 8 │
│ │ B_preGC │ 9 │
│ │ DC_CCR7+ │ 10 │
│ │ DC_cDC1 │ 11 │
│ │ DC_cDC2 │ 12 │
│ │ DC_pDC │ 13 │
│ │ Endo │ 14 │
│ │ FDC │ 15 │
│ │ ILC │ 16 │
│ │ Macrophages_M1 │ 17 │
│ │ Macrophages_M2 │ 18 │
│ │ Mast │ 19 │
│ │ Monocytes │ 20 │
│ │ NK │ 21 │
│ │ NKT │ 22 │
│ │ T_CD4+ │ 23 │
│ │ T_CD4+_TfH │ 24 │
│ │ T_CD4+_TfH_GC │ 25 │
│ │ T_CD4+_naive │ 26 │
│ │ T_CD8+_CD161+ │ 27 │
│ │ T_CD8+_cytotoxic │ 28 │
│ │ T_CD8+_naive │ 29 │
│ │ T_TIM3+ │ 30 │
│ │ T_TfR │ 31 │
│ │ T_Treg │ 32 │
│ │ VSMC │ 33 │
└─────────────────────┴──────────────────┴─────────────────────┘
extra_categorical_covs State Registry
┏━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓
┃ Source Location ┃ Categories ┃ scvi-tools Encoding ┃
┡━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩
│ adata.obs['Method'] │ 3GEX │ 0 │
│ │ 5GEX │ 1 │
│ │ │ │
└─────────────────────┴────────────┴─────────────────────┘
3.2 训练模型
estimate the reference cell type signatures，这一步骤没有GPU将非常耗时
代码语言：javascript代码运行次数：0运行复制## 如果你有GPU超算，需要在前面导入cell2location模块的时候设置
import os
# this line should go before importing cell2location
os.environ["THEANO_FLAGS"] = 'device=cuda,floatX=float32,force_device=True'
# 训练模型
mod.train(max_epochs=250)
# 我这里使用GPU，运行了18m24s
3.3 检查模型是否需要更多的训练
可以绘制训练期间的 ELBO 损失历史，删除前20个 epoch。这个图在训练结束时应该有一个下降的趋势并趋于平稳。如果 max_epochs 仍在减小，则需要增加max_epochs：
代码语言：javascript代码运行次数：0运行复制mod.plot_history(20)
plt.savefig(results_folder+'01-mod.plot_history.png') # 保存为PNG格式
训练结果图：
横轴为epochs，纵轴为ELBO loss，这个图在训练结束时应该有一个下降的趋势并趋于平稳。
3.4 保存训练后估计的细胞丰度
输出训练的结果：
代码语言：javascript代码运行次数：0运行复制# In this section, we export the estimated cell abundance (summary of the posterior distribution).
adata_ref = mod.export_posterior(
adata_ref, sample_kwargs={'num_samples': 1000, 'batch_size': 2500}
)
# Save model
mod.save('./reference_signatures/', overwrite=True)
adata_file = "./reference_signatures/sc.h5ad"
adata_ref.write(adata_file)
adata_file保存的数据后续还可以重新加载回来：
代码语言：javascript代码运行次数：0运行复制adata_file = "./reference_signatures/sc.h5ad"
adata_ref = sc.read_h5ad(adata_file)
mod = cell2location.models.RegressionModel.load('./reference_signatures/', adata_ref)3.4 检查评估结果
评估标准：
1.Reconstruction accuracy图：这个二维直方图的大部分观测值应该沿着有噪声的对角线
2.图2：由于批处理效应，估计的表达特征与每个聚类的平均表达特征不同。对于不受批处理影响的scRNA-seq数据集(本数据集)，可以使用聚类平均表达来代替使用模型估计signature特征。当这个图与对角线图非常不同时(例如，y轴上的值非常低，到处都是密度)，这表明特征估计存在问题。
代码语言：javascript代码运行次数：0运行复制mod.plot_QC()评估图：
3.5 提取参考细胞类型特征
提取参考细胞类型特征为pd.DataFrame：空间转录组数据的细胞类型映射只需要估计得到的signatures
代码语言：javascript代码运行次数：0运行复制# export estimated expression in each cluster
if 'means_per_cluster_mu_fg' in adata_ref.varm.keys():
inf_aver = adata_ref.varm['means_per_cluster_mu_fg'][[f'means_per_cluster_mu_fg_{i}'
for i in adata_ref.uns['mod']['factor_names']]].copy()
else:
inf_aver = adata_ref.var[[f'means_per_cluster_mu_fg_{i}'
for i in adata_ref.uns['mod']['factor_names']]].copy()
inf_aver.columns = adata_ref.uns['mod']['factor_names']
inf_aver.iloc[0:12, 0:12]这个就是用单细胞数据得到的细胞类型signature了，用于空转数据的解卷积。
4.运行Cell2location
利用单细胞signature矩阵，对空转数据进行注释，首先提取单细胞与空转数据共同基因：
代码语言：javascript代码运行次数：0运行复制# find shared genes and subset both anndata and reference signatures
intersect = np.intersect1d(adata_vis.var_names, inf_aver.index)
adata_vis = adata_vis[:, intersect].copy()
inf_aver = inf_aver.loc[intersect, :].copy()
# prepare anndata for cell2location model
cell2location.models.Cell2location.setup_anndata(adata=adata_vis, batch_key="sample")4.1 定义参数值
在进行空转数据注释之前，需要定义两个比较重要的参数
N_cells_per_location：每个spot 的细胞数，这个可以通过配对的 histology images进行评估，或者 通过squdipy 进行细胞分割，然后估计细胞数、或者通过 10X loupe browser 手动数
detection_alpha：建议用 20，因为从已发表的研究来看，人体组织在实验中受到的技术影响较大。如果是使用的小鼠等受技术影响低的数据集，那可以用 200。当然了，最好是两个都试试。
更多详细的说明参考：the flow diagram and the note，.pdf
本教程：N_cells_per_location=30，detection_alpha=20
# create and train the model
mod = cell2location.models.Cell2location(
adata_vis, cell_state_df=inf_aver,
# the expected average cell abundance: tissue-dependent
# hyper-prior which can be estimated from paired histology:
N_cells_per_location=30,
# hyperparameter controlling normalisation of
# within-experiment variation in RNA detection:
detection_alpha=20
)
mod.view_anndata_setup()
4.2 Training cell2location
此步骤也非常耗时：在使用GPU的情况下，总共运行了42m 34.3s
代码语言：javascript代码运行次数：0运行复制mod.train(max_epochs=30000,
# train using full data (batch_size=None)
batch_size=None,
# use all data points in training because
# we need to estimate cell abundance at all locations
train_size=1
)
# plot ELBO loss history during training, removing first 100 epochs from the plot
mod.plot_history(1000)
plt.legend(labels=['full data training']);
plt.savefig('03-mod.plot_history.png') # 保存为PNG格式
横轴为epochs，纵轴为ELBO loss，这个图在训练结束时应该有一个下降的趋势并趋于平稳。
4.3 保存结果
训练完之后，将每个 spot 估计的细胞丰度，保存到空转数据中。并且保存模型：
代码语言：javascript代码运行次数：0运行复制# In this section, we export the estimated cell abundance (summary of the posterior distribution).
adata_vis = mod.export_posterior(
adata_vis, sample_kwargs={'num_samples': 1000, 'batch_size': mod.adata.n_obs}
)
mod.save("./", overwrite=True)
# 生成 model.pt
# Save anndata object with results
adata_file = "./sp.h5ad"
adata_vis.write(adata_file)
adata_file
数据后续又可以被加载回来：
代码语言：javascript代码运行次数：0运行复制adata_file = "./sp.h5ad"
adata_vis = sc.read_h5ad(adata_file)
mod = cell2location.models.Cell2location.load("./", adata_vis)
4.4 mapping质量评估
检查重建精度以评估映射是否存在任何问题，图应该大致是对角线的。
代码语言：javascript代码运行次数：0运行复制mod.plot_QC()
plt.savefig('04-mod.plot_QC.png') # 保存为PNG格式
结果图：
5.在空间坐标上可视化评估的细胞丰度
现在终于到最后一步啦，可以看看空间切片上不同spot中的细胞类型。
5.1 查看具体的每个spot上的细胞丰度：
预测的结果中有这几种：means_cell_abundance_w_sf, stds_cell_abundance_w_sf, q05_cell_abundance_w_sf, q95_cell_abundance_w_sf
q05_cell_abundance_w_sf：表示每个位置上每种细胞类型的预期细胞数量，包括可能的小数（fractional values），即细胞类型的密度- 5% 分位数：这个值是后验分布的 5% 分位数，意味着模型对这个值有较高的置信度，认为至少有这么多细胞存在
adata_vis.obsm
adata_vis.obsm['q05_cell_abundance_w_sf']
adata_vis.obsm['q95_cell_abundance_w_sf']
5.2 在空间上可视化
选择 q05_cell_abundance_w_sf的结果：
# add 5% quantile, representing confident cell abundance, 'at least this amount is present',
# to adata.obs with nice names for plotting
adata_vis.obs[adata_vis.uns['mod']['factor_names']] = adata_vis.obsm['q05_cell_abundance_w_sf']
# select one slide
from cell2location.utils import select_slide
slide = select_slide(adata_vis, 'V1_Human_Lymph_Node')
# plot in spatial coordinates
with mpl.rc_context({'axes.facecolor': 'black','figure.figsize': [4.5, 5]}):
sc.pl.spatial(slide, cmap='magma',
# show first 8 cell types
color=['B_Cycling', 'B_GC_LZ', 'T_CD4+_TfH_GC', 'FDC','B_naive', 'T_CD4+_naive', 'B_plasma', 'Endo'],
ncols=4, size=1.3,img_key='hires',
# limit color scale at 99.2% quantile of cell abundance
vmin=0, vmax='p99.2'
)
plt.savefig('05-sc.pl.spatial.png') # 保存为PNG格式
这里选择了8种细胞进行可视化：
也可以选择多个细胞类型展示在一张切片上：这里选择了三种细胞 T_CD4+_naive，B_naive，FDC 进行共定位
代码语言：javascript代码运行次数：0运行复制# Now we use cell2location plotter that allows showing multiple cell types in one panel
from cell2location.plt import plot_spatial
# select up to 6 clusters
clust_labels = ['T_CD4+_naive', 'B_naive', 'FDC']
clust_col = ['' + str(i) for i in clust_labels] # in case column names differ from labels
clust_col
slide = select_slide(adata_vis, 'V1_Human_Lymph_Node')
with mpl.rc_context({'figure.figsize': (15, 15)}):
fig = plot_spatial(
adata=slide,
# labels to show on a plot
color=clust_col, labels=clust_labels,
show_img=True,
# 'fast' (white background) or 'dark_background'
style='fast',
# limit color scale at 99.2% quantile of cell abundance
max_color_quantile=0.992,
# size of locations (adjust depending on figure size)
circle_diameter=6,
colorbar_position='right'
)
plt.savefig('06-sc.pl.spatial.png') # 保存为PNG格式
结果：
这个跟其他空间工具解卷积结果来看，对于同一个spot如果有多种细胞类型，这里展示的是一个 颜色的叠加而不是饼图，嗯，我感觉还是饼图可能会更直观一点，这里如果颜色选取有强弱很容易造成视觉上的错觉吧。
软件版本
代码语言：javascript代码运行次数：0运行复制cell2location.__version__
# '0.1.4'
发布评论