细胞状态密度分布是表征表型景观的关键指标，对揭示分化、再生和疾病机制具有重要意义。Mellon是2024年发表在Nature Methods，通过单细胞或者空间单细胞转录组对细胞密度进行计算的一种工具。本文提供R版本。

1. 格式转换

library(zellkonverter)
load('PATH/scdata1.Rdata') #V4版本
writeH5AD(as.SingleCellExperiment(scdata1), "scdata1_merge.h5ad")

2.Mellon安装

pip install mellon 或者
conda install -c conda-forge mellon

依赖包：
scanpy
palantir:

pip install git+https://github.com/dpeerlab/Palantir

3.加载包

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)

4.matplotlib设置

#setting plot preferences
#%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"

5.读取数据

ad = sc.read("/PATH/scdata1_merge.h5ad") 

ad.obsm["X_umap"] = np.array(ad.obsm["UMAP"])  #for seurat V4(V5 not confirmed)

#umap plot (validation)
sc.pl.scatter(ad, basis="umap", color="merge_celltype", show=False)
plt.savefig("./1.celltype_umap.pdf", bbox_inches='tight') 
plt.close()

6. 数据预处理

ad.obsm["X_pca"] = ad.obsm["PCA"]
dm_res = palantir.utils.run_diffusion_maps(ad, pca_key="X_pca", n_components=20)

ad.obsm['X_diffmap'] = dm_res['EigenVectors']
ad.uns['diffmap_evals'] = dm_res['EigenValues']

ms_data = palantir.utils.determine_multiscale_space(ad)
ad.obsm['DM_EigenVectors_multiscaled'] = ms_data

7.密度计算

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])
)

可视化

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

plt.savefig("./2.umap_density.pdf", bbox_inches='tight') 
plt.close()

##plot3
sc.pl.violin(ad, "mellon_log_density", "merge_celltype", rotation=45,show=False)
plt.tight_layout()
plt.savefig("./3.celltype_density.pdf", bbox_inches='tight') 
plt.close()

8.亚群拟时序

#提取亚群
selected_cells = ad.obs['merge_celltype'].isin(['MC0', 'MC1', 'MC2', 'MC3', 'MC4', 'MC5'])
ad_subset = ad[selected_cells, :].copy()

##设置start and terminal cell
start_cell = ad_subset.obs[ad_subset.obs["nmf_cluster"] == "3"].index[0]  ##一个cell
terminal_cells = ad_subset.obs[ad_subset.obs["nmf_cluster"] == "MC2"].index.tolist() ##optional

pr_res = palantir.core.run_palantir(
    ad_subset,early_cell=start_cell, terminal_states=terminal_cells)

ad_subset.obs['palantir_pseudotime'] = pr_res.pseudotime
ad_subset.obsm['palantir_fate_probabilities'] = pr_res.branch_probs

fate_probs = ad_subset.obsm["palantir_fate_probabilities"]

##plot 
import pandas as pd
import matplotlib.pyplot as plt

ct_colors = pd.Series( 
    {
        '1': '#E686C1',  # 红色
        '2': '#7B4B33',  # 绿色
        '3': '#F28C26',  # 蓝色
        '4': '#2E7AC5',  # 黄色
        '5': '#A6322E',  # 品红
#        'MC5': '#DAD761',  # 青色
    },
    index=ad_subset.obs["nmf_cluster"].values.categories
)

plt.figure(figsize=[8, 3])
plt.scatter(
    ad_subset.obs["palantir_pseudotime"],
    ad_subset.obs["mellon_log_density"],
    s=10,
    color=ct_colors[ad_subset.obs["nmf_cluster"]],
)

unique_celltypes = np.unique(ad_subset.obs["nmf_cluster"])
handles = [
    plt.Line2D([0], [0], marker='o', color='w', markerfacecolor=ct_colors[ct], markersize=8)
    for ct in unique_celltypes
]
plt.legend(
    handles,
    unique_celltypes,
    title="Cell Type",
    bbox_to_anchor=(1.05, 1),
    loc='upper left',
    frameon=False
)

#  坐标轴标签
plt.xlabel("Pseudotime")
plt.ylabel("Mellon Log Density")

# 保存输出
plt.tight_layout()
plt.savefig("4.melan_time.pdf", bbox_inches='tight', dpi=300)
plt.close()

与monocle3类似，可以查看单个基因随时序的变化，只不过这里又增加了density维度（即，可以查看gene expression与cell density的关系）。

##GAS5
import pandas as pd
import numpy as np

# 创建一个只有一个 lineage（如 'Melan'）的布尔值 DataFrame
palantir_lineage_df = pd.DataFrame(
    {'Melan': np.ones(ad_subset.n_obs, dtype=bool)},
    index=ad_subset.obs_names
)
ad_subset.obsm['palantir_lineage_cells'] = palantir_lineage_df


plt.figure(figsize=[8, 3])

palantir.plot.plot_branch(
    ad_subset,
    branch_name="Melan",
    position="mellon_log_density",
    color="GAS5",
    color_layer="local_variability",
    masks_key="palantir_lineage_cells",
    s=100,
)

plt.savefig("4.melan_time_GAS5.pdf", bbox_inches='tight', dpi=300)
plt.close()

: