ixxmu
commented
2 months ago 整合单细胞和空转数据多种方法之Cell2location by 生信随笔
除了我们上期介绍过的CellTrek算法【整合单细胞和空转数据多种方法之CellTrek】，还有非常多的算法用于整合单细胞和空转数据，如何快速系统的了解更多的主流整合算法呢？这么多种算法我们又应该选择哪种呢？答案是查阅综述+大量实践！
例如，2022年发表在Nat Methods上的这篇空转Benchmarking文章（DOI: 10.1038/s41592-022-01480-9），来自瞿昆老师团队。文章系统比较了16个空间转录组数据分析算法：
结果显示，cell2location、spatialDWLS、RCTD等算法具有更优异的去卷积性能，能够更准确地拆分空间转录组数据的细胞类型组成。
另外，来自2022年的Briefings in Bioinformatics综述文章（DOI: 10.1093/bib/bbac245）测试了10种整合算法：
“Taken together patterns observed from internal and external reference, our results suggest that RCTD, cell2location and stereoscope are among the most robust to batch effects between reference scRNA-seq and target ST data.”
同样表明RCTD和cell2location算法依旧"遥遥领先"！
整合算法很多，我们按照排名一个一个来，这次我们先介绍cell2location以及代码实战。
Cell2location发表于2022年的Nature Biotechnology，题为《Cell2location maps fine-grained cell types in spatial transcriptomics》。该方法依赖于贝叶斯模型，旨在解析空间转录组数据中更为精细的细胞类型的，并创建不同组织的细胞组成图谱。
本文将介绍：
Cell2location算法的工作流程； 使用Python代码实战。
一. Cell2location算法的工作流程
Cell2location采用贝叶斯层次框架。首先使用外部单细胞RNA测序数据作为参考数据，估计细胞类型特异性基因表达特征。然后，通过负二项分布（negative binomial regression）对观察到的空间表达计数矩阵建模，其中均值参数取决于参考细胞类型的特征，过度分散参数使用指数-伽马复合先验建模，旨在使大多数基因具有低过度分散。基因特异性技术灵敏度和基因-位置特异性的加性偏移被包括在均值参数的一部分，每个都使用单独的层次伽马先验进行建模。Cell2location进一步使用层次伽马先验对细胞类型特征的回归权重进行建模，并将回归权重分解为多个潜在组的贡献，这可以解释为具有共享细胞类型丰度特征的斑点，旨在跨具有相似细胞组成的位置借用强度。最后，cell2location采用变分贝叶斯推断来近似后验分布，并相应地生成参数估计。总的来说，Cell2location提供了一种可靠的途径来处理模型中的不确定因素，可以解释细胞类型丰度的线性依存关系。相比于其他的联合分析工具，Cell2location体现了更高的运算效率。
注意Cell2location运行非常慢，建议使用GPU进行加速！
二. Cell2location代码实现
Step0.环境部署
github在https://github.com/BayraktarLab/cell2location#Installation
官方教程在：https://cell2location.readthedocs.io/en/latest/
可以使用Conda安装Cell2location：
conda create -y -n cell2loc python=3.9
conda activate cell2loc
pip install cell2location[tutorials]
如果需要衔接Jupyter notebook使用的话，需要在cell2loc这个环境里安装一下几个插件，然后就能在Jupyter notebook里选择cell2loc这个环境：
mamba install -y nb_conda_kernels ipykernel
python -m ipykernel install --user --name cell2loc --display-name cell2loc
Step1.加载示例数据及预处理
示例数据我们使用上期的CellTrek算法的示例小鼠皮层数据【整合单细胞和空转数据多种方法之CellTrek】。运行流程图：
Python包加载：
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
注意：cell2location需要count原始数据作为input数据，不能用标准化后的。参考：https://github.com/BayraktarLab/cell2location/issues/249
1.1 加载示例空转数据：
adata_st = sc.read_h5ad('Step1.brain_st_cortex.h5ad')
adata_st
AnnData object with n_obs × n_vars = 1075 × 31053
obs: 'in_tissue', 'array_row', 'array_col', 'orig.ident', 'nCount_Spatial', 'nFeature_Spatial', 'slice', 'region', 'id', 'Spatial_snn_res.0.8', 'seurat_clusters', 'keep_idx'
var: 'vst.mean', 'vst.variance', 'vst.variance.expected', 'vst.variance.standardized', 'vst.variable'
uns: 'seurat_clusters_colors', 'spatial'
obsm: 'X_pca', 'X_umap', 'spatial'
检查数据：
sc.pl.spatial(adata_st, color='seurat_clusters')
线粒体编码基因与空间定位无关，因为它们的表达代表单细胞和细胞核数据中的技术产物，而不是线粒体的生物丰度。然而，这些基因在每个位置组成了15-40%的mRNA。因此，为了避免噪音，我们强烈建议去除线粒体基因：
mt-小鼠MT-人
# find mitochondria-encoded (MT) genes
adata_st.var['MT_gene'] = [gene.startswith('mt-') for gene in adata_st.var.index]
# remove MT genes for spatial mapping (keeping their counts in the object)
adata_st.obsm['MT'] = adata_st[:, adata_st.var['MT_gene'].values].X.toarray()
adata_st = adata_st[:, ~adata_st.var['MT_gene'].values]
adata_st
View of AnnData object with n_obs × n_vars = 1075 × 31040
obs: 'in_tissue', 'array_row', 'array_col', 'orig.ident', 'nCount_Spatial', 'nFeature_Spatial', 'slice', 'region', 'id', 'Spatial_snn_res.0.8', 'seurat_clusters', 'keep_idx'
var: 'vst.mean', 'vst.variance', 'vst.variance.expected', 'vst.variance.standardized', 'vst.variable', 'MT_gene'
uns: 'seurat_clusters_colors', 'spatial'
obsm: 'X_pca', 'X_umap', 'spatial', 'MT'
过滤掉一些线粒体基因。
1.2 加载示例单细胞数据：
adata_ref = sc.read_h5ad('Step1.brain_sc.h5ad')
adata_ref
AnnData object with n_obs × n_vars = 4785 × 34617
obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'sample_id', 'sample_type', 'organism', 'donor', 'sex', 'age_days', 'eye_condition', 'genotype', 'driver_lines', 'reporter_lines', 'brain_hemisphere', 'brain_region', 'brain_subregion', 'injection_label_direction', 'injection_primary', 'injection_secondary', 'injection_tract', 'injection_material', 'injection_exclusion_criterion', 'facs_date', 'facs_container', 'facs_sort_criteria', 'rna_amplification_set', 'library_prep_set', 'library_prep_avg_size_bp', 'seq_name', 'seq_tube', 'seq_batch', 'total_reads', 'percent_exon_reads', 'percent_intron_reads', 'percent_intergenic_reads', 'percent_rrna_reads', 'percent_mt_exon_reads', 'percent_reads_unique', 'percent_synth_reads', 'percent_ecoli_reads', 'percent_aligned_reads_total', 'complexity_cg', 'genes_detected_cpm_criterion', 'genes_detected_fpkm_criterion', 'tdt_cpm', 'gfp_cpm', 'class', 'subclass', 'cluster', 'confusion_score', 'cluster_correlation', 'core_intermediate_call', 'id', 'cell_type', 'dbscan_filter'
var: 'vst.mean', 'vst.variance', 'vst.variance.expected', 'vst.variance.standardized', 'vst.variable'
obsm: 'X_pca', 'X_umap'
检查数据：
sc.settings.set_figure_params(dpi=100, dpi_save=300, figsize=(4, 4))
sc.pl.umap(adata_ref, color='seurat_clusters')
注意：在估计参考细胞类型的signature之前，作者建议进行宽松一点的基因选择。作者更推荐于这种方式，而不是标准的高变异基因选择，因为作者的程序保留了罕见基因的标记，同时去除了大多数无信息的基因。
默认参数cell_count_cutoff=5，cell_percentage_cutoff2=0.03，nonz_mean_cutoff=1.12比较合适，但是用户可以增加截断值以排除更多的基因。为了保留罕见细胞类型的标记基因，作者建议将cell_count_cutoff设置为较低的值，例如5，但是cell_percentage_cutoff2和nonz_mean_cutoff可以适当增加，可以控制在8,000至16,000个基因。
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()
在这个二维直方图中，橙色矩形突出显示了基于表达该基因的细胞数量（Y轴）和检测到该基因的细胞的平均RNA计数（X轴）的组合而排除的基因。
Step2. 参考细胞类型特征估计(NB回归)
根据scRNA-seq数据估计特征，考虑批效应，使用负二项回归模型。这个步骤非常慢，最好有GPU加速：
首先，为回归模型准备一个数据对象:
# prepare anndata for the regression model
cell2location.models.RegressionModel.setup_anndata(adata=adata_ref,
batch_key='orig.ident', # 10X reaction / sample / batch
labels_key='cell_type' # cell type, covariate used for constructing signatures
#categorical_covariate_keys=['Method'] # multiplicative technical effects (platform, 3' vs 5', donor effect)
)
# create the regression model
from cell2location.models import RegressionModel
mod = RegressionModel(adata_ref)
# view anndata_setup as a sanity check
mod.view_anndata_setup()
现在我们训练模型以估计参考细胞类型的特征。
请注意，为了在您的数据上实现收敛（即获得损失的稳定），用户可能需要增加 max_epochs=250（请参见下文）。
还请注意，这里我们使用的 batch_size=2500 要比 scvi-tools 默认值大得多，并对所有数据中的所有细胞进行训练（train_size=1）- 这两个参数都是默认值。
mod.train(max_epochs=250, use_gpu=True,train_size=1)
确定模型是否需要更多训练。
在这里，我们绘制训练过程中的 ELBO 损失历史，从图中去掉前 20 个 epochs。这个图应该呈下降趋势，并在训练结束时趋于平稳。如果仍在下降，请增加 max_epochs。
mod.plot_history(20)
上面这幅图就提示我们需要增加
max_epochs了，官方教程的图是这样的，在训练结束时趋于平稳：
# 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, 'use_gpu': True}
)
# Save model
mod.save("./reference_signatures/", overwrite=True)
# Save anndata object with results
adata_ref.write("./sc_cell2location.h5ad")
用户可以直接计算后验分布的5%、50%和95%分位数，而不是使用1000个样本（或其他任何分位数）从分布中抽取。这样可以加快在大型数据集上的应用并减少内存需求，但无法用这种方式计算后验均值和标准差。
下面使用的是官方的代码报错了：
adata_ref = mod.export_posterior(
adata_ref, use_quantiles=True,
# choose quantiles
add_to_obsm=["q05","q50", "q95", "q0001"],
sample_kwargs={'batch_size': 2500, 'use_gpu': True}
)
我检查了这个函数：
add_to_obsm似乎被add_to_varm替换了，更新一下作者的代码，运行就成功了：
adata_ref = mod.export_posterior(
adata_ref, use_quantiles=True,
# choose quantiles
add_to_varm=["q05","q50", "q95", "q0001"],
sample_kwargs={'batch_size': 2500, 'use_gpu': True}
)
接下来的
mod.plot_QC也报错了，我这里跳过：
#报错了，跳过
#mod.plot_QC()
#模型和输出的h5ad，可以像这样加载:
adata_ref = sc.read_h5ad("./sc_cell2location.h5ad")
mod = cell2location.models.RegressionModel.load("./reference_signatures/", adata_ref)
提取参考细胞类型签名作为pd.DataFrame。
负二项回归模型的所有参数都被导出到参考 anndata 对象中，但是对于空间映射，我们只需要每个基因在每个细胞类型中的估计表达。在这里，我们从标准输出中提取这些信息：
# 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:5, 0:15]
Step3. Cell2location: spatial mapping
寻找共享基因，准备数据。基于参考signature，对空转数据取子集:
# find shared genes and subset both anndata and reference signatures
intersect = np.intersect1d(adata_st.var_names, inf_aver.index)
adata_st = adata_st[:, intersect].copy()
inf_aver = inf_aver.loc[intersect, :].copy()
# prepare anndata for cell2location model
cell2location.models.Cell2location.setup_anndata(adata=adata_st, batch_key="orig.ident")
重要提示： 要使用cell2location空间映射模型，用户需要指定两个超参数（N_cells_per_location和detection_alpha）。有关设置这些超参数及其影响的详细指导，请参阅流程图和说明。
这个cell2location好多重要提示。。。
选择超参数 N_cells_per_location！
将预期的每个组织位置的细胞丰度N_cells_per_location调整到实际组织情况是很有用的。这个值可以从成对的组织学图像中估计，如上面的说明所述。请将本教程中呈现的值（N_cells_per_location=30）更改为您的组织中观察到的值。
选择超参数 detection_alpha！
为了提高在幻灯片/批次内RNA检测敏感性存在较大技术变异性的数据集上的准确性和灵敏性 - 用户需要放宽对每个位置规范化的正则化（使用detection_alpha=20）。当用户观察到每个位置的总RNA计数的空间分布与基于组织学检查的预期细胞数不匹配时，说明在用户的样本中存在较大的RNA检测敏感性技术变异性。
作者最初选择高正则化（detection_alpha=200）作为默认值，因为作者在论文中使用的小鼠大脑和人类淋巴结数据集具有较低的技术效应，并且使用高正则化强度提高了每个位置的总估计细胞丰度与组织学定量核数之间的一致性（请参见cell2location论文的图S8F）。然而，在许多合作中，作者发现在人类组织的Visium实验中存在技术效应。这促使了将detection_alpha的新默认值设为20以及建议在用户的数据上测试这两个设置（detection_alpha=20和detection_alpha=200）。
# create and train the model
mod = cell2location.models.Cell2location(
adata_st, 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()
Training cell2location（这个步骤非常慢）:
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,
use_gpu=True,
)
# plot ELBO loss history during training, removing first 100 epochs from the plot
mod.plot_history(1000)
plt.legend(labels=['full data training']);
这里我遇到报错了:
ValueError: Expected value argument (Tensor of shape (2500, 20925)) to be within the support (IntegerGreaterThan(lower_bound=0)) of the distribution GammaPoisson(), but found invalid values:
tensor([[0.0000, 0.9725, 0.9725, ..., 0.0000, 0.0000, 0.0000],
[0.0000, 0.0000, 0.0000, ..., 0.0000, 0.0000, 0.0000],
[0.0000, 0.0000, 0.0000, ..., 0.0000, 0.0000, 0.0000],
...,
[0.0000, 0.0000, 0.0000, ..., 0.0000, 0.0000, 0.0000],
[0.0000, 0.0000, 0.0000, ..., 0.0000, 0.0000, 0.0000],
[0.0000, 0.0000, 0.0000, ..., 0.0000, 0.0000, 0.0000]],
device='cuda:0')
查了github：
https://github.com/BayraktarLab/cell2location/issues/128 https://github.com/BayraktarLab/cell2location/issues/120
原来是我的空转数据不是raw count数据，我重新把空转数据读入scanpy...这里跳过这个debug步骤，再次强调：cell2location的单细胞和空转需要raw count数据作为input！
另外，遇到报错不用紧张，多查查github或者是中英文帖子，如果一个算法不是很小众的话，基本会有答案的！
这里为了节省时间，我没有跑30000次迭代（在GPU的加速下需要运行50分钟），而是选择了300次，在GPU的加速下运行1分钟不到。
Step4. 导出预测的细胞丰度结果
# In this section, we export the estimated cell abundance (summary of the posterior distribution).
adata_st = mod.export_posterior(
adata_st, sample_kwargs={'num_samples': 1000, 'batch_size': mod.adata.n_obs, 'use_gpu': True}
)
# Save model
mod.save("./cell2location_map/", overwrite=True)
# Save anndata object with results
adata_st.write("./cell2location_map/st_cell2location_res.h5ad")
模型和输出h5ad，可以像这样重新加载：
adata_st = sc.read_h5ad("./cell2location_map/st_cell2location_res.h5ad")
mod = cell2location.models.Cell2location.load("./cell2location_map/", adata_st)
mod.plot_QC()
#fig = mod.plot_spatial_QC_across_batches()
图形呈对角线状说明质量较好。
导出预测的细胞丰度结果为csv文件：
作者建议使用后验分布的5%分位数的结果，表示模型具有高置信度的细胞丰度值(即“至少存在这个数量”)。结果在
adata_st.obsm下：
adata_st.obsm
#AxisArrays with keys: spatial, MT, means_cell_abundance_w_sf, stds_cell_abundance_w_sf, q05_cell_abundance_w_sf, q95_cell_abundance_w_sf
pd.DataFrame(adata_st.obsm['q05_cell_abundance_w_sf']).to_csv("./cell2location_map/st_cell2location_res.csv")
这个结果可能对Python党来说没啥用，作为R语言老用户，我觉得可以读到Seurat对象里进行可视化。
Step5. 可视化
终于到可视化环节了，
cell2location运行确实有点慢...
作者建议使用后验分布的5%分位数的结果，表示模型具有高置信度的细胞丰度值(即“至少存在这个数量”):
# 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']
adata_st.obsm['q05_cell_abundance_w_sf']
我们这里只有一张切片，就不用跑这个步骤了slide = select_slide(adata_vis, 'V1_Human_Lymph_Node')
# 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(adata_st, cmap='magma',
# show first 8 cell types
color=['L6 CT', 'L5 IT', 'L6b', 'L2/3 IT', 'L4', 'Astro'],
ncols=4, size=1.3,
img_key='lowres',
# limit color scale at 99.2% quantile of cell abundance
vmin=0, vmax='p99.2'
)
我不是很懂皮层，可以参考了Seurat官网的结果，可以看到cell2location预测还是比较准确的：
# 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 = ['L6 CT', 'L5 IT', 'L6b']
clust_col = ['' + str(i) for i in clust_labels] # in case column names differ from labels
#slide = select_slide(adata_st, 'V1_Human_Lymph_Node')
with mpl.rc_context({'figure.figsize': (10, 10)}):
fig = plot_spatial(
adata=adata_st,
# labels to show on a plot
color=clust_col, labels=clust_labels,
show_img=True,
img_key='lowres',
# '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'
)
# 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 = ['L6 CT', 'L5 IT', 'L6b']
clust_col = ['' + str(i) for i in clust_labels] # in case column names differ from labels
#slide = select_slide(adata_st, 'V1_Human_Lymph_Node')
with mpl.rc_context({'figure.figsize': (10, 10)}):
fig = plot_spatial(
adata=adata_st,
# labels to show on a plot
color=clust_col, labels=clust_labels,
show_img=False,
img_key='lowres',
# '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'
)
Step6. 下游分析
6.1 用Leiden聚类识别离散组织区域
通过使用由cell2location估计的细胞丰度对位置进行聚类，我们识别在细胞组成方面存在差异的组织区域。
通过使用每种细胞类型的估计细胞丰度对Visium点进行聚类，我们找到组织区域。我们构建一个表示在估计的细胞丰度中位置相似性的K最近邻（KNN）图，然后应用Leiden聚类。KNN邻居的数量应根据数据集的大小和解剖学定义的区域大小进行调整（例如，与大脑的大小相比，海马区域相对较小，因此可能被大量邻居掩盖）。可以在一系列的KNN邻居和Leiden聚类分辨率下执行此操作，直到获得与组织解剖结构匹配的聚类。
聚类是在所有Visium截面/批次中共同完成的，因此区域身份是直接可比的。当存在多个批次之间的强技术效应时（在本例中并非如此），原则上可以使用sc.external.pp.bbknn来在KNN构建过程中考虑这些效应。
下图的seurat_clusters是来自seurat的无监督聚类结果，可以看到除了分辨率设置高低的区别以外，scanpy与Seurat的亚群大致是对应的（包含关系）。
# compute UMAP using KNN graph based on the cell2location output
sc.tl.umap(adata_st, min_dist = 0.3, spread = 1)
# show regions in UMAP coordinates
with mpl.rc_context({'axes.facecolor': 'white',
'figure.figsize': [5, 5]}):
sc.pl.umap(adata_st, color=['seurat_clusters','region_cluster'], size=30,
color_map = 'RdPu', ncols = 2, legend_loc='on data',
legend_fontsize=20)
# plot in spatial coordinates
with mpl.rc_context({'axes.facecolor': 'black',
'figure.figsize': [4.5, 5]}):
sc.pl.spatial(adata_st, color=['seurat_clusters','region_cluster'],
size=1.3, img_key='lowres', alpha=0.5)
6.2 使用矩阵分解(NMF)识别细胞区室/组织区
在这里，我们使用cell2location映射结果来识别细胞类型的空间共存，以更好地理解组织结构并预测细胞间的相互作用。我们对来自cell2location的细胞类型丰度估计进行了非负矩阵分解（NMF）。类似于将NMF应用于常规scRNA-seq的已建立的优势，加性NMF分解将空间细胞类型丰度文件组成组件，捕获共定位的细胞类型。这种基于NMF的分解自然地考虑了多个细胞类型和微环境可以在同一Visium位置共存的事实，同时在组织区域之间共享信息。
提示：
在实践中，最好对一系列因子R=5,..,30进行NMF训练，并选择合适的R。
如果要找到一些最明显的细胞区，使用较少的因子。如果要找到非常强的共定位信号并假设大多数细胞类型不共定位，则使用大量因子（>30）。
下面作者展示了如何执行这个分析。作者将NMF包装成一个pipeline：
from cell2location import run_colocation
res_dict, adata_st = run_colocation(
adata_st,
model_name='CoLocatedGroupsSklearnNMF',
train_args={
'n_fact': np.arange(11, 13), # IMPORTANT: use a wider range of the number of factors (5-30)
'sample_name_col': 'orig.ident', # columns in adata_vis.obs that identifies sample
'n_restarts': 3 # number of training restarts
},
# the hyperparameters of NMF can be also adjusted:
model_kwargs={'alpha': 0.01, 'init': 'random', "nmf_kwd_args": {"tol": 0.000001}},
export_args={'path': './CoLocatedComb/'}
)
对于每个因子编号，模型都会生成以下文件夹输出列表：
cell_type_fractions_heatmap/”：细胞类型（行）在NMF组件（列）上估计的NMF权重的点图；cell_type_fractions_mean/：用于点图的数据factor_markers/：列出每个NMF因子最具特异性的前10个细胞类型的表格models/：保存的NMF模型predictive_accuracy/：2D直方图绘图，显示NMF如何解释cell2location输出spatial/：空间坐标中位置上的NMF权重location_factors_mean/：用于空间坐标图的数据stability_plots/：在训练重新启动之间NMF权重的稳定性
提示：
NMF模型的输出，如因子加载，存储在adata.uns[f"mod_coloc_n_fact{n_fact}"]中，格式与adata.uns['mod']中的主要cell2location结果类似。
# Here we plot the NMF weights (Same as saved to `cell_type_fractions_heatmap`)
res_dict['n_fact12']['mod'].plot_cell_type_loadings()
6.3 估计空间数据中每个基因的细胞类型特异性表达
在空间坐标中绘制基因的细胞类型特异性表达。
# Compute expected expression per cell type
expected_dict = mod.module.model.compute_expected_per_cell_type(
mod.samples["post_sample_q05"], mod.adata_manager
)
# Add to anndata layers
for i, n in enumerate(mod.factor_names_):
adata_st.layers[n] = expected_dict['mu'][i]
# Save anndata object with results
#adata_st.write("./sp.h5ad")
遇到报错了，解决方案：
https://github.com/BayraktarLab/cell2location/issues/302 https://github.com/BayraktarLab/cell2location/issues/143
运行上述代码前需要先运行：
mod.samples = adata_st.uns['mod']
然后运行：
# Compute expected expression per cell type
expected_dict = mod.module.model.compute_expected_per_cell_type(
mod.samples["post_sample_q05"], mod.adata_manager
)
# Add to anndata layers
for i, n in enumerate(mod.factor_names_):
adata_st.layers[n] = expected_dict['mu'][i]
# Save anndata object with results
#adata_st.write("./sp.h5ad")
下面的可视化部分也有bug，主要是plot_genes_per_cell_type这个绘图函数没有封装进cell2location 【ModuleNotFoundError: No module named 'tutorial_utils'】，解决方案https://github.com/BayraktarLab/cell2location/issues/280；
原因是plot_genes_per_cell_type这个函数适合用来运行作者的示例数据，但是不适合我们自己的数据，所以作者没把这个函数封装进python包。
这里我对plot_genes_per_cell_type函数改造了一下：
def plot_genes_per_cell_type(slide, genes, ctypes):
slide.var['SYMBOL'] = slide.var.index
n_genes = len(genes)
n_ctypes = len(ctypes)
fig, axs = plt.subplots(
nrows=n_genes, ncols=n_ctypes + 1, figsize=(4.5 * (n_ctypes + 1) + 2, 5 * n_genes + 1), squeeze=False
)
# axs = axs.reshape((n_genes, n_ctypes+1))
# plots of every gene
for j in range(n_genes):
# limit color scale at 99.2% quantile of gene expression (computed across cell types)
quantile_across_ct = np.array(
[
np.quantile(slide.layers[n][:, slide.var["SYMBOL"] == genes[j]].toarray(), 0.992)
for n in slide.uns["mod"]["factor_names"]
]
)
quantile_across_ct = np.partition(quantile_across_ct.flatten(), -2)[-2]
sc.pl.spatial(
slide,
cmap="magma",
color=genes[j],
# layer=ctypes[i],
#gene_symbols="SYMBOL",
ncols=4,
size=1.3,
#img_key="hires",
# limit color scale at 99.2% quantile of gene expression
vmin=0,
vmax="p99.2",
ax=axs[j, 0],
show=False,
)
# plots of every cell type
for i in range(n_ctypes):
sc.pl.spatial(
slide,
cmap="magma",
color=genes[j],
layer=ctypes[i],
gene_symbols="SYMBOL",
ncols=4,
size=1.3,
#img_key="hires",
# limit color scale at 99.2% quantile of gene expression
vmin=0,
vmax=quantile_across_ct,
ax=axs[j, i + 1],
show=False,
)
axs[j, i + 1].set_title(f"{genes[j]} {ctypes[i]}")
return fig, axs
然后运行就成功了：
# list cell types and genes for plotting
ctypes = ['L6 CT', 'L5 IT', 'L6b']
genes = ["Hpca", "Plp1"]
with mpl.rc_context({'axes.facecolor': 'black'}):
# select one slide
# slide = select_slide(adata_vis, 'V1_Human_Lymph_Node')
plot_genes_per_cell_type(adata_st, genes,ctypes);
我不太懂皮层，所以随便找了几个基因做可视化展示。
除了这些可视化内容，cell2location还有一些个性化绘图内容，有机会我们再做介绍！
- END -
https://mp.weixin.qq.com/s/Y1rBQh-G7R_EaEdNMENHdA