Workflow to Explore Trajectories from Root Cluster

[marcouderzo]

Hi,
I am trying to find the possible trajectories from the root node, but I am not entirely sure I'm using the tool right. In the tutorial, it seems that you already know what to expect, since you're listing the root node and the target one.

My workflow has been the following:

1. Custom preprocessing, Leiden clustering and cluster annotation (so, skipped all the stSME/Louvain code)
2. Running pseudotime()

adata.uns["iroot"] = st.spatial.trajectory.set_root(adata,use_label="clusters",cluster=28,use_raw=False)
st.spatial.trajectory.pseudotime(adata,eps=50,use_rep="X_pca",use_label="clusters")

st.pl.trajectory.pseudotime_plot(adata_pseudo, use_label="clusters", edge_alpha = 0.01, node_size=0.1, dpi=300, spot_size=2)

The plot doesn't really show anything useful, not even in umap space.

3. Now, I don't know if I should run pseudotimespace_global() or pseudotimespace_local() for my use case. I am confused on why I should specify a list_clusters=[] argument to construct a specific trajectory. If I don't, it will construct a trajectory with all clusters which doesn't start from my root at all. What is the actual use case of the two methods, and which should I choose if I want to explore possible trajectories between clusters from the root node to either all or a subset of target clusters. I guess what you're showing on the pseudotime_space tutorial using _global makes sense if you already know what you're looking for, right?

