How to compare genotypes in Multivelo and assess whether differences are statistically significant?

[yamihn]

Hi, and thank you for developing and maintaining Multivelo — it's a great tool!

I am currently analyzing single-cell data and using Multivelo to infer dynamics across different genotypes. I would like to compare the results between genotypes (e.g., velocity profiles, switch times, gene-level dynamics) and determine whether the differences observed are meaningful or simply due to stochastic variability.

Could you please advise:

What is the recommended way to compare genotypes in Multivelo?

Are there built-in functions or recommended downstream tools to assess whether differences are statistically significant and not just stochastic?

Is there a best practice for visualizing these comparisons?

Thanks in advance for your help!

Best regards

[danielee0707]

Hi yamihn, thank you for your interest in using our tool. For comparing velocity parameters across different conditions (e.g., samples with different genotypes), I recommend taking a look at our new method https://github.com/welch-lab/MultiVeloVAE, which has innate sample integration and batch correction capability and makes cross-sample comparisons simpler. It has functions to perform statistical tests as well. There are example notebooks in that repo on how to use it; feel free to try it out, and we hope to publish it soon. Thanks.

[yamihn]

Thank you for your suggestion and for developing MultiVeloVAE!
I have a question: is it possible to use MultiVeloVAE with a multiomic dataset (e.g., scRNA + scATAC) that has already been integrated using Seurat?
Specifically, can the integrated data objects or matrices produced by Seurat be used as input, or does MultiVeloVAE require starting from raw, non-integrated data? Reading the notebook, it seems that this could be possible, but for MultiVelo I have been having trouble converting the peak annotation and peak linkage files to the required format. Do you have any suggestions or best practices for this conversion, especially when starting from Seurat-integrated data?

Thank you very much for your help!

[danielee0707]

It is possible to use MultiVeloVAE on already integrated data -- basically treating them as one sample and training without batch correction. But I still recommend following the single-sample preprocessing steps to process each sample and concatenate them together later for MultiVeloVAE training (see MultiVeloVAE preprocess examples). The Seurat-based integration doesn't provide guarantees on retaining lineage information across modalities needed during velocity inference.

To merge peak-level accessibility to genes from multiple samples, you can use Signac's merging functions https://stuartlab.org/signac/articles/merging. Although in this case, you may need to recompute feature lineages inside Signac or Cicero and aggregate towards each gene by writing some custom scripts.

[yamihn]

ok, thank you for your help!

[yamihn]

"Although in this case, you may need to recompute feature lineages inside Signac or Cicero and aggregate towards each gene by writing some custom scripts."
Hi again, do you had any advice to generate the bedpe file for linkages?

I run cicero and I generated the peak-peak links with the accessibility score, but the genes are missing, and I went through this kind of error:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[109], line 1
----> 1 adata_atac = mv.aggregate_peaks_10x(adata,
      2                                     'atac_peak_annotation.tsv',
      3                                     'feature_linkage.multivelo.bedpe')
      5 # Ora adata_atac_all è l’oggetto AnnData combinato pronto all’uso!

File ~/miniconda3/lib/python3.12/site-packages/multivelo/auxiliary.py:258, in aggregate_peaks_10x(adata_atac, peak_annot_file, linkage_file, peak_dist, min_corr, gene_body, return_dict, parallel, n_jobs)
    255 gene2 = t3.split('_')
    256 # one of the peaks is in promoter, peaks belong to the
    257 # same gene or are close in distance
--> 258 if (((gene1[1] == "promoter") !=
    259     (gene2[1] == "promoter")) and
    260     ((gene1[0] == gene2[0]) or
    261      (float(tmp[11]) < peak_dist))):
    263     if gene1[1] == "promoter":
    264         gene = gene1[0]

IndexError: list index out of range

This is the head of my self generated bedpe file, where between <> I put the peaks coordinates since I don't have genes info

chr1	10006926	10007771	chr1	9755379	9756277	<chr1:10006926-10007771><chr1:9755379-9756277>	0.0	.	.	.	.	peak-peak
chr1	10006926	10007771	chr1	9768497	9769396	<chr1:10006926-10007771><chr1:9768497-9769396>	0.0	.	.	.	.	peak-peak
chr1	10006926	10007771	chr1	9772615	9773452	<chr1:10006926-10007771><chr1:9772615-9773452>	0.0	.	.	.	.	peak-peak
chr1	10006926	10007771	chr1	9797715	9798617	<chr1:10006926-10007771><chr1:9797715-9798617>	0.0	.	.	.	.	peak-peak
chr1	10006926	10007771	chr1	9847887	9848822	<chr1:10006926-10007771><chr1:9847887-9848822>	0.0	.	.	.	.	peak-peak
chr1	10006926	10007771	chr1	9908349	9909194	<chr1:10006926-10007771><chr1:9908349-9909194>	0.0	.	.	.	.	peak-peak
chr1	10006926	10007771	chr1	9943665	9944562	<chr1:10006926-10007771><chr1:9943665-9944562>	0.0	.	.	.	.	peak-peak
chr1	10006926	10007771	chr1	9947417	9948246	<chr1:10006926-10007771><chr1:9947417-9948246>	0.0	.	.	.	.	peak-peak
chr1	10006926	10007771	chr1	9966992	9967496	<chr1:10006926-10007771><chr1:9966992-9967496>	0.0	.	.	.	.	peak-peak
chr1	10006926	10007771	chr1	9967638	9968418	<chr1:10006926-10007771><chr1:9967638-9968418>	0.0	.	.	.	.	peak-peak
chr1	10006926	10007771	chr1	9981399	9982330	<chr1:10006926-10007771><chr1:9981399-9982330>	0.0	.	.	.	.	peak-peak
chr1	10006926	10007771	chr1	9988243	9989086	<chr1:10006926-10007771><chr1:9988243-9989086>	0.0	.	.	.	.	peak-peak
chr1	10006926	10007771	chr1	9989758	9990620	<chr1:10006926-10007771><chr1:9989758-9990620>	0.0	.	.	.	.	peak-peak

