新鲜出炉！顶刊单细胞算法 Mellon，从单个细胞层面估计细胞状态密度

最新推荐文章于 2024-09-13 16:27:15 发布

生信碱移

最新推荐文章于 2024-09-13 16:27:15 发布

阅读量641

点赞数 15

文章标签：算法人工智能前端 python 大数据数据可视化

本文链接：https://blog.csdn.net/weixin_49214410/article/details/141218202

版权

生信碱移

Mellon估计细胞密度

分化是一个动态过程，是所有多细胞生物发育和功能的基础。细胞状态密度是细胞分布的一种表示，受基本生物学过程的影响，体现着细胞分化过程中的关键转变。一般来讲，增殖可以增加处于某种状态的细胞数量，从而导致高细胞状态密度，而细胞凋亡会降低细胞状态密度（图1b）。当细胞变化不大时，它们可以汇聚到检查点，从而导致高细胞状态密度。相反，如过在罕见的瞬时细胞中（比如分化的中间阶段）所见的转录加速会导致降低细胞状态密度。简而言之，快速、协调的转录加速将导致低密度的瞬时状态区域，后者将连接高密度细胞区域，这已被证明是从人类到植物的不同环境中发育进程的标志。

▲ 图1：正常的细胞密度（图a）与实际生理分化过程中的细胞密度（图b）。

单细胞研究强调了细胞状态密度异质性的重要性。尽管如此，目前在单细胞数据中进行高维密度估计的方法通常会产生嘈杂的结果，并且难以提供具有生物学意义的解释。为此，来自福瑞德·哈金森癌症研究中心的研究者开发了一款单细胞算法工具 Mellon，旨在从单细胞数据中稳健地估计细胞状态密度，该算法于 2024 年 6 月 18 日发表于 Nature Methods [IF: 48.0] 杂志。

▲ DOI: 10.1038/s41592-024-02302-w

从原理上讲，① Mellon 首先从单细胞数据出发，使用扩散图（diffusion maps, 默认）或其他降维方法（如UMAP）将数据映射到一个高维表示空间。② 随后，将进一步计算每个细胞的最近邻距离，并利用这些距离推测局部密度。③ 随后，Mellon 通过构建高斯过程模型，将相似细胞状态之间的密度进行平滑连接和估计，推断出整个高维细胞状态空间的连续密度函数。④ 最后再使用推断的密度函数，对每个细胞的状态密度进行估计。⑤ 不仅如此，作者还在上述基础上设计了一种基因变化分析程序，以确定在低密度区域中表达变化显著的基因。根据基因在某状态下相对于其邻居细胞表达的绝对差异，对在低密度和相邻高密度区域的细胞中对基因进行排名，并用密度的倒数加权局部变异性。这使得在低密度区域中表达变化较大的基因排名会更高（换句话说，也就是对基因能够给出一个重要性的评估，这个评估是基于细胞发育转换过程中基因的短期变化）。

▲ 原理概要：Mellon 依赖于细胞邻域距离和密度之间的内在关系，假设最近邻距离的分布与细胞状态密度相关。

由于连续密度函数的构建，使得 Mellon 能够捕捉细胞状态空间的连续特征，并且从不同类型的单细胞多模态数据中推断细胞状态密度。

Mellon 算法整体通过 Python 实现，基于scanpy库单细胞对象实现了简便的细胞密度估计。下面简单介绍一下它的使用，更多详细信息可以看看 github 仓：

https://github.com/settylab/Mellon

1、加载 Mellon 与相关依赖库

用户可以通过以下代码加载Mellon与相关依赖库：

import numpy as np

import matplotlib
import matplotlib.pyplot as plt

from sklearn.cluster import k_means
import palantir
import mellon
import scanpy as sc

import warnings
from numba.core.errors import NumbaDeprecationWarning

# 忽略特定类型的警告消息
warnings.simplefilter("ignore", category=NumbaDeprecationWarning)

绘图偏好的设置：

