How are genes selected for velocity estimation

[mxposed]

Hello, I've got a trivial question:
which genes are used for velocity estimation, leading to future cell state and then transition probabilities?

Are those the highly-variable genes from filter_genes_dispersion? Are they selected based on spliced counts dispersions only or unspliced too? Are they filtered additionally to remove non-informative genes?

[WeilerP]

@mxposed, by default, parameter inference is only run on variables marked as velocity_genes. To identify these, you first filter variables based on counts (filter_genes or via the arguments min_counts, min_counts_u, min_cells, min_cells_u, min_shared_counts, min_shared_cells in filter_and_normalize). After normalization, highly dispersive genes are selected. We store the spliced count matrix in adata.X, i.e. we calculate dispersion on spliced counts only (see here). Finally, this set of genes is further reduced by estimating if they are informative or not. This is based on, for example, their R squared score of the steady state velocity model exceeding a pre-defined threshold. See here for the full definition.
When working with the dynamical model, you can then also filter genes based on their likelihood (e.g. done in latent_time here).

[mxposed]

Thank you! Additionally I appreciate links to the code

[mxposed]

I read the source carefully and explored a bit with the scv.datasets.pancreas() and wanted to write down more details about the selection of velocity genes:

1. As Philipp wrote above, counts are first filtered and normalized, then HVGs are selection (2000 in the pancreas dataset tutorial)
2. Then during recover_dynamics call the steady-state velocity model is fit mode="deterministic"
3. If fit_r2 from this is below default 0.01 then gene is excluded from velocity_genes (n=622)
4. If gamma from this is below default 0.01 then gene is excluded from velocity_genes (n=107). Union of step 3 & 4 is n=705
5. Then dynamical model is fit for remaining genes
6. During velocity call with mode="dynamical" the genes with fit_likelihood <= min_likelihood (default 0.001) are removed from velocity genes (n=17)
7. After that genes with too-high and too-low unspliced_std / spliced_std ratio, called scaling are removed from velocity (n=254)

In step 7 the too-low value is the minimum of 10-th percentile and 0.03, and the too-high value is the maximum of 90-th percentile and 3.

In total this explains how 2000 HVGs produced 1024 velocity genes.

Questions & comments:

1. There's no option to change default values in steps 3 and 4
2. I understand the general logic for step 7, but why 10th and 90th percentiles and values 0.03 and 3?
3. Do I get it right that the boolean velocity_genes mask is affecting the velocity_graph call later on? And so the transition probabilities?
4. Do I get it right that if I want to adjust filtering from step 7, I can just directly assign velocity_genes after the call to sc.tl.velocity?
5. What's the meaning of r2_adjusted here?

[WeilerP]

@mxposed, I'll try to answer your questions and comments individually.

1. There's no option to change default values in steps 3 and 4

Totally agree that it should be possible to manually define thresholds in steps three and four. Regarding step three: You could run the steady-state model yourself and use a custom threshold to define velocity genes. scv.tl.recover_dynamics should then no longer run the steady-state model but use the provided column directly. Currently working on refactoring the code to boost performance and make it more readable and adjustable.

2. I understand the general logic for step 7, but why 10th and 90th percentiles and values 0.03 and 3?

I have honestly no idea how the values were chosen. I assume they were found heuristically by looking at different datasets. @VolkerBergen would know more and could hopefully shed some light on it.

3. Do I get it right that the boolean velocity_genes mask is affecting the velocity_graph call later on? And so the transition probabilities?

I believe so, yes. But would have to check again myself. They already affect the downstream analysis since models are only inferred for those genes (or a subset of them).

4. Do I get it right that if I want to adjust filtering from step 7, I can just directly assign velocity_genes after the call to sc.tl.velocity?

Yes, you could but off the top of my head, I wouldn't do so. I'd manually assign/tweak the velocity gene assignment prior to calling scv.tl.recover_dynamics. I'm afraid that you'd otherwise change the ground truth in the two steps.

5. What's the meaning of r2_adjusted here?

Sorry, cannot offer any insight on this (besides simply stating what the code is doing). @VolkerBergen would again know more.

[mxposed]

Thank you very much for your detailed answers!