The RNA velocity are opposite in 'dynamical' and 'stochastic' modes.

Question

[Dan]

Hi

When I run velocity using both 'dynamical' and 'stochastic' modes:

import scvelo as scv
import scanpy as sc

scv.pp.moments(adata)
scv.tl.velocity(adata, mode='stochastic')
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding_stream(adata, basis='umap', 
                                 arrow_size = 1, 
                                 density = 2,
                                 size = 50, 
                                 arrow_style = '->',
                                 color = 'leiden_r1',
                                 alpha = 0.2,
                                 dpi = 300,
                                 legend_loc = 'on data',
                                 integration_direction = 'both',
                                 arrow_color = 'k',
                                 figsize = (6,6))

scv.tl.recover_dynamics(adata)
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding_stream(adata, basis="umap", 
                                 arrow_size = 1, 
                                 density = 2,
                                 size = 50, 
                                 arrow_style = '->',
                                 color = 'leiden_r1',
                                 alpha = 0.2,
                                 dpi = 300,
                                 legend_loc = 'on data',
                                 integration_direction = 'both',
                                 arrow_color = 'k',
                                 figsize = (6,6))

Why the arrow directions are opposite in 'dynamical' and 'stochastic' modes? Which one is correct?

Thanks a lot

[Chris]

Hello Dan. I also run this analysis. May I know the command you run to get the loom file? I got error when try velocyto.

[Dan]

import os
filenames = []
for filename in os.listdir('velocyto'):
    filenames.append('velocyto/'+filename)

import loompy
loompy.combine(filenames, 'velocyto/TY_velocyto.loom', key="Accession")
ldata = scv.read('velocyto/TY_velocyto.loom', cache=False)

import pickle

# Getting back the objects:
with open('adata .pkl', 'rb') as f: 
    adata = pickle.load(f)
adata .obs['sampleindex'] = adata .obs[["sample", "index"]].apply(":".join, axis=1)
adata .obs['sampleindex'] = [si+'x' for si in adata .obs['sampleindex']]
adata .obs.set_index('sampleindex', drop=True, inplace=True)

adata2 = scv.utils.merge(adata, ldata)
adata = adata2

[Chris]

Thank you. Seem you already have the loom file. I am trying to get it with velocyto. This is the command I ran:

velocyto run10x -m Hg38_rmsk.gtf WT/outs /cellranger/reference/refdata-cellranger-GRCh38-3.0.0/genes/genes.gtf

ERROR - This is an older version of cellranger, cannot check if the output are ready, make sure of this yourself. ERROR - Can not locate the barcodes.tsv file!
line 91, in run10x.
bcfile = bcmatches[0].
IndexError: list index out of range.

[Dan]

I think you should use this command: velocyto run10x -m Hg38_rmsk.gtf WT /cellranger/reference/refdata-cellranger-GRCh38-3.0.0/genes/genes.gtf you should use the root directory of cellranger count, not the outs subdirectory.

[Chris]

I tried it and got the same error unfortunately.

[Dan]

you should use the root directory of cellranger count, not the outs subdirectory.

[Chris]

I also tried to unzip the barcodes.tsv.gz file or added the full path to the WT folder. After upload the bam file and rename to possorted_genome_bam.bam, I got:

FileNotFoundError: [Errno 2] No such file or directory: 'samtools': 'samtools' My data is multiome, not single cell.

Answer 1

[james.hawley]

Conflicting directions in RNA velocity analysis is not uncommon. Gorin _et al._, 2022 describes the following when looking at velocity directions un embeddings:

the directions reported in embeddings were unrepresentative of the actual velocity magnitudes in high-dimensional space, as well as severely distorted relative to the linear projection. These discrepancies are a potential cause for concern.

even if we have “true” velocity directions, they only contain a limited amount of highly local information. As expected from the fair performance in Fig 9eiii, the inferred linear embedding (Fig 10b) was globally and locally (Fig 10f) faithful: the model precisely matches the assumptions of the parameter inference workflow. However, the estimates were rapidly distorted upon applying the nonlinear embedding procedure (Fig 10c), rotating many cell-specific directions and suggesting transitions from the green reverting population to the light blue perturbed population, whereas, the true trajectory is from light blue to green. The results of the Boolean procedure were slightly more faithful to the linear projection (Fig 10f) but otherwise qualitatively similar (Fig 10d). This is the method’s best-case performance.

Unfortunately, I can't point to the precise reason why this may be the case for the stochastic and dynamic modes of scvelo. But this may be an effect of the embedding itself from the high-dimensional velocity estimates. It might not be the case that one mode is "correct" and the other is "incorrect". Both of them might be incorrect, too, depending on what your broader question is and how the data was processed to get to this point.