# 绘图偏好的设置
%matplotlib inline
matplotlib.rcParams["figure.figsize"] = [4, 4]
matplotlib.rcParams["figure.dpi"] = 125
matplotlib.rcParams["image.cmap"] = "Spectral_r"
# no bounding boxes or axis:
matplotlib.rcParams["axes.spines.bottom"] = "on"
matplotlib.rcParams["axes.spines.top"] = "off"
matplotlib.rcParams["axes.spines.left"] = "on"
matplotlib.rcParams["axes.spines.right"] = "off"

2、示例数据加载

演示的 scRNA-seq 数据集是 Mellon 手稿中所述的预处理 T 细胞耗竭的骨髓数据ad，为Scanpy库的 Anndata 对象格式（已经包含了常规分析后的数据内容）：

ad_url = "https://fh-pi-setty-m-eco-public.s3.amazonaws.com/mellon-tutorial/preprocessed_t-cell-depleted-bm-rna.h5ad"
ad = sc.read("data/preprocessed_t-cell-depleted-bm-rna.h5ad", backup_url=ad_url)
ad
#AnnData object with n_obs × n_vars = 8627 × 17226
#   obs: 'sample', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'batch', 'DoubletScores', 'n_counts', 'leiden', 'phenograph', 'log_n_counts', 'celltype', 'palantir_pseudotime', 'selection', 'NaiveB_lineage', 'mellon_log_density', 'mellon_log_density_clipped'
#    var: 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'PeakCounts'
#    uns: 'DMEigenValues', 'DM_EigenValues', 'NaiveB_lineage_colors', 'celltype_colors', 'custom_branch_mask_columns', 'hvg', 'leiden', 'mellon_log_density_predictor', 'neighbors', 'pca', 'sample_colors', 'umap'
#    obsm: 'DM_EigenVectors', 'X_FDL', 'X_pca', 'X_umap', 'branch_masks', 'chromVAR_deviations', 'palantir_branch_probs', 'palantir_fate_probabilities', 'palantir_lineage_cells'
#   varm: 'PCs', 'geneXTF'
#   layers: 'Bcells_lineage_specific', 'Bcells_primed', 'MAGIC_imputed_data'
#    obsp: 'DM_Kernel', 'DM_Similarity', 'connectivities', 'distances', 'knn'

可视化一下：

sc.pl.scatter(ad, basis="umap", color="celltype")

3、细胞密度计算

①预处理：

%%time
dm_res = palantir.utils.run_diffusion_maps(ad, pca_key="X_pca", n_components=20)
#Determing nearest neighbor graph...
#CPU times: user 1min 35s, sys: 4.27 s, total: 1min 39s
#Wall time: 1min 30s

② 细胞密度计算：

%%time
# 细胞密度计算
model = mellon.DensityEstimator()
log_density = model.fit_predict(ad.obsm["DM_EigenVectors"])

predictor = model.predict
# 更新数据对象
ad.obs["mellon_log_density"] = log_density
ad.obs["mellon_log_density_clipped"] = np.clip(
    log_density, *np.quantile(log_density, [0.05, 1])
)
#[2023-07-09 08:16:44,570] [INFO    ] Computing nearest neighbor distances.
#[2023-07-09 08:16:45,799] [INFO    ] Using covariance function Matern52(ls=0.007897131365355458).
#[2023-07-09 08:16:45,802] [INFO    ] Computing 5,000 landmarks with k-means clustering.
#[2023-07-09 08:16:48,405] [INFO    ] Doing low-rank Cholesky decomposition for 8,627 samples and 5,000 landmarks.
#[2023-07-09 08:16:53,261] [INFO    ] Using rank 5,000 covariance representation.
#[2023-07-09 08:16:54,303] [INFO    ] Running inference using L-BFGS-B.
#[2023-07-09 08:17:11,188] [INFO    ] Computing predictive function.
#CPU times: user 49.5 s, sys: 48.5 s, total: 1min 37s
#Wall time: 28.4 s