[danielee0707]

The <> pairs should contain gene annotations for each peak. Some example lines of the bedpe file from the 10X mouse brain example look like these

chr1	3119414	3121804	chr1	3903399	3903929	<_intergenic><Gm37381_distal>	0.3832	.	.	8.7036	783055	peak-peak
chr1	3669432	3672890	chr1	4098158	4098427	<Xkr4_promoter><Rp1_distal;Gm37381_distal>	0.3518	.	.	6.1079	427131	peak-peak

You can use Homer http://homer.ucsd.edu/homer/ngs/annotation.html or Signac https://stuartlab.org/signac/reference/annotation to map the regions to nearby genes and fill in the <> values. The peak annotation strategy used by CellRanger is described here https://www.10xgenomics.com/support/software/cell-ranger-arc/latest/analysis/outputs/peak-annotations.

It's also very important that you copy the coaccessibility scores from Cicero (coaccess, between 0 and 1) to the field after <><> instead of setting all values to 0.0, to be used later as a filter when selecting linkages inside the mv.aggregate_peaks_10x function.

[yamihn]

Thank you very much for all your explanations. I will proceed to try MultiveloAE from cellranger matrices and will keep you informed of the results.
Best regards,
Olga

[yamihn]

Hi again! As you suggested, I’m following the MultiVeloAE pipeline starting from preprocessing each sample, concatenating them, and integrating.
But—sorry if I missed it—in this pipeline is there a step where a common set of peaks is selected across all matrices during integration? I have 8 samples, 2 time points, and 2 experimental conditions.

Final question: after integration (using mvv_hspc_integration.ipynb), how can I check which genes follow model 1 or model 2?

Thanks again!

[danielee0707]

You don't need to find a common set of peaks during preprocessing - simply aggregate the peaks towards genes in each sample individually and concatenate ATAC anndata objects together before integration. See this preprocessing notebook.

In the new approach, we generalized the idea of decoupling from Model 1 and Model 2 to continuous cell-wise factors. In short, a positive decoupling value is like Model 2, and a negative decoupling value is like Model 1. You can find decoupling patterns following mvv_hspc_integration.ipynb

adata.layers[f'{key}_decoupling'] = adata.layers[f'{key}_kc'] - adata.layers[f'{key}_rho']
adata.layers[f'{key}_coupling'] = adata.layers[f'{key}_kc'] + adata.layers[f'{key}_rho'] - 1

And you can also test for genes with differential decoupling patterns in a group of cells (such as a cell type) following this mouse brain notebook. Note that if you run integration, you need to modify a line to

df_dd = vv.model.differential_dynamics(adata, adata_atac, model=model, batch_key='batch', idx1=a, idx2=b, mode='change', test_decoupling=True, save_raw=True)

To let it know there is a batch label.

These will be explained clearly in the manuscript. But in the meantime, let me know if you have any questions. Thanks.

[yamihn]

Hi @danielee0707 i have another question for you. A gene could be decoupled in One cluster compared tò another (as I saw in differential decoupling). But I would like tò see the whole umap colored by intensity of coupling / decoupling for a given gene. Is It possible?
Thank you

[yamihn]

Looking into things more deeply, I was wondering how to compare genotypes. Let me explain: I don’t want to perform a differential decoupling analysis — I just want to see how the velocities change based on genotype.
Should I start from the beginning by selecting only the WT cells and running the analysis, then proceed with the KO?
How can I compare the results to ensure that any observed differences are genuin? Furthermore, I think the umap will be differents...
Any help will be welcome!
Thanks

[danielee0707]

The definition of velocity is a little ambiguous here. Do you mean the streamplot that is displayed as arrows on UMAPs, or the actual ds/dt value matrix, or perhaps the fate probabilities computed from cellular transitions? Looking at velocity alone usually isn't that intuitive because it's an instantaneous hidden quantity, whereas changes of cell type transition computed from the velocities are often more biologically relevant. But if comparing the numeric differences of velocities between the samples is exactly what you're looking for, you could use some simple statistics like the Wilcoxon rank-sum or other pseudo-bulk-based approaches -- basically like a differential expression test, after running MultiVeloVAE for each sample. But if you wish to see the cellular transition differences on UMAPs, you could run MultiVeloVAE in integration mode and compute a joint UMAP streamplot.

Because the comparison is between two different samples with no overlap, there is no guarantee that any numeric difference isn't from the technical part of the batch effect, even after batch correction (which can alter true signals too). The risk is always there.

[yamihn]

Thank you. I have another question. Is it possible to run DTW?
can you please drive me through this?
Moreover, about the model, is it possible to reload it and work with it? It was saved as

model.save_model(model_path)
model.save_anndata(data_path, file_name="out.h5ad")

but how can I reload it properly?