4. st.pl.cluster_plot(adata,use_label="clusters", show_trajectories=True,show_subcluster=False) generates a plot with all clusters (of which half of them have the same color and I can't seem to change the cmap to have more than 20 colors), and besides that, only a couple of arrows appear. st.pl.trajectory.tree_plot(data) generates an empty canvas.

Could you please help me understand what's missing or what I am doing wrong? Thank you!

Best,

Marco Uderzo
Bioinformatics Engineer (IE) - Team "Spatio-Temporal Logic of Adult Neurogenesis"
Institute of Functional Genomics (IGFL) - CNRS UMR5242
ENS de Lyon

[newmana]

With many pseudotime methods you don't know what is the start and what is the end of the trajectory. We have gene expression, the graph of connected nodes and distance - none of that tells you which are the start and which are the end - you just have a chain - the start/end is domain specific knowledge - for our example we've picked clusters that represent one step in breast cancer (DCIS) to the next step (IDC).

Monocle 3's tutorial goes through this as well:
https://cole-trapnell-lab.github.io/monocle3/docs/trajectories/#order-cells

I'd say you don't have to know what your clusters are representing but it certainly helps - you can create any number of clusters (2, 10, 100, etc) depending on your parameters (like resolution) - many of which won't necessarily be meaningful.

My understanding of local vs global - local has a single cluster and you're focusing on trajectory within a cluster and global is looking across clusters - that's why you're getting all clusters to all clusters when you don't define any specific ones.

It's a bit worrying you can't see anything useful in UMAP space. If you can start with marker genes, pathologists annotations or cell typing to help you determine which clusters of cells are meaningful.

[newmana]

And if you want all of the trajectories it's inferred you can query it out of adata.uns['available_paths']

[marcouderzo]

Thank you for your explaination. There is another issue which is very strange to me:

adata_gd75.uns["iroot"] = st.spatial.trajectory.set_root(adata_gd75,use_label="clusters",cluster=21,use_raw=False)
st.spatial.trajectory.pseudotime(adata_gd75,eps=50,use_rep="X_pca",use_label="clusters")

lineage = []
for pth in adata_gd75.uns['available_paths'].values():
    if pth[0] == 21: # root
        lineage.append(pth)
print(lineage) 

[[21, 2], [21, 3, 2], [21, 4, 2], [21, 14, 2], [21, 16, 2], [21, 2, 3], [21, 3], [21, 2, 4], [21, 4], [21, 14, 4], [21, 16, 4], [21, 2, 11], [21, 3, 11], [21, 4, 11], [21, 14, 11], [21, 2, 14], [21, 4, 14], [21, 14], [21, 2, 16], [21, 4, 16], [21, 16], [21], [21, 14, 24], [21, 2, 30], [21, 4, 30], [21, 14, 30], [21, 32]]

The path ['21', '2', '11'] exists, so I'll run pseudotimespace_global() on that path

st.spatial.trajectory.pseudotimespace_global(adata_gd75,use_label="clusters", list_clusters=['21', '2', '11'])

Start to construct the trajectory: 11 -> 21 -> 2
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[112], line 1
----> 1 st.spatial.trajectory.pseudotimespace_global(adata_gd75,use_label="clusters", list_clusters=['21', '2', '11'])

File [/opt/conda/envs/spatialdatahd-stlearn/lib/python3.10/site-packages/stlearn/spatial/trajectory/pseudotimespace.py:65](https://192.33.153.151/opt/conda/envs/spatialdatahd-stlearn/lib/python3.10/site-packages/stlearn/spatial/trajectory/pseudotimespace.py#line=64), in pseudotimespace_global(adata, use_label, use_rep, n_dims, list_clusters, model, step, k)
     59 else:
     60     raise ValueError(
     61         "Please choose the right model! Available models: 'mixed', 'spatial' "
     62         + "and 'gene_expression' "
     63     )
---> 65 global_level(
     66     adata,
     67     use_label=use_label,
     68     list_clusters=list_clusters,
     69     w=w,
     70     use_rep=use_rep,
     71     n_dims=n_dims,
     72 )
     74 return adata

File [/opt/conda/envs/spatialdatahd-stlearn/lib/python3.10/site-packages/stlearn/spatial/trajectory/global_level.py:120](https://192.33.153.151/opt/conda/envs/spatialdatahd-stlearn/lib/python3.10/site-packages/stlearn/spatial/trajectory/global_level.py#line=119), in global_level(adata, list_clusters, w, use_label, use_rep, n_dims, return_graph, verbose)
    118 H_sub: networkx.Graph = H.edge_subgraph(edge_list)
    119 if not nx.is_connected(H_sub.to_undirected()):
--> 120     raise ValueError(
    121         "The chosen clusters are not available to construct the spatial "
    122         + "trajectory! Please choose other path."
    123     )
    124 H_sub = nx.DiGraph(H_sub)
    125 prepare_root = []

ValueError: The chosen clusters are not available to construct the spatial trajectory! Please choose other path.

Why is ['21', '2', '11'] becoming 11 -> 21 -> 2 and therefore failing to construct the spatial trajectory?

[newmana]

So it's a bit confusing:

• pseudotime - is classic pseudotime based only on GE - it uses adata.uns['available_paths'] - these are classic results.
• pseudotimespace_global - is an optimal path of pseudotime plus distance - it uses adata.uns['PTS_graph'] - these are an optimal graph based result.

['21', '2', '11'] going to 11 -> 21 -> 2 is the DPT (diffusion pseudotime) ordered version of the pseudotime result and then it checks whether they are spatially connected (and if they're not throws an exception).

Going back to your question 3 - "I guess what you're showing on the pseudotime_space tutorial using _global makes sense if you already know what you're looking for, right?" - this is essentially correct - you have to know what you're looking for.

If you wanted to use it for exploration - you'd have to try setting different root nodes, run pseudo time calls, check they are feasible and then rank them somehow.

[marcouderzo]

Okay so, if I understood correctly, the following workflow should work for exploration:

1. Choose root node (e.g. 21), run pseudotime()
2. For each path from adata.uns['available_paths'] being [root_node , node_x, node_y], run pseudotimespace_global(list_clusters=path).
3. If the DPT ordered version (like ['21', '2', '11'] going to 11 -> 21 -> 2) is spatially connected (so if it doesn't throw exceptions), plot it and save the "11 -> 21 -> 2" version to check the biological feasibility based on differentiation plausibility.
4. If the DPT ordered version does not make sense, discard both "['21', '2', '11']" and "11 -> 21 -> 2" paths.

Is there a way to see all the DPT ordered versions of "available_paths" in a similar way as looking at adata.uns['available_paths']? Meaning, can I extract or print out a list of DPT-refined paths, so that I can do something like:

lineage = []
for pth in adata.uns['DPT_refined_paths'].values():
    if pth[0] == 21: # root
        lineage.append(pth)
print(lineage) 

Thank you very much for the help!

[newmana]

This is a good idea but it doesn’t do that - I’ll keep this issue open as a reminder that this is a change to make.

I like the idea of having exploratory and confirmatory tools.

[marcouderzo]

Thank you! Having a way to check all paths refined by pseudotimespace_global() would be indeed awesome to have implemented!