# 更新的anndata对象增加了 obs['log_density'] 和 obs['log_density_clipped'] 两列
ad
#AnnData object with n_obs × n_vars = 8627 × 17226
#    obs: 'sample', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'batch', 'DoubletScores', 'n_counts', 'leiden', 'phenograph', 'log_n_counts', 'celltype', 'palantir_pseudotime', 'selection', 'NaiveB_lineage', 'mellon_log_density', 'mellon_log_density_clipped'
#    var: 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'PeakCounts'
#    uns: 'DMEigenValues', 'DM_EigenValues', 'NaiveB_lineage_colors', 'celltype_colors', 'custom_branch_mask_columns', 'hvg', 'leiden', 'mellon_log_density_predictor', 'neighbors', 'pca', 'sample_colors', 'umap'
#    obsm: 'DM_EigenVectors', 'X_FDL', 'X_pca', 'X_umap', 'branch_masks', 'chromVAR_deviations', 'palantir_branch_probs', 'palantir_fate_probabilities', 'palantir_lineage_cells'
#    varm: 'PCs', 'geneXTF'
#    layers: 'Bcells_lineage_specific', 'Bcells_primed', 'MAGIC_imputed_data'
#    obsp: 'DM_Kernel', 'DM_Similarity', 'connectivities', 'distances', 'knn'

可视化一下密度分布：

sc.pl.scatter(
    ad, color=["mellon_log_density", "mellon_log_density_clipped"], basis="umap"
)

4、针对结果的进一步分析

① 展示不同细胞类型的密度估计：

fig, (ax1, ax2) = plt.subplots(1, 2, width_ratios=[3, 2], figsize=[12, 4])
sc.pl.violin(ad, "mellon_log_density", "celltype", rotation=45, ax=ax1, show=False)
sc.pl.scatter(ad, color="celltype", basis="umap", ax=ax2, show=False)
plt.show()

