Question about data interpretation

[sneddonucsf]

Hello,

Thank you for the great package!
I am trying to use CellRank on my dataset of developing human fetal endocrine cells. Other pseudotime methods (Slingshot) has constructed the lineages as such:

So 4 terminally differentiated hormone+ cells, and 4 progenitor populations (clusters 2.0.0 -> 2.1.1) with varying degrees of differentiation potency. Because RNA Velocity and CytoScape do not work for my data, I am attempting to use the pseudotime kernel on my data with diffusion pseudotime from Scanpy. First, I've converted my Seurat object to an Adata object and then calculated psedotime with the following code:

adata.uns['iroot'] = np.flatnonzero(adata.obs['annotation']  == 'cEP')[0]
sc.tl.dpt(adata, n_branchings = 0)

This gives me the following pseudotime measurements:

As you can see, the pseudotime values for the beta cells is low compared to the rest of the hormone+ cells. For reference, PAGA looks as the following:

I've then followed through the rest of the pseudotime kernel/estimators:

pk = cr.tl.kernels.PseudotimeKernel(adata)
pk.compute_transition_matrix(threshold_scheme='soft')
scv.pl.velocity_embedding_stream(
    adata, color="annotation", vkey= "T_fwd", basis="umap", legend_loc="right", add_outline = True)

This generally agrees with our Slingshot lineage reconstruction

from cellrank.tl.estimators import GPCCA
g_fwd = GPCCA(pk)
g_fwd.compute_schur(n_components=20)
g_fwd.plot_spectrum(real_only=True)

This gives me a eigenvalue of 3, but calculating the macrostates with this value doesn't quite capture what I think are true macrostates:

So I calculated 5 macrostates, essentially the beginning of my trajectories (cEP) and the terminally differentiated cells types:

g_fwd.set_terminal_states_from_macrostates(["Alpha", "Beta", "Epsilon", "Delta"])
g_fwd.plot_coarse_T()

The coarse plot looks:

with the cEP and Beta clusters being denoted as early states. I would not expect beta cells to be denoted as such. When I calculate the absorption probabilities, I get:

g_fwd.plot_terminal_states(same_plot=True)
g_fwd.compute_absorption_probabilities()
g_fwd.compute_absorption_probabilities(time_to_absorption="all")
g_fwd.plot_absorption_probabilities()

I have the following questions:
1.) Can the pseudotime kernel/estimator results be used for other downstream applications such as directed PAGA, circular projection, plotting lineages, etc.?
2.) In the absorption probabilities plot, is the interpretation that cells that are colored as the terminal state cluster are "fated" to become that cell type?
3.) The absorption probabilities for Beta cells is low as compared to the other terminal states - I'm guessing that this is due to the lower pseudotime values as compared to the other hormone+ cells. Would you recommend trying to get better pseudotime reconstruction and redoing the kernel/estimators? How can I infer lineages?
Thank you very much!

[Marius1311]

Hi @sneddonucsf, thanks for using Cellrank. First question from my side; if you have a slingshot pseudotime and it works well on your data, then why don't you just use that for the pseudotime kernel? You may use any pseudotime method that works for your data!

With respect to your specific questions...

1. yes
2. yes; you can view this more specifically by visualizing one trajectory per plot, check the API (I think same_plot=False)
3. yes, just use your slingshot pseudotime if you're happy with it; that should improve your fate probabilities.

[sneddonucsf]

Hi @Marius1311, thank you for the reply! I've tweaked my pseudotime calculation and now the data looks much better:

Sorry, I have some more questions:
1.) For the macrostate computation, how many "states" would one expect in biologically? I know that you can calculate this and look at the eigenvalues, but in my case the value was 3, which would not give me the expected terminal states. Is it valid to increase this number based on biological interpretation? And would we expect the most initial progenitor (for example, the Ngn3_Low population in the mouse pancreas tutorial) to be a macrostate? Just wondering if I should increase the value to capture all of my terminal states as well as the most initial progenitor.
2.) For directed PAGA with pseudotime data, is it required to have the initial_state? I know that in the directed PAGA tutorial uses cellrank.tl.initial_states(), but I believe this is calculated with velocity information?

Thanks again!

[Marius1311]

Hi!
Re 1: I would say it depends what your primary interest is. If this is a well studied system and you don't expect to find novel initial/terminal states (which I believe is the case), then I would argue it's fine to set the initial and terminal states in a semi-automatic fashion by computing macrostates and manually declaring some of them to be initial/terminal. This enables you to compute fate probabilities which in your case is probably most relevant - it allows you to zoom in on regulation, which genes may be implicated in Alpha vs. Beta differentiation for example?
Re 2: I'm not exactly sure how to do this and I think I've never tried this, please have a look trough closed issues whether anyone tried this before and if not feel free to open a new one specifically for this.

[sneddonucsf]

@Marius1311 That makes sense, in this case we know what to expect for the initial and terminal states, it's more which intermediate states lead to which terminal (so, constructing the lineage). Thank you very much for your help!