为什么要进行伪时序分析？ 在常规单细胞 RNA-seq 实验中，每个细胞通常只在一个时间点被测序。测序本身会破坏细胞，因此不能像活细胞成像那样直接追踪同一个细胞从早期状态到晚期状态的真实轨迹。
但是，很多生物过程并不是同步发生的。例如造血分化、胚胎发育、上皮细胞成熟、免疫细胞激活、肿瘤细胞状态转变等，同一个样本中可能同时捕获到处在不同进程阶段的细胞。伪时序分析（pseudotemporal ordering）的目的，就是利用这些 snapshot 数据中的转录组连续变化，把细胞按照某个生物过程的相对进展顺序排列起来。
简单理解就是：
离散捕获的单细胞状态 → 推断连续变化过程 → 给每个细胞一个相对顺序 这个相对顺序称为 pseudotime。在 root cell 选择合理、轨迹结构成立的前提下，pseudotime 小的细胞通常被解释为更接近起始状态，pseudotime 大的细胞通常被解释为更接近终末状态。
需要注意的是，伪时序不是实际时间。它的单位不是小时、天或发育阶段，而是细胞在表达空间中沿着某条推断轨迹前进的相对位置。因此，伪时序适合回答：
某些 marker gene 是否沿着状态转换逐渐上升或下降。 不同样本或条件中的细胞是否分布在不同的伪时序阶段。 它不适合直接回答：
任意两个 UMAP 上相邻的 cluster 是否存在发育关系。
伪时序分析适合什么数据？ 伪时序分析依赖一个关键假设：数据中确实存在一个主要的、连续的生物过程，并且这个过程能够在转录组空间中被捕捉到。
适合做伪时序的常见场景包括：
组织成熟过程，例如上皮细胞从 progenitor 到 mature cell。 免疫激活或耗竭过程，例如 naive T cell 到 effector/exhausted T cell。 肿瘤细胞状态连续变化，例如 EMT、增殖、缺氧或治疗反应状态。 不太适合直接做伪时序的情况包括：
样本中主要是彼此无关的成熟细胞类型，例如普通 PBMC 中把 T cells、B cells、monocytes、NK cells 全部放在一起强行排序。 UMAP 上看起来有一条线，但这条线主要由 total_counts、pct_counts_mt、细胞周期或批次效应驱动。 没有明确起始状态，也没有 marker gene、采样时间点或实验设计支持根细胞选择。 目标过程本身是环形或周期性，例如细胞周期，这类过程不应该简单用单根线性 pseudotime 解释。 因此，伪时序分析通常不应该在整个对象上直接运行，而是应该先基于细胞注释、marker gene 和研究问题，选出一个生物学上相关的 lineage 或 cell state subset。
伪时序的基本思想 不同伪时序方法的实现方式不同，但通常包含几个共同步骤：
选择细胞群 -> 构建低维表示或邻居图 -> 推断轨迹拓扑 -> 指定起始细胞或起始状态 -> 计算每个细胞到起点的相对距离 -> 解释基因表达和细胞状态随 pseudotime 的变化 其中最关键的是两个问题。
第一个问题是轨迹拓扑。也就是细胞状态之间应该被解释成线性、分叉、树状、环状，还是更复杂的图结构。不同工具能处理的拓扑不同。例如 DPT 适合基于扩散距离计算相对顺序，PAGA 更适合先观察 cluster 之间的全局连接结构，Monocle3 会学习 principal graph，Slingshot 则通过 cluster 和 principal curves 推断 lineage。
第二个问题是根细胞。多数伪时序方法需要指定起始细胞或起始状态。这个起点不能完全靠算法自动猜，最好来自已知 marker、时间点、实验处理、细胞类型层级或其他独立证据。如果根细胞选错，pseudotime 方向可能整体反过来，甚至得到完全不合理的排序。
常见方法怎么选择？ 当前常见方法可以粗略分成几类。
PAGA PAGA（partition-based graph abstraction）不是单纯的 pseudotime 计算工具，而是用来抽象 cluster 之间的全局拓扑结构。它会在细胞邻居图的基础上，估计不同 cluster 之间是否存在可靠连接。
PAGA 适合在正式计算伪时序前回答：
PAGA 可以和 DPT 搭配使用：先用 PAGA 检查拓扑，再基于合理的 root cell 计算 DPT pseudotime。
DPT DPT（diffusion pseudotime）基于扩散图和随机游走距离，从指定 root cell 出发给每个细胞计算相对伪时序。在 Scanpy 中可以直接使用 sc.tl.diffmap() 和 sc.tl.dpt() 。
DPT 的优点是容易接入 Scanpy 流程，适合作为 Python 工作流中的基础伪时序方法。它的限制是对 root cell、邻居图和连通性比较敏感；如果数据存在多个分离分支或 disconnected graph，需要谨慎解释。
在 Scanpy 中，DPT 更适合用于计算 pseudotime；分支结构和全局拓扑通常更建议结合 PAGA 判断，而不是单独依赖 DPT 的分支检测参数。
Monocle3 Monocle3 是 R 生态中常用的轨迹推断工具。它会在低维空间中学习 principal graph，并沿图结构计算 pseudotime。Monocle3 对分叉结构和轨迹可视化支持较好，常用于发育分化问题。
需要注意的是，Monocle3 中每个 partition 可以形成独立轨迹。如果某个 partition 没有指定 root node，相关细胞可能得到 infinite pseudotime。因此对多分支或多 partition 数据，需要明确每条轨迹的起点。
Slingshot Slingshot 也是 R 生态中成熟的轨迹推断方法。它通常基于低维坐标和 cluster label，先构建 cluster 间的 minimum spanning tree，再拟合 principal curves，最后给每条 lineage 上的细胞分配 pseudotime。
Slingshot 适合已经有较好聚类结果，并且希望在一个或多个 lineage 上平滑拟合轨迹的场景。它允许加入起始 cluster 或终末 cluster 作为先验信息。
读取已经注释后的数据 伪时序分析通常从已经完成注释的对象开始。这里仍然沿用前面教程的文件结构。对于 PBMC 示例数据，它更适合展示代码结构，而不是强行解释为真实发育轨迹。
import os import numpy as np import pandas as pd import scanpy as sc import seaborn as sns import matplotlib.pyplot as plt from IPython.display import display 读取注释后的对象：
adata = sc.read_h5ad( "/home/data/t090639/project/pbmc_10k/data/processed/07_pbmc_10k_annotated.h5ad" ) adata 检查对象中是否有后续需要的信息：
print(adata.obs.columns) print(adata.obsm.keys()) print(adata.obsp.keys()) print(adata.layers.keys()) print(adata.raw is not None ) 理想情况下，对象中应至少包含：
adata.obs["cell_type"] ：细胞类型注释。 adata.obsm["X_pca"] ：PCA 坐标。 adata.obsm["X_umap"] ：UMAP 坐标。 adata.obsp["connectivities"] 和 adata.obsp["distances"] ：邻居图。 adata.raw ：完整 Log(CP10k+1) 表达矩阵。 adata.layers["counts"] ：原始 UMI count。 如果是多样本整合后的对象，还可能包含：
adata.obsm["X_pca_harmony"] ：Harmony 校正后的 PCA 表示。 adata.obs["sample"] 、 adata.obs["batch"] 、 adata.obs["condition"] ：样本和实验设计信息。 选择用于轨迹分析的细胞群 轨迹推断前要先选择一个生物学上合理的细胞群。不要把所有注释后的细胞类型直接放在一起排序。
例如，如果研究对象是造血分化，可以选择 HSC、progenitor 和相关终末细胞；如果研究对象是 T cell 激活或耗竭，可以选择 T cell 相关亚群；如果研究对象是上皮成熟，可以选择同一上皮谱系内的 progenitor、intermediate 和 mature cells。
这里给一个通用模板。实际分析时需要把 lineage_cell_types 替换成当前项目中的相关细胞类型。
cell_type_key = "cell_type" if cell_type_key not in adata.obs.columns: raise ValueError( f"当前对象中没有找到
{cell_type_key} ，请先完成细胞类型注释。" ) lineage_cell_types = [ "CD14+ monocytes" , "FCGR3A+ monocytes" , "Dendritic cells" , ] lineage_cell_types = [ x for x in lineage_cell_types if x in adata.obs[cell_type_key].astype(str).unique() ] if len(lineage_cell_types) == 0 : raise ValueError( "没有找到可用于轨迹分析的细胞类型，请重新设置 lineage_cell_types。" ) adata_traj = adata[ adata.obs[cell_type_key].astype(str).isin(lineage_cell_types) ].copy() adata_traj 如果 adata_traj.n_obs 很少，或者 lineage_cell_types 为空，说明当前对象并不适合直接使用这组细胞类型，需要重新选择谱系。
选择后先查看这些细胞在 UMAP 上的位置：
sc.pl.umap( adata_traj, color=[cell_type_key], legend_loc= "right margin" ) 这一步只是检查选择的细胞是否在表达空间中形成相对连续或相邻的结构。不能因为 UMAP 上看起来连成一条线，就直接证明它们存在真实发育关系。
重新构建轨迹子集的邻居图 把细胞 subset 出来之后，建议针对这个子集重新运行 PCA、邻居图和 UMAP。原因是：全数据对象中的邻居图是为所有细胞类型构建的，子集内部的局部结构可能需要重新估计。
如果当前对象有 adata.raw ，可以先从完整 Log(CP10k+1) 表达矩阵恢复子集，然后重新选择高变基因和 PCA。这里应从已经 subset 后的 adata_traj.raw 恢复，而不是直接对原始 adata.raw 按细胞名切片；部分 anndata 版本中，后者可能导致 .raw.X 和 .obs 行数不一致。
if adata_traj.raw is not None : obs_keep = adata_traj.obs.copy() obsm_keep = { key: adata_traj.obsm[key].copy() for key in adata_traj.obsm.keys() } adata_traj = adata_traj.raw.to_adata() adata_traj.obs = obs_keep for key, value in obsm_keep.items(): adata_traj.obsm[key] = value 重新筛选高变基因：
sc.pp.highly_variable_genes( adata_traj, n_top_genes= 2000 , flavor= "seurat" ) adata_traj.raw = adata_traj.copy() adata_traj = adata_traj[:, adata_traj.var[ "highly_variable" ]].copy() Scaling 和 PCA：
sc.pp.scale(adata_traj, max_value= 10 ) sc.pp.pca( adata_traj, n_comps= 50 , svd_solver= "arpack" , random_state= 0 ) 选择 PC 数：
sc.pl.pca_variance_ratio(adata_traj, n_pcs= 50 , log= True ) n_pcs = 30 构建邻居图并重新运行 UMAP：
sc.pp.neighbors( adata_traj, n_neighbors= 15 , n_pcs=n_pcs, random_state= 0 ) sc.tl.umap( adata_traj, min_dist= 0.3 , spread= 1.0 , random_state= 0 ) 画图检查：
plot_cols = [cell_type_key] for col in [ "total_counts" , "n_genes_by_counts" , "pct_counts_mt" , "sample" , "batch" ]: if col in adata_traj.obs.columns: plot_cols.append(col) sc.pl.umap( adata_traj, color=plot_cols, ncols= 2 , legend_loc= "right margin" )
如果子集 UMAP 主要被 QC 指标或批次驱动，需要先回到质控、整合或细胞选择步骤，而不是继续计算 pseudotime。
使用 PAGA 检查轨迹拓扑 在计算 pseudotime 前，可以先用 PAGA 查看 cluster 或 cell type 之间的连接结构。
如果子集中还没有合适的聚类标签，可以重新运行 Leiden：
sc.tl.leiden( adata_traj, resolution= 0.5 , key_added= "traj_leiden_res0_5" , flavor= "igraph" , n_iterations= 2 , random_state= 0 ) 选择用于 PAGA 的分组列：
paga_group_key = "traj_leiden_res0_5" # 如果细胞类型注释足够细，也可以改成 cell_type_key。 # 但如果 cell_type 过粗，PAGA 会把内部异质性合并掉， # 因此通常先用较细的 Leiden cluster 检查拓扑。 # paga_group_key = cell_type_key adata_traj.obs[paga_group_key] = adata_traj.obs[paga_group_key].astype( "category" ) 运行 PAGA：
sc.tl.paga( adata_traj, groups=paga_group_key ) 查看 PAGA 图：
sc.pl.paga( adata_traj, color=paga_group_key, threshold= 0.03 )
PAGA 图中的边表示不同 group 之间在邻居图上的连接可信度。边越粗，说明两个 group 在细胞图中连接越强。这个连接强度不是 p 值，也不是统计意义上的因果关系，更不是直接证明发育方向。
也可以用 PAGA 初始化 UMAP，让二维图更尽量保留全局拓扑：
sc.tl.umap( adata_traj, init_pos= "paga" , random_state= 0 ) sc.pl.umap( adata_traj, color=[paga_group_key], legend_loc= "right margin" ) 如果 PAGA 显示多个 group 之间几乎没有连接，或者图结构和已有生物学知识完全冲突，就不应该继续强行解释为一条连续轨迹。
选择 root cell 计算 DPT pseudotime 前必须指定 root cell。这个 root cell 应该代表生物过程的起始状态。
常见 root 选择依据包括：
已知 progenitor 或 stem-like cell type。 lineage tracing 或其他独立证据支持的起点。 PAGA 或 diffusion map 中与早期状态一致的端点。 不要只因为某个细胞在 UMAP 最左边或最上边，就把它当作 root。UMAP 坐标轴本身没有生物学方向。
根据已知细胞类型选择 root 如果已有明确的早期细胞类型，可以从该细胞类型中选择一个代表细胞作为 root。
root_cell_type = "CD14+ monocytes" root_mask = adata_traj.obs[cell_type_key].astype(str).eq(root_cell_type) if root_mask.sum() == 0 : raise ValueError( f"当前数据中没有找到 root cell type: {root_cell_type} " ) root_ix = np.flatnonzero(root_mask.to_numpy())[ 0 ] adata_traj.uns[ "iroot" ] = int(root_ix) 这种写法选择的是该细胞类型中的第一个细胞。更严谨的做法是结合 marker score 或 diffusion component，选择更接近早期状态中心的细胞。
根据 marker gene score 选择 root 如果起始状态有一组 marker gene，可以先计算 root score，再选择 score 最高的细胞作为 root。
root_marker_genes = [ "LYZ" , "S100A8" , "S100A9" ] use_raw = adata_traj.raw is not None gene_names = adata_traj.raw.var_names if use_raw else adata_traj.var_names root_marker_genes = [gene for gene in root_marker_genes if gene in gene_names] if len(root_marker_genes) == 0 : raise ValueError( "没有在当前对象中找到 root marker genes。" ) sc.tl.score_genes( adata_traj, gene_list=root_marker_genes, score_name= "root_score" , use_raw=use_raw ) root_scores = adata_traj.obs[ "root_score" ].to_numpy().copy() if "root_cell_type" in globals() and root_cell_type
in adata_traj.obs[cell_type_key].astype(str).unique(): candidate_mask = adata_traj.obs[cell_type_key].astype(str).eq(root_cell_type).to_numpy() root_scores[~candidate_mask] = -np.inf root_ix = int(np.nanargmax(root_scores)) adata_traj.uns[ "iroot" ] = root_ix adata_traj.obs[ "is_root_cell" ] = False adata_traj.obs.iloc[root_ix, adata_traj.obs.columns.get_loc( "is_root_cell" )] = True 检查 root score 分布：
sc.pl.umap( adata_traj, color=[ "root_score" , "is_root_cell" , cell_type_key], ncols= 3 , legend_loc= "right margin" ) 如果 root score 最高的细胞是明显低质量细胞、doublet 或不在预期起始群中，需要重新定义 root marker 或先过滤异常细胞。
使用 DPT 计算 pseudotime DPT 需要先计算 diffusion map。
sc.tl.diffmap( adata_traj, n_comps= 15 ) 查看 diffusion map：
diffmap_colors = [cell_type_key] if "root_score" in adata_traj.obs.columns: diffmap_colors.append( "root_score" ) sc.pl.diffmap( adata_traj, color=diffmap_colors, components=[ "1,2" , "2,3" ], ncols= 2 )
diffusion components 可以帮助观察非线性连续结构。某些数据中，起始细胞会位于某个 diffusion component 的端点。但这仍然需要结合 marker gene 和生物学知识判断。
确认已经设置 root：
adata_traj.uns[ "iroot" ] 计算 DPT：
sc.tl.dpt( adata_traj, n_dcs= 10 ) adata_traj.obs[ "dpt_pseudotime" ] = ( adata_traj.obs[ "dpt_pseudotime" ] .replace([np.inf, -np.inf], np.nan) ) 计算完成后，结果会写入：
adata_traj.obs[ "dpt_pseudotime" ].head() 可视化 pseudotime：
sc.pl.umap( adata_traj, color=[ "dpt_pseudotime" , cell_type_key], color_map= "gnuplot" , ncols= 2 , legend_loc= "right margin" ) 也可以查看不同细胞类型或 cluster 的 pseudotime 分布：
sc.pl.violin( adata_traj, keys= "dpt_pseudotime" , groupby=cell_type_key, rotation= 45 ) 如果发现本应处于早期的细胞整体 pseudotime 很高，或者终末状态比起始状态更接近 0，通常说明 root 选择不合理，或者轨迹拓扑不符合当前方法假设。
检查 pseudotime 是否合理 伪时序结果不能只看一张 UMAP。至少需要做几类检查。
检查是否符合细胞类型层级 按细胞类型统计 pseudotime：
pt_summary = ( adata_traj.obs .groupby(cell_type_key, observed= True )[ "dpt_pseudotime" ] .agg([ "count" , "median" , "mean" , "min" , "max" ]) .sort_values( "median" ) ) pt_summary 如果已知早期群的 median pseudotime 不在前面，或已知终末群分布非常混乱，需要重新检查 root、cell subset、邻居图和批次效应。
检查是否被 QC 指标驱动 qc_cols = [ col for col in [ "total_counts" , "n_genes_by_counts" , "pct_counts_mt" ] if col in adata_traj.obs.columns ] if len(qc_cols) > 0 : sc.pl.umap( adata_traj, color=[ "dpt_pseudotime" ] + qc_cols, ncols= 2 , color_map= "viridis" ) display( adata_traj.obs[[ "dpt_pseudotime" ] + qc_cols] .corr(method= "spearman" ) ) 如果 pseudotime 和 total_counts 、 n_genes_by_counts 或 pct_counts_mt 高度相关，需要警惕它可能主要反映技术梯度，而不是生物过程。
检查是否被批次或样本驱动 sample_cols = [col for col in [ "sample" , "batch" , "condition" ] if col in adata_traj.obs.columns] for col in sample_cols: print(col) sc.pl.violin( adata_traj, keys= "dpt_pseudotime" , groupby=col, rotation= 45 ) 如果某个样本几乎全部处于早期，而另一个样本几乎全部处于晚期，需要结合实验设计判断这是真实生物差异，还是样本处理、批次或细胞组成差异。
对于病例/对照或处理组/对照组，不能直接把每个细胞当作独立重复来比较 pseudotime。更稳妥的做法是按样本统计 pseudotime 分布、阶段比例或 pseudobulk 表达，再在样本层面做统计。
检查 root 选择的敏感性 伪时序结果对 root cell 可能很敏感。可以选择几个候选 root，分别计算 pseudotime，再比较结果是否稳定。
candidate_root_cells = {} for group in adata_traj.obs[cell_type_key].astype(str).unique(): group_mask = adata_traj.obs[cell_type_key].astype(str).eq(group) candidate_root_cells[group] = adata_traj.obs_names[group_mask][ 0 ] candidate_root_cells 对每个候选 root 计算 DPT：
pt_keys = [] for root_name, root_cell in candidate_root_cells.items(): adata_tmp = adata_traj.copy() root_ix = np.flatnonzero(adata_tmp.obs_names == root_cell)[ 0 ] adata_tmp.uns[ "iroot" ] = int(root_ix) sc.tl.dpt( adata_tmp, n_dcs= 10 ) safe_root_name = root_name.replace(
" " , "_" ).replace( "+" , "pos" ) pt_key = f"dpt_root_ {safe_root_name} " adata_traj.obs[pt_key] = adata_tmp.obs[ "dpt_pseudotime" ].to_numpy() pt_keys.append(pt_key) 比较不同 root 得到的 pseudotime：
adata_traj.obs[pt_keys].corr(method= "spearman" ) 如果不同 root 的结果只是方向相反，说明轨迹可能相对稳定，但方向需要由生物学决定。如果不同 root 得到完全不同的排序，说明当前数据的轨迹结构不稳定，不应过度解释。
查看基因随 pseudotime 的变化 伪时序分析最常见的下游任务，是观察基因表达是否沿 pseudotime 连续变化。
先准备一组要查看的基因。实际项目中应该换成和目标过程相关的 marker、转录因子、通路基因或差异基因。
trend_genes = [ "LYZ" , "S100A8" , "S100A9" , "FCGR3A" , "MS4A7" ] use_raw = adata_traj.raw is not None gene_names = adata_traj.raw.var_names if use_raw else adata_traj.var_names trend_genes = [gene for gene in trend_genes if gene in gene_names] trend_genes 先在 UMAP 上查看表达：
if len(trend_genes) > 0 : sc.pl.umap(
color=[ "dpt_pseudotime" ] + trend_genes, use_raw=use_raw, ncols= 3 , vmax= "p99" , color_map= "viridis" ) 再按 pseudotime 分箱，画平均表达趋势：
if len(trend_genes) > 0 : pt_df = sc.get.obs_df( adata_traj, keys=[ "dpt_pseudotime" ] + trend_genes, use_raw=use_raw ) pt_df = pt_df.replace([np.inf, -np.inf], np.nan).dropna(subset=[ "dpt_pseudotime" ]) pt_df[ "pt_bin" ] = pd.qcut( pt_df[ "dpt_pseudotime" ], q= 20 , duplicates= "drop" ) trend_df = ( pt_df .groupby( "pt_bin" , observed= True )[[ "dpt_pseudotime" ] + trend_genes] .mean() .reset_index() ) trend_long = trend_df.melt( id_vars=[ "pt_bin" , "dpt_pseudotime" ], value_vars=trend_genes, var_name= "gene" , value_name= "expression" ) plt.figure(figsize=( 6 , 4 )) sns.lineplot( data=trend_long, x= "dpt_pseudotime" , y= "expression" , hue= "gene" , marker= "o" ) plt.xlabel( "DPT pseudotime" ) plt.ylabel( "Mean log-normalized expression" ) plt.tight_layout() 这种图适合探索表达趋势，但不等同于严格统计检验。原因是 pseudotime 本身是从同一份表达矩阵估计出来的，如果再直接对每个细胞做基因和 pseudotime 的显著性检验，p 值往往会过于乐观。
更严谨的基因趋势分析可以考虑：
在 R 中用
tradeSeq 对 Slingshot 或其他 lineage 的 pseudotime 拟合 GAM。 在 Monocle3 中用 graph_test() 或 fit_models() 分析随轨迹变化的基因。 对多样本数据，在样本层面汇总阶段、细胞比例或 pseudobulk 表达，再做统计比较。 分支结构怎么处理？ 如果 PAGA、Monocle3 或 Slingshot 显示存在分叉，就不能简单把所有细胞压成一条线性 pseudotime。
分叉场景下需要先回答：
某个基因是在所有 lineage 中共同变化，还是只在某条分支中变化。 在 Scanpy 中，PAGA 可以帮助观察分支拓扑，但它本身不直接给每条 lineage 拟合平滑曲线。如果主要问题是分支差异和 lineage-specific gene dynamics，通常可以考虑 Slingshot + tradeSeq 或 Monocle3。
多样本和条件比较中的注意事项 如果要比较不同条件下的 pseudotime 分布，不能简单把所有细胞当成独立样本。单细胞层面的细胞数很多，但真正的生物学重复通常是 donor、mouse、patient 或 sample。
更稳妥的分析方式包括：
每个样本分别统计 pseudotime 分布的中位数、均值或分位数。 把 pseudotime 分成 early、middle、late 阶段后，按样本统计各阶段细胞比例。 在每个样本内按 lineage 和阶段做 pseudobulk，再比较基因表达。 示例：按样本统计 pseudotime 阶段比例。
if "sample" in adata_traj.obs.columns: adata_traj.obs[ "pseudotime_stage" ] = pd.qcut( adata_traj.obs[ "dpt_pseudotime" ], q= 3 , labels=[ "early" , "middle" , "late" ] ) stage_by_sample = pd.crosstab( adata_traj.obs[ "sample" ], adata_traj.obs[ "pseudotime_stage" ], normalize= "index" ) stage_by_sample 如果每个 condition 只有一个样本，不能可靠地区分条件差异和样本差异。这和差异表达、细胞比例分析一样，是实验设计限制，不是轨迹算法能解决的问题。
保存 保存完成伪时序分析后的子集对象：
out = "/home/data/t090639/project/pbmc_10k/data/processed/09_pbmc_10k_pseudotime_dpt.h5ad" adata_traj.uns[ "trajectory_method" ] = "DPT" adata_traj.uns[ "trajectory_root_index" ] = int(adata_traj.uns[ "iroot" ]) adata_traj.uns[ "trajectory_group_key" ] = paga_group_key adata_traj.write_h5ad(out) os.path.exists(out) 保存后的对象中应至少包含：
adata_traj.obs["dpt_pseudotime"] ：DPT 伪时序。 adata_traj.obs["root_score"] ：如果使用 marker score 选择 root。 adata_traj.obsm["X_diffmap"] ：diffusion map 坐标。 adata_traj.uns["iroot"] ：DPT 使用的 root cell 位置。 adata_traj.uns["paga"] ：PAGA 拓扑结果。 adata_traj.raw ：用于基因表达趋势展示的 Log(CP10k+1) 表达矩阵。 常见问题 能不能直接对所有细胞做伪时序？ 一般不建议。伪时序应该围绕一个明确的生物过程，而不是把所有细胞类型强行排序。成熟 T cells、B cells、monocytes 和 NK cells 之间通常不是一条连续分化轨迹。
UMAP 上看起来像一条线，是否就可以做伪时序？ 不一定。UMAP 是二维可视化，会扭曲全局距离和方向。伪时序应基于邻居图、diffusion map、PAGA、principal graph 或 principal curves 等模型，并结合 marker gene、root 选择和生物学知识验证。
pseudotime 数值能不能解释成真实时间？ 不能。pseudotime 是相对顺序，不是小时、天或真实发育时间。即使有真实时间点，也只能说 pseudotime 是否和时间点大致一致，不能直接把 pseudotime 换算成实际时间。
根细胞选错会怎么样？ 可能导致方向反转，也可能导致完全错误的排序。根细胞选择应尽量来自已知早期 marker、时间点、实验设计或其他独立证据，并通过多个 root 的敏感性分析检查稳定性。
整合后的数据能不能做伪时序？ 可以用整合后的 embedding 构建邻居图和轨迹拓扑，但需要检查整合是否改变了目标生物过程。如果条件、时间点或发育阶段本身和批次高度相关，过度校正可能会削弱甚至扭曲真实轨迹。
基因表达趋势不应该从整合 embedding 中解释，而应回到未被整合算法改写的表达矩阵。多样本条件比较仍应回到样本层面，避免把细胞当成独立生物学重复。
能不能用 pseudotime 找差异基因？ 可以探索基因随 pseudotime 的变化，但显著性检验要谨慎。因为 pseudotime 是从表达矩阵估计出来的，直接逐细胞做基因与 pseudotime 的检验容易低估不确定性。正式分析应考虑 GAM、tradeSeq、Monocle3 的图相关检验，或基于样本重复的 pseudobulk 策略。
DPT、Monocle3、Slingshot 结果不一致怎么办？ 这很常见。不同方法的假设、低维表示、拓扑模型和 root 处理方式都不同。应优先判断哪个结果更符合独立生物学证据，而不是简单选择图最好看的方法。