② 比较细胞的估计密度与伪时间分析结果。注意，示例数据展示的是 palantir 算法的伪时间估计，这在数据中体现在细胞的伪时间估计ad.obs['palantir_pseudotime] 以及细胞所属的发育轨迹分支 ad.obsm['palantir_lineage_cells']（展示如下）：

ad.obs["palantir_pseudotime"], ad.obsm["palantir_lineage_cells"]
#(IM-1393_BoneMarrow_TcellDep_1_multiome#GTGAGCGAGTCTCACC-1    0.138621
# IM-1393_BoneMarrow_TcellDep_1_multiome#GAGTCAAAGTCCTTCA-1    0.336803
# IM-1393_BoneMarrow_TcellDep_1_multiome#TGTGCGCAGTCGCTAG-1    0.684445
# IM-1393_BoneMarrow_TcellDep_1_multiome#ATATGTCCAATGCCTA-1    0.311772
# IM-1393_BoneMarrow_TcellDep_1_multiome#CTTAGTTTCGCTAGTG-1    0.704443
#                                                                ...
# IM-1393_BoneMarrow_TcellDep_2_multiome#ATCCGTGAGGGATTAG-1    0.309793
# IM-1393_BoneMarrow_TcellDep_2_multiome#ATTTAGGTCAGGTTTA-1    0.306938
# IM-1393_BoneMarrow_TcellDep_2_multiome#GATGGACAGATAAAGC-1    0.311873
# IM-1393_BoneMarrow_TcellDep_2_multiome#GAAGGCTAGCTATATG-1    0.312016
# IM-1393_BoneMarrow_TcellDep_2_multiome#AGACAATAGGCTCATG-1    0.311489
# Name: palantir_pseudotime, Length: 8627, dtype: float64,
#                                                     NaiveB    Ery    pDC  \
# IM-1393_BoneMarrow_TcellDep_1_multiome#GTGAGCGA...   False  False  False
# IM-1393_BoneMarrow_TcellDep_1_multiome#GAGTCAAA...   False   True  False
# IM-1393_BoneMarrow_TcellDep_1_multiome#TGTGCGCA...   False   True  False
# IM-1393_BoneMarrow_TcellDep_1_multiome#ATATGTCC...   False  False  False
# IM-1393_BoneMarrow_TcellDep_1_multiome#CTTAGTTT...   False   True  False
# ...                                                    ...    ...    ...
# IM-1393_BoneMarrow_TcellDep_2_multiome#ATCCGTGA...   False  False  False
# IM-1393_BoneMarrow_TcellDep_2_multiome#ATTTAGGT...   False  False  False
# IM-1393_BoneMarrow_TcellDep_2_multiome#GATGGACA...   False  False  False
# IM-1393_BoneMarrow_TcellDep_2_multiome#GAAGGCTA...   False  False  False
# IM-1393_BoneMarrow_TcellDep_2_multiome#AGACAATA...   False  False  False
#
#                                                      Mono
# IM-1393_BoneMarrow_TcellDep_1_multiome#GTGAGCGA...   True
# IM-1393_BoneMarrow_TcellDep_1_multiome#GAGTCAAA...  False
# IM-1393_BoneMarrow_TcellDep_1_multiome#TGTGCGCA...  False
# IM-1393_BoneMarrow_TcellDep_1_multiome#ATATGTCC...   True
# IM-1393_BoneMarrow_TcellDep_1_multiome#CTTAGTTT...  False
# ...                                                   ...
# IM-1393_BoneMarrow_TcellDep_2_multiome#ATCCGTGA...   True
# IM-1393_BoneMarrow_TcellDep_2_multiome#ATTTAGGT...   True
# IM-1393_BoneMarrow_TcellDep_2_multiome#GATGGACA...   True
# IM-1393_BoneMarrow_TcellDep_2_multiome#GAAGGCTA...   True
# IM-1393_BoneMarrow_TcellDep_2_multiome#AGACAATA...   True

# [8627 rows x 4 columns])

看完了这个伪时间分析的数据格式，给大伙演示一下如何基于palantir.plot.plot_branch方法可视化细胞密度与伪时间估计得分：

palantir.plot.plot_branch(
    ad,
    branch_name="NaiveB", # 只看NaiveB这个分支轨迹
    position="mellon_log_density", # mellon估计的细胞密度
    color="celltype",  # 将颜色指定为细胞类型
    masks_key="palantir_lineage_cells", 
    s=100,
)
plt.show()

当然，还可以显示指定基因的表达水平：

# Plot imputed expressino of EBF1 in the comparison between pseudotime and log density
palantir.plot.plot_branch(
    ad,
    branch_name="NaiveB",
    position="mellon_log_density",
    color="EBF1",# 将颜色指定基因表达量
    color_layer="MAGIC_imputed_data",
    masks_key="palantir_lineage_cells",
    s=100,
)
plt.show()

③ 可视化指定基因表达的局部变异性，局部变异性通过比较指定基因在细胞与其相邻细胞的状态，提供了每个细胞状态下基因表达变化的量度。用户可以使用palantir.utils.run_local_variability方法计算基因表达的局部变异性，并进行可视化分析：

# 计算基因表达的局部变异性
palantir.utils.run_local_variability(ad)
# 可视化
palantir.plot.plot_branch(
    ad,
    branch_name="NaiveB",
    position="mellon_log_density",
    color="EBF1",
    color_layer="local_variability",
    masks_key="palantir_lineage_cells",
    s=100,
)
plt.show()

作者在原文对这个结果进行了解读，给各位老铁参考参考：

EBF1作为B细胞分化的主要调控因子，具有最高的变异分数。同时，其上调阶段精准定位于干细胞和分化细胞之间的低密度区域。

④ 指定细胞的子集，进一步通过 UMAPs 降维图探索细胞密度变化：

# 提取NaiveB轨迹下的细胞
bcell_lineage_cells = ad.obs_names[ad.obsm["palantir_lineage_cells"]["NaiveB"]]
# 通过bcell_lineage_cells指定细胞子集
sc.pl.embedding(
    ad[bcell_lineage_cells],
    basis="umap",
    color=["mellon_log_density_clipped", "celltype"],
)