MultiVeloVAE

https://github.com/welch-lab/multivelovae
Admin
MultiVeloVAE is a deep learning package for velocity inference from multi-lineage, multi-omic, and...
Tokens:32,029
Snippets:246
Trust Score:8.2
License:BSD-3-Clause
Update:2 months ago
Show doc for...
Tokens:
Raw
### Install MultiVeloVAE via pip

Source: https://github.com/welch-lab/multivelovae/blob/main/README.md

Use this command to install the package from PyPI.

```bash
pip install multivelovae
```

--------------------------------

### Importing Libraries

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_hspc_starsolo.ipynb

Initial setup of required data analysis and visualization libraries.

```python
import csv
import numpy as np
import pandas as pd
import scanpy as sc
import scvelo as scv
import matplotlib.pyplot as plt
import multivelovae as vv
```

--------------------------------

### Install MultiVeloVAE via pip

Source: https://github.com/welch-lab/multivelovae/blob/main/docs/usage.md

Use this command to install the package in your Python environment.

```console
$ pip install multivelovae
```

--------------------------------

### Initialize Environment and Imports

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_mouse_brain_scRNA.ipynb

Import necessary libraries for data analysis, visualization, and deep learning.

```python
import numpy as np
import scanpy as sc
import scvelo as scv
import torch
import matplotlib.pyplot as plt
import umap
import os
import multivelovae as vv
from datetime import datetime
import seaborn as sns
import pandas as pd
```

--------------------------------

### Import required libraries

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_mix_integration.ipynb

Initializes the environment with necessary data analysis and machine learning libraries.

```python
import numpy as np
import anndata as ad
import scanpy as sc
import scvelo as scv
import torch
from sklearn.decomposition import PCA
import umap
import os
import multivelovae as vv
from datetime import datetime
```

--------------------------------

### Import required libraries

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_macrophage_integration.ipynb

Initializes the environment by importing necessary data science and multivelovae modules.

```python
import numpy as np
import scanpy as sc
import scvelo as scv
import torch
import matplotlib.pyplot as plt
import umap
import os
import multivelovae as vv
from datetime import datetime
import seaborn as sns
import pandas as pd
import plotnine as p9
```

--------------------------------

### Importing Libraries

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_hspc_macrophage_merge.ipynb

Initial imports for data manipulation, single-cell analysis, and multivelovae integration.

```python
import numpy as np
import pandas as pd
import scanpy as sc
import scvelo as scv
import multivelovae as vv
```

--------------------------------

### Model initialization output

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_EB_multiome.ipynb

Console output indicating the configuration parameters used during model initialization.

```text
Latent dimension set to 16.
Learning rate set to 9.1e-4 based on data sparsity.
Early stop threshold set to 1.6.
Using Gaussian Prior.
Initializing using the steady-state and dynamical models.
```

--------------------------------

### Initialize and train MultiVeloVAE

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_mouse_brain_scRNA.ipynb

Configures the VAEChrom model with specific parameters, trains it, and saves the resulting model and AnnData object.

```python
torch.manual_seed(2022)
np.random.seed(2022)
model = vv.VAEChrom(adata,
                    device='cuda:0',
                    plot_init=False,
                    gene_plot=gene_plot,
                    cluster_key="clusters",
                    figure_path=figure_path,
                    embed="umap")

model.train(plot=False,
            gene_plot=gene_plot,
            figure_path=figure_path,
            embed="umap")

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

--------------------------------

### Initialize and train VAEChrom model

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_macrophage_integration.ipynb

Configures the VAEChrom model with specific batch and regression parameters, trains the model, and saves the output.

```python
torch.manual_seed(2022)
np.random.seed(2022)
model = vv.VAEChrom(adata,
                    adata_atac,
                    batch_key='batch',
                    ref_batch=0,
                    batch_hvg_key='highly_variable',
                    var_to_regress=['total_unspliced3'],
                    device='cuda:0',
                    plot_init=False,
                    gene_plot=gene_plot,
                    cluster_key="leiden",
                    figure_path=figure_path,
                    embed="umap")

model.train(plot=False,
            gene_plot=gene_plot,
            figure_path=figure_path,
            embed="umap")

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

--------------------------------

### Prepare knockout datasets

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_mix_integration.ipynb

Create copies of the layers and set the target gene expression to zero for knockout simulation.

```python
c1 = adata_atac.layers['Mc'].copy()
u1 = adata.layers['Mu'].copy()
s1 = adata.layers['Ms'].copy()
c2 = adata_atac.layers['Mc'].copy()
u2 = adata.layers['Mu'].copy()
s2 = adata.layers['Ms'].copy()
```

```python
c1[:, np.where(adata.var_names == 'SPI1')[0][0]] = 0
u1[:, np.where(adata.var_names == 'SPI1')[0][0]] = 0
s1[:, np.where(adata.var_names == 'SPI1')[0][0]] = 0
c2[:, np.where(adata.var_names == 'GATA1')[0][0]] = 0
u2[:, np.where(adata.var_names == 'GATA1')[0][0]] = 0
s2[:, np.where(adata.var_names == 'GATA1')[0][0]] = 0
```

```python
spi1_ko_set = model.prepare_dataset(c1, u1, s1, only_full=True)
gata1_ko_set = model.prepare_dataset(c2, u2, s2, only_full=True)
```

--------------------------------

### Placeholder for Seurat integration

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_hspc_starsolo.ipynb

Indicates the step for running external Seurat WNN analysis.

```python
### Run seurat_wnn.R ###
```

--------------------------------

### Configuring Analysis Settings

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_hspc_macrophage_merge.ipynb

Sets global verbosity, display options, and figure parameters for scvelo and pandas.

```python
scv.settings.verbosity = 3
scv.settings.presenter_view = True
scv.set_figure_params('scvelo')
pd.set_option('display.max_columns', 100)
pd.set_option('display.max_rows', 200)
np.set_printoptions(suppress=True)
```

--------------------------------

### Reperform Preprocessing Steps

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_hspc_macrophage_merge.ipynb

Applies log transformation, regression of cell cycle differences, scaling, and moment computation for RNA velocity analysis.

```python
# Reperform several preprocessing steps
sc.pp.log1p(adata_)
sc.pp.regress_out(adata_, ['CC_difference'])
sc.pp.scale(adata_)
scv.pp.moments(adata_, n_pcs=30, n_neighbors=30)
sc.pp.log1p(adata2_)
sc.pp.regress_out(adata2_, ['CC_difference'])
sc.pp.scale(adata2_)
scv.pp.moments(adata2_, n_pcs=30, n_neighbors=30)
```

--------------------------------

### Training logs

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_EB_multiome.ipynb

Detailed logs showing ELBO values and timing for training stages.

```text
Output:
2613 out of 3138 = 83.3% genes have good ellipse fits.
KS-test result: [1. 1. 1. 1. 1. 1. 1.]
Assigning cluster 6 to repressive.
Initial induction: 2532, repression: 606 out of 3138.
-------------------------- Train a MultiVeloVAE -------------------------
*********      Creating Training and Validation Datasets      *********
Total Number of Iterations Per Epoch: 12, test iteration: 22
*********                      Finished.                      *********
*********                      Stage  1                       *********
Epoch 1: Train ELBO = -21489.391, Test ELBO = -250868.052		Total Time =   0 h :  0 m :  0 s
Epoch 50: Train ELBO = 14445.737, Test ELBO = 14370.251		Total Time =   0 h :  0 m :  7 s
Epoch 100: Train ELBO = 16401.478, Test ELBO = 16270.854		Total Time =   0 h :  0 m : 15 s
Epoch 150: Train ELBO = 16905.745, Test ELBO = 16748.806		Total Time =   0 h :  0 m : 22 s
Epoch 200: Train ELBO = 17068.374, Test ELBO = 16926.925		Total Time =   0 h :  0 m : 29 s
Epoch 250: Train ELBO = 17170.108, Test ELBO = 17017.487		Total Time =   0 h :  0 m : 36 s
Epoch 300: Train ELBO = 17220.520, Test ELBO = 17036.000		Total Time =   0 h :  0 m : 44 s
Epoch 350: Train ELBO = 17260.815, Test ELBO = 17087.239		Total Time =   0 h :  0 m : 51 s
Epoch 400: Train ELBO = 17289.872, Test ELBO = 17104.871		Total Time =   0 h :  0 m : 58 s
Epoch 450: Train ELBO = 17318.591, Test ELBO = 17122.996		Total Time =   0 h :  1 m :  6 s
Epoch 500: Train ELBO = 17335.009, Test ELBO = 17131.167		Total Time =   0 h :  1 m : 13 s
Epoch 550: Train ELBO = 17353.692, Test ELBO = 17138.741		Total Time =   0 h :  1 m : 20 s
Epoch 600: Train ELBO = 17363.631, Test ELBO = 17144.883		Total Time =   0 h :  1 m : 28 s
Epoch 650: Train ELBO = 17378.854, Test ELBO = 17153.178		Total Time =   0 h :  1 m : 35 s
*********     Stage 1: Early Stop Triggered at epoch 655.     *********
*********     Retrieving best model from iteration 7723.      *********
*********                      Stage  2                       *********
Cell-wise KNN estimation.
Using 148 latent neighbors to select ancestors.
```

--------------------------------

### Initialize MultiVeloVAE Model

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_integration.ipynb

Sets the model identifier key for subsequent operations.

```python
key = "vae"
```

--------------------------------

### Define output paths

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_macrophage_integration.ipynb

Sets base directory paths for saving models, figures, and data.

```python
dataset = "8489-MV-1-9060-MV-3"
model_path_base = f"checkpoints/{dataset}"
figure_path_base = f"figures/{dataset}"
data_path_base = f"data/{dataset}"
```

--------------------------------

### Import MultiVeloVAE

Source: https://github.com/welch-lab/multivelovae/blob/main/docs/api.md

Standard import statement for accessing the library functionality.

```default
import multivelovae as vv
```

--------------------------------

### Visualizing QC Metrics

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_hspc_starsolo.ipynb

Creating violin plots to visualize QC metrics across cells.

```python
sc.pl.violin(adata_, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt', 'pct_counts_rb', 'pct_counts_hla', 'pct_counts_hb', 'pct_counts_in_top_20_genes'], jitter=0.4, multi_panel=True)
```

--------------------------------

### Set Output Paths

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_mouse_brain_scRNA.ipynb

Construct final paths for figures, models, and data using the current date.

```python
figure_path = figure_path_base+'/'+date
model_path = model_path_base+'/'+date
data_path = data_path_base
```

--------------------------------

### Data Preparation and Subsetting

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_mix_integration.ipynb

Assigns cluster labels and creates subsets for scRNA cells.

```python
adata.obs['clusters'] = adata.obs['leiden']
```

```python
# Subset to scRNA cells only for later use
adata_sub = adata[adata.obs['batch']=='scRNA',:].copy()
adata_atac_sub = adata_atac[adata.obs['batch']=='scRNA',:].copy()
```

--------------------------------

### Perform Leiden clustering and visualize

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_EB_starsolo.ipynb

Runs Leiden clustering at a specified resolution and visualizes the resulting clusters on a UMAP.

```python
sc.tl.leiden(adata, resolution=0.7)
sc.pl.umap(adata, color='leiden')
```

--------------------------------

### Define Path Variables

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_mouse_brain_scRNA.ipynb

Set base paths for saving model checkpoints, figures, and data.

```python
model_path_base = f"checkpoints/{dataset}"
figure_path_base = f"figures/{dataset}"
data_path_base = f"data/{dataset}"
```

--------------------------------

### Training progress bar

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_EB_multiome.ipynb

Progress bar display during the initial training phase.

```text
  0%|          | 0/3138 [00:00<?, ?it/s]
```

--------------------------------

### Visualize QC metrics

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_hspc_starsolo.ipynb

Generate violin plots and scatter plots to assess quality control metrics across cells.

```python
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt', 'pct_counts_rb', 'pct_counts_hla', 'pct_counts_hb', 'pct_counts_in_top_20_genes'], jitter=0.4, multi_panel=True)
```

```python
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts', color='pct_counts_mt')
sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt', color='pct_counts_in_top_20_genes')
sc.pl.scatter(adata, x='total_counts', y='pct_counts_rb', color='pct_counts_mt')
sc.pl.scatter(adata, x='pct_counts_mt', y='pct_counts_rb', color='n_genes_by_counts')
```

--------------------------------

### Define Date Variable

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_mouse_brain_scRNA.ipynb

Capture the current date for file naming conventions.

```python
now = datetime.now()
date = now.strftime("%m_%d_%Y")
```

--------------------------------

### Initialize AnnData for perturbation

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_mix_integration.ipynb

Create new AnnData objects to store perturbation results.

```python
adata_spi1 = ad.AnnData(np.zeros((adata.n_obs, adata.n_vars*3)), obs=adata.obs)
adata_spi1.uns = adata.uns
adata_spi1.obsm = adata.obsm
adata_spi1.obsp = adata.obsp
```

```python
adata_gata1 = ad.AnnData(np.zeros((adata.n_obs, adata.n_vars*3)), obs=adata.obs)
adata_gata1.uns = adata.uns
adata_gata1.obsm = adata.obsm
adata_gata1.obsp = adata.obsp
```

--------------------------------

### Calculate QC metrics

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_hspc_starsolo.ipynb

Compute library size, fraction of spliced/unspliced counts, and gene set percentages for quality control.

```python
# Define some extra metrics
counts_u = adata.layers['unspliced'].sum(1).A1
counts_s = adata.layers['spliced'].sum(1).A1
adata.obs['total_unspliced'] = counts_u
adata.obs['total_spliced'] = counts_s
adata.obs['log1p_total_unspliced'] = np.log1p(counts_u)
adata.obs['log1p_total_spliced'] = np.log1p(counts_s)
adata.obs['fraction_u'] = counts_u / (counts_s + counts_u)
```

```python
adata.var['mt'] = adata.var_names.str.startswith('MT-')
adata.var['rb'] = adata.var_names.str.startswith(('RPS', 'RPL'))
adata.var['hla'] = adata.var_names.str.startswith('HLA-')
adata.var['hb'] = np.in1d(adata.var_names, hb_genes)
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt', 'rb', 'hla', 'hb'], percent_top=[20], log1p=True, inplace=True)
```

--------------------------------

### Refinement logs

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_EB_multiome.ipynb

Detailed logs showing refinement rounds and final ELBO results.

```text
Output:
Percentage of Invalid Sets: 0.004
Average Set Size: 17
Finished. Actual Time:   0 h :  0 m :  3 s
*********             Velocity Refinement Round 1             *********
Epoch 678: Train ELBO = 17250.644, Test ELBO = 17073.383		Total Time =   0 h :  1 m : 42 s
*********     Round 1: Early Stop Triggered at epoch 678.     *********
*********     Retrieving best model from iteration 8120.      *********
Cell-wise KNN estimation.
Finished. Actual Time:   0 h :  0 m :  0 s
*********             Velocity Refinement Round 2             *********
Epoch 694: Train ELBO = 17221.109, Test ELBO = 17035.956		Total Time =   0 h :  1 m : 46 s
*********     Round 2: Early Stop Triggered at epoch 694.     *********
*********     Retrieving best model from iteration 8307.      *********
Change in x0: 0.110
Cell-wise KNN estimation.
Finished. Actual Time:   0 h :  0 m :  0 s
*********             Velocity Refinement Round 3             *********
Epoch 706: Train ELBO = 17217.978, Test ELBO = 17019.916		Total Time =   0 h :  1 m : 49 s
*********     Round 3: Early Stop Triggered at epoch 706.     *********
*********     Retrieving best model from iteration 8450.      *********
Change in x0: 0.069
*********      Stage 2: Early Stop Triggered at round 3.      *********
Final: Train ELBO = 17217.978,	Test ELBO = 17019.916
*********      Finished. Total Time =   0 h :  1 m : 49 s     *********
Computing velocity.
Selected 648 velocity genes.
Writing anndata output to file.
```

--------------------------------

### Compute Differential Dynamics with Raw Save

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_macrophage_integration.ipynb

Calculates differential dynamics and saves the raw output.

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

--------------------------------

### Plot phase portraits with chromatin accessibility by batch

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_integration.ipynb

Visualizes phase portraits comparing RNA and chromatin accessibility data colored by batch.

```python
scv.pl.scatter(adata, basis=top_genes[:20], x='Mu', y='Mc', color='batch', ncols=5, frameon=False, fontsize=20, size=30)
scv.pl.scatter(adata, basis=top_genes[:20], x=f'{key}_uhat', y=f'{key}_chat', color='batch', ncols=5, frameon=False, fontsize=20, size=30)
```

--------------------------------

### Configure plotting environment

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_macrophage_integration.ipynb

Sets the seaborn theme for subsequent visualizations.

```python
sns.set_theme(style="white", palette='tab10')
```

--------------------------------

### Configure plotnine settings

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_mouse_brain_multiome.ipynb

Set global DPI and figure size for plotnine visualizations.

```python
import plotnine as p9
p9.options.dpi = 100
p9.options.figure_size = (5, 4)
```

--------------------------------

### Visualize inferred rate parameters

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_integration.ipynb

Generates a 2x2 grid of histograms to compare rate parameters across different conditions. Requires an AnnData object with pre-computed keys.

```python
fig, axs = plt.subplots(2, 2, figsize=(15, 12))
axs[0, 0].hist(adata.var[f'{key}_alpha_c_0'], bins=100, alpha=0.6);
axs[0, 0].hist(adata.var[f'{key}_alpha_c_1'], bins=100, alpha=0.6);
axs[0, 0].set_title('alpha_c');
a1 = adata.var[f'{key}_alpha_0'].values
a1 = a1[a1 < np.percentile(a1, 99.5)]
a2 = adata.var[f'{key}_alpha_1'].values
a2 = a2[a2 < np.percentile(a2, 99.5)]
axs[0, 1].hist(a1, bins=100, alpha=0.6);
axs[0, 1].hist(a2, bins=100, alpha=0.6);
axs[0, 1].set_title('alpha');
axs[1, 0].hist(adata.var[f'{key}_beta_0'], bins=100, alpha=0.6);
axs[1, 0].hist(adata.var[f'{key}_beta_1'], bins=100, alpha=0.6);
axs[1, 0].set_title('beta');
axs[1, 1].hist(adata.var[f'{key}_gamma_0'], bins=100, alpha=0.6);
axs[1, 1].hist(adata.var[f'{key}_gamma_1'], bins=100, alpha=0.6);
axs[1, 1].set_title('gamma');
```

--------------------------------

### Plot inferred rate parameter distributions

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_macrophage_integration.ipynb

Visualizes histograms for alpha, beta, and gamma rate parameters across datasets.

```python
fig, axs = plt.subplots(1, 4, figsize=(30, 5))
axs[0].hist(adata.var[f'{key}_alpha_c_0'], bins=100, alpha=0.6);
axs[0].hist(adata.var[f'{key}_alpha_c_1'], bins=100, alpha=0.6);
axs[0].set_title('alpha_c');
a1 = adata.var[f'{key}_alpha_0'].values
a1 = a1[a1 < np.percentile(a1, 99.5)]
a2 = adata.var[f'{key}_alpha_1'].values
a2 = a2[a2 < np.percentile(a2, 99.5)]
axs[1].hist(a1, bins=100, alpha=0.6);
axs[1].hist(a2, bins=100, alpha=0.6);
axs[1].set_title('alpha');
axs[2].hist(adata.var[f'{key}_beta_0'], bins=100, alpha=0.6);
axs[2].hist(adata.var[f'{key}_beta_1'], bins=100, alpha=0.6);
axs[2].set_title('beta');
axs[3].hist(adata.var[f'{key}_gamma_0'], bins=100, alpha=0.6);
axs[3].hist(adata.var[f'{key}_gamma_1'], bins=100, alpha=0.6);
axs[3].set_title('gamma');
```

--------------------------------

### Save and Load AnnData Objects

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_macrophage_starsolo.ipynb

Save preprocessed AnnData objects to disk and reload them for further analysis.

```python
adata.write_h5ad('9060-MV-3_adata_prepro.h5ad')
adata_atac.write_h5ad('9060-MV-3_adata_atac_prepro.h5ad')
```

```python
adata = sc.read_h5ad('9060-MV-3_adata_prepro.h5ad')
adata_atac = sc.read_h5ad('9060-MV-3_adata_atac_prepro.h5ad')
```

--------------------------------

### Save and Load AnnData Objects

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_EB_starsolo.ipynb

Save preprocessed AnnData objects to disk and reload them for further analysis.

```python
adata.write_h5ad('8964-KL-1_adata_prepro.h5ad')
adata_atac.write_h5ad('8964-KL-1_adata_atac_prepro.h5ad')
```

```python
adata = sc.read_h5ad('8964-KL-1_adata_prepro.h5ad')
adata_atac = sc.read_h5ad('8964-KL-1_adata_atac_prepro.h5ad')
```

--------------------------------

### Test knockout models

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_mix_integration.ipynb

Run the model test function for the knockout datasets.

```python
# Get PU.1 KO results
out1 = model.test(spi1_ko_set, covar=model.var_to_regress.detach().cpu().numpy())
chat1, uhat1, shat1, _, _, _, z1, std_z1, t1, std_t1, t_1, kc1, rho1 = out1
```

```python
# Get GATA1 KO results
out2 = model.test(gata1_ko_set, covar=model.var_to_regress.detach().cpu().numpy())
chat2, uhat2, shat2, _, _, _, z2, std_z2, t2, std_t2, t_2, kc2, rho2 = out2
```

--------------------------------

### Generate dynamic gene expression plots

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_macrophage_integration.ipynb

Creates dynamic plots for marker genes using MultiVeloVAE, supporting batch correction and coloring by metadata.

```python
# Dynamic plots of several marker genes
vv.dynamic_plot(adata, adata_atac, gene_plot, color_by='leiden', batch_correction=True)
```

```python
vv.dynamic_plot(adata, adata_atac, gene_plot, color_by='batch', batch_correction=True)
```

--------------------------------

### Check CUDA Availability

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_mouse_brain_scRNA.ipynb

Verify if a GPU is available for torch computations.

```python
torch.cuda.is_available()
```

--------------------------------

### Plot phase portraits by batch

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_integration.ipynb

Visualizes gene phase portraits colored by batch labels to assess integration.

```python
## Phase portraits compare before and after integration
scv.pl.scatter(adata, basis=top_genes[:20], color='batch', ncols=5, frameon=False, fontsize=20, size=30)
scv.pl.scatter(adata, basis=top_genes[:20], x=f'{key}_shat', y=f'{key}_uhat', color='batch', ncols=5, frameon=False, fontsize=20, size=30)
```

--------------------------------

### Plot phase portraits for top genes by batch

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_macrophage_integration.ipynb

Generates scatter plots and phase portraits for the top 20 genes, colored by batch.

```python
scv.pl.scatter(adata, basis=top_genes[:20], color='batch', ncols=5, frameon=False, fontsize=20, size=30)
scv.pl.scatter(adata, basis=top_genes[:20], x=f'{key}_shat', y=f'{key}_uhat', color='batch', ncols=5, frameon=False, fontsize=20, size=30)
```

--------------------------------

### Visualize Clusters and Prior Time

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_mouse_brain_scRNA.ipynb

Create a UMAP scatter plot colored by clusters and prior time.

```python
os.makedirs(figure_path_base, exist_ok=True)
scv.pl.scatter(adata, basis='umap', color=['clusters', 'tprior'], legend_loc='on data')
```

--------------------------------

### Detect outliers based on QC metrics

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_hspc_starsolo.ipynb

Identify outliers by applying automated thresholds to multiple QC metrics and storing the result in the observation metadata.

```python
adata_.obs["outlier"] = (
    vv.is_outlier(adata_, "log1p_total_counts", 3, 4)
    | vv.is_outlier(adata_, "log1p_n_genes_by_counts", 3, 4)
    | vv.is_outlier(adata_, "pct_counts_in_top_20_genes", 4, 4)
    | vv.is_outlier(adata_, "pct_counts_mt", 4, 3)
    | vv.is_outlier(adata_, "pct_counts_rb", 3, 4)
)
adata_.obs.outlier.value_counts()
```

--------------------------------

### Load Seurat results

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_hspc_starsolo.ipynb

Loads nearest neighbor indices, distances, and cell names generated by Seurat.

```python
nn_idx = np.loadtxt("seurat/nn_idx.txt", delimiter=',')
nn_dist = np.loadtxt("seurat/nn_dist.txt", delimiter=',')
nn_cell_names = list(pd.read_csv("seurat/nn_cell_names.txt", header=None)[0])
nn_idx.shape
```

--------------------------------

### Load STARsolo counts

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_macrophage_starsolo.ipynb

Read 10x Genomics formatted count matrices and assign spliced and unspliced layers.

```python
adata = sc.read_10x_mtx('starsolo-filtered/', var_names='gene_symbols', make_unique=True, cache=True, gex_only=True)
adata.layers['unspliced'] = sc.read_10x_mtx('starsolo-filtered/', prefix='unspliced_', var_names='gene_symbols', make_unique=True, cache=True, gex_only=True).X
adata.layers['spliced'] = sc.read_10x_mtx('starsolo-filtered/', prefix='spliced_', var_names='gene_symbols', make_unique=True, cache=True, gex_only=True).X
adata.obs_names = [x + '-1' for x in adata.obs_names]
```

--------------------------------

### Calculating Quality Control Metrics

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_hspc_starsolo.ipynb

Computing QC metrics including mitochondrial, ribosomal, HLA, and hemoglobin content.

```python
# Compute QC metrics
adata_.var['mt'] = adata_.var_names.str.startswith('MT-')
adata_.var['rb'] = adata_.var_names.str.startswith(('RPS', 'RPL'))
adata_.var['hla'] = adata_.var_names.str.startswith('HLA-')
adata_.var['hb'] = np.in1d(adata_.var_names, hb_genes)
sc.pp.calculate_qc_metrics(adata_, qc_vars=['mt', 'rb', 'hla', 'hb'], percent_top=[20], log1p=True, inplace=True)
```

--------------------------------

### Load Dataset

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_mouse_brain_scRNA.ipynb

Define the dataset name and load the processed AnnData object.

```python
dataset = "VV_braindev"
adata = sc.read_h5ad('Braindev_post.h5ad')
```

--------------------------------

### Map colors to cell types

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_macrophage_integration.ipynb

Creates a dictionary mapping cell types to their assigned colors.

```python
colors = {x:y for x,y in zip(adata.obs['leiden'].cat.categories, adata.uns['leiden_colors'])}
```

--------------------------------

### Generate dynamic plots

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_integration.ipynb

Creates dynamic plots for marker genes with batch correction applied.

```python
# Dynamic plots of several marker genes colored by batch labels
vv.dynamic_plot(adata, adata_atac, gene_plot, color_by='batch', batch_correction=True)
```

--------------------------------

### Visualize Coupling Factors

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_integration.ipynb

Generates scatter plots for chat, uhat, and shat using the coupling factor.

```python
# Couping factors of each modality on dynamic plots
scv.pl.scatter(adata, basis=['AZU1', 'HBD', 'HDC', 'LYZ', 'PF4', 'SPINK2'], x=f'{key}_time', y=f'{key}_chat', color=f'{key}_coupling', dpi=300, legend_loc='none', frameon=False, fontsize=30, color_map='coolwarm', size=40,
               norm=mcolors.CenteredNorm(halfrange=1), colorbar=False, wspace=0.1)
scv.pl.scatter(adata, basis=['AZU1', 'HBD', 'HDC', 'LYZ', 'PF4', 'SPINK2'], x=f'{key}_time', y=f'{key}_uhat', color=f'{key}_coupling', dpi=300, legend_loc='none', frameon=False, fontsize=30, color_map='coolwarm', size=40,
               norm=mcolors.CenteredNorm(halfrange=1), colorbar=False, wspace=0.1)
scv.pl.scatter(adata, basis=['AZU1', 'HBD', 'HDC', 'LYZ', 'PF4', 'SPINK2'], x=f'{key}_time', y=f'{key}_shat', color=f'{key}_coupling', dpi=300, legend_loc='none', frameon=False, fontsize=30, color_map='coolwarm', size=40,
               norm=mcolors.CenteredNorm(halfrange=1), colorbar=False, wspace=0.1)
```

--------------------------------

### Plot phase portraits with chromatin accessibility by leiden

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_integration.ipynb

Visualizes phase portraits comparing RNA and chromatin accessibility data colored by leiden clusters.

```python
scv.pl.scatter(adata, basis=top_genes[:20], x='Mu', y='Mc', color='leiden', ncols=5, frameon=False, fontsize=20, size=30)
scv.pl.scatter(adata, basis=top_genes[:20], x=f'{key}_uhat', y=f'{key}_chat', color='leiden', ncols=5, frameon=False, fontsize=20, size=30)
```

--------------------------------

### Load processed datasets

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_mouse_brain_multiome.ipynb

Reads the RNA and ATAC AnnData objects from disk.

```python
dataset = "10XMouse"
adata = sc.read_h5ad('adata_postpro.h5ad')
adata_atac = sc.read_h5ad('adata_atac_postpro.h5ad')
```

--------------------------------

### Visualize phase portraits

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_EB_multiome.ipynb

Plot scatter plots for top genes to visualize phase portraits and model predictions.

```python
scv.pl.scatter(adata, basis=top_genes[:20], color='leiden', ncols=5)
```

```python
scv.pl.scatter(adata, basis=top_genes[:20], x=f'{key}_shat', y=f'{key}_uhat', color='leiden', ncols=5)
```

```python
adata.layers['Mc'] = adata_atac.layers['Mc'].copy()
```

```python
scv.pl.scatter(adata, x='Mu', y='Mc', basis=top_genes[:20], color='leiden', ncols=5)
```

```python
scv.pl.scatter(adata, basis=top_genes[:20], x=f'{key}_uhat', y=f'{key}_chat', color='leiden', ncols=5)
```

--------------------------------

### Plot phase portraits by leiden cluster

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_integration.ipynb

Visualizes gene phase portraits colored by leiden clusters.

```python
scv.pl.scatter(adata, basis=top_genes[:20], color='leiden', ncols=5, frameon=False, fontsize=20, size=30)
scv.pl.scatter(adata, basis=top_genes[:20], x=f'{key}_shat', y=f'{key}_uhat', color='leiden', ncols=5, frameon=False, fontsize=20, size=30)
```

--------------------------------

### Generate Dynamics Plots by Batch Labels

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_mix_integration.ipynb

Plots gene dynamics with velocity arrows, colored by batch labels for selected genes.

```python
# Dynamics plots for selected genes with velocity arrows, colored by batch labels
batch_raw = adata.obs['batch'].cat.categories.to_numpy()
batch_colors = {batch_raw[i]: adata.uns['batch_colors'][i] for i in range(len(batch_raw))}
np.random.seed(2022)
for gene in gene_plot:
    vv.plot_vel(adata.obs[f'{key}_time'],
                adata[:,gene].layers[f'{key}_chat'],
                adata[:,gene].layers[f'{key}_uhat'],
                adata[:,gene].layers[f'{key}_shat'],
                adata[:,gene].layers[f'{key}_velocity_c'],
                adata[:,gene].layers[f'{key}_velocity_u'],
                adata[:,gene].layers[f'{key}_velocity'],
                cell_labels=adata.obs['batch'].to_numpy(),
                cell_type_colors=batch_colors,
                title=gene,
                frame_on=False,
                axis_on=False,
                legend=False)
```

--------------------------------

### Visualize data with UMAP

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_integration.ipynb

Creates scatter plots on UMAP embeddings colored by leiden clusters and batch.

```python
os.makedirs(figure_path_base, exist_ok=True)
scv.pl.scatter(adata, basis='umap', color='leiden', legend_loc='right margin')
scv.pl.scatter(adata, basis='umap', color='batch', legend_loc='right margin')
```

--------------------------------

### Compare batch-corrected and original reconstructions

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_mix_integration.ipynb

Visualizes the impact of batch correction on gene reconstructions across different color groupings.

```python
# Phase portraits of genes, original, reconstructed and batch-corrected, and reconstructed and not-batch-corrected
scv.pl.scatter(adata, basis=top_genes[:20], color='leiden', ncols=5)
scv.pl.scatter(adata, basis=top_genes[:20], x=f'{key}_shat', y=f'{key}_uhat', color='leiden', ncols=5)
scv.pl.scatter(adata, basis=top_genes[:20], x=f'{key}_shat_batch', y=f'{key}_uhat_batch', color='leiden', ncols=5)
```

```python
scv.pl.scatter(adata, basis=top_genes[:20], color='batch', ncols=5)
scv.pl.scatter(adata, basis=top_genes[:20], x=f'{key}_shat', y=f'{key}_uhat', color='batch', ncols=5)
scv.pl.scatter(adata, basis=top_genes[:20], x=f'{key}_shat_batch', y=f'{key}_uhat_batch', color='batch', ncols=5)
```

--------------------------------

### Analyze Decoupling and Coupling Factors

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_mouse_brain_multiome.ipynb

Calculates decoupling and coupling layers and visualizes them using scatter plots.

```python
# Decouping and coupling factors on UMAPs
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
scv.pl.scatter(adata, basis='umap', color=np.mean(adata[:, adata.var[f'{key}_velocity_genes']].layers[f'{key}_decoupling'], axis=1), legend_loc='right margin')
scv.pl.scatter(adata, basis='z_umap', color=np.mean(adata[:, adata.var[f'{key}_velocity_genes']].layers[f'{key}_decoupling'], axis=1), legend_loc='right margin')
```

--------------------------------

### Visualize QC metrics on UMAP

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_hspc_macrophage_merge.ipynb

Uses scanpy to plot various QC metrics and scores on a UMAP embedding.

```python
sc.pl.umap(adata_concat, color=['n_counts', 'n_Mc', 'n_Mu', 'n_Ms',
                                'S_score', 'G2M_score', 'total_unspliced', 'total_spliced',
                                'fraction_u', 'total_unspliced2', 'total_spliced2', 'fraction_u2',
                                'total_unspliced3', 'total_spliced3', 'fraction_u3',
                                'total_counts', 'pct_counts_mt', 'pct_counts_rb', 'pct_counts_hla',
                                'pct_counts_hb', 'pct_counts_in_top_20_genes', 'CC_difference'])
```

--------------------------------

### Define comparison groups

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_macrophage_integration.ipynb

Creates boolean masks for differential dynamics analysis.

```python
a = np.isin(adata.obs['leiden'], ['M1 Macrophage', 'M2 Macrophage'])
b = np.isin(adata.obs['leiden'], ['DC'])
```

--------------------------------

### Loading and filtering gene expression data

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_EB_starsolo.ipynb

Loading GEX data, filtering cells and genes, and computing QC metrics.

```python
adata_ = sc.read_10x_mtx('filtered_feature_bc_matrix/', var_names='gene_symbols', cache=True, gex_only=True)
adata_.var_names_make_unique()
```

```python
plt.hist(adata_.X.sum(1), bins=100, range=(0, 80000));
```

```python
sc.pp.filter_cells(adata_, min_genes=500)
sc.pp.filter_genes(adata_, min_cells=10)
```

```python
# Compute QC metrics
adata_.var['mt'] = adata_.var_names.str.startswith('MT-')
adata_.var['rb'] = adata_.var_names.str.startswith(('RPS', 'RPL'))
adata_.var['hla'] = adata_.var_names.str.startswith('HLA-')
adata_.var['hb'] = np.in1d(adata_.var_names, hb_genes)
sc.pp.calculate_qc_metrics(adata_, qc_vars=['mt', 'rb', 'hla', 'hb'], percent_top=[20], log1p=True, inplace=True)
```

--------------------------------

### Generate phase portraits

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_mouse_brain_multiome.ipynb

Create scatter plots showing reconstructed relationships between spliced, unspliced, and chromatin data.

```python
# Phase portraits of recontructed unspliced vs spliced for genes of interest, colored by cell type, latent time, and transcription state
scv.pl.scatter(adata, basis=['Robo2', 'Satb2', 'Gria2', 'Grin2b'], x=f'{key}_shat', y=f'{key}_uhat', color='leiden', legend_loc='none', frameon=False, fontsize=20)
scv.pl.scatter(adata, basis=['Robo2', 'Satb2', 'Gria2', 'Grin2b'], x=f'{key}_shat', y=f'{key}_uhat', color=f'{key}_time', legend_loc='none', frameon=False, fontsize=20, color_map='gnuplot')
scv.pl.scatter(adata, basis=['Robo2', 'Satb2', 'Gria2', 'Grin2b'], x=f'{key}_shat', y=f'{key}_uhat', color=f'{key}_rho', legend_loc='none', frameon=False, fontsize=20, color_map='coolwarm')
```

```python
# Phase portraits of recontructed chromatin vs unspliced for genes of interest, colored by cell type, latent time, and transcription state
scv.pl.scatter(adata, basis=['Robo2', 'Satb2', 'Gria2', 'Grin2b'], x=f'{key}_uhat', y=f'{key}_chat', color='leiden', legend_loc='none', frameon=False, fontsize=20)
scv.pl.scatter(adata, basis=['Robo2', 'Satb2', 'Gria2', 'Grin2b'], x=f'{key}_uhat', y=f'{key}_chat', color=f'{key}_time', legend_loc='none', frameon=False, fontsize=20, color_map='gnuplot')
scv.pl.scatter(adata, basis=['Robo2', 'Satb2', 'Gria2', 'Grin2b'], x=f'{key}_uhat', y=f'{key}_chat', color=f'{key}_kc', legend_loc='none', frameon=False, fontsize=20, color_map='coolwarm')
```

--------------------------------

### Filter genes and visualize distributions

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/preprocess/example_prepro_macrophage_starsolo.ipynb

Filter genes by minimum cell count and plot count distributions.

```python
sc.pp.filter_genes(adata, min_cells=10)
```

```python
plt.hist(adata.X.sum(1), bins=100);
```

```python
plt.hist(adata.layers['unspliced'].sum(1), bins=100, alpha=0.5);
plt.hist(adata.layers['spliced'].sum(1), bins=100, alpha=0.5);
```

--------------------------------

### Visualize cell type distribution

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_macrophage_integration.ipynb

Creates a stacked histogram to visualize cell type distribution across batches.

```python
fig, ax = plt.subplots(1, 1, figsize=(25, 6))
sns.histplot(
    data=adata_copy.obs,
    x="leiden",
    hue="batch",
    multiple="stack",
    discrete=True,
    shrink=.7
);
sns.move_legend(ax, "upper left", bbox_to_anchor=(1, 1))
ax.tick_params(axis='x', rotation=35)
ax.set_xlabel('');
```

--------------------------------

### Plot dynamic VC, VU, and VS

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_macrophage_integration.ipynb

Plots dynamic VC, VU, and VS values for top differential velocity genes.

```python
# Dynamic plots of VC, VU, and VS of top differential velocity genes in macrophages
vv.dynamic_plot(adata_sub, adata_atac_sub, df_dd_.sort_values(by=['bayes_factor_v', 'log2_diff_v'], ascending=[False, False])[df_dd_['log2_diff_v']>0].index[:10], color_by='leiden', by='velocity')
```

--------------------------------

### Plot log-transformed mean distributions

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_macrophage_integration.ipynb

Visualizes the distribution of log-transformed mean values for C, U, and S metrics.

```python
plt.hist(np.log(df_dd_mac['mean_c1']), bins=100, alpha=0.7);
plt.hist(np.log(df_dd_mac['mean_u1']), bins=100, alpha=0.7);
plt.hist(np.log(df_dd_mac['mean_s1']), bins=100, alpha=0.7);
```

--------------------------------

### Visualize Expression Distributions

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_macrophage_integration.ipynb

Plots histograms of log-transformed mean expression values for different components.

```python
plt.hist(np.log(df_dd_dc['mean_c1']), bins=100, alpha=0.7);
plt.hist(np.log(df_dd_dc['mean_u1']), bins=100, alpha=0.7);
plt.hist(np.log(df_dd_dc['mean_s1']), bins=100, alpha=0.7);
```

--------------------------------

### Inspect AnnData Object

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_mouse_brain_scRNA.ipynb

Display the structure and contents of the loaded AnnData object.

```python
adata
```

--------------------------------

### Visualize latent time distribution

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_macrophage_integration.ipynb

Displays a violin plot of latent time across different batches.

```python
# Latent time distribution of two datasets
sns.violinplot(y="batch", x=f'{key}_time', data=adata.obs, orient='h');
```

--------------------------------

### Define comparison groups

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_mouse_brain_multiome.ipynb

Define lists of cell types for differential dynamics comparisons.

```python
# Define comparison pairs
a_list = [['V-SVZ'],
          ['Upper Layer'],
          ['Upper Layer', 'Deeper Layer'],
          ['Deeper Layer']]
b_list = [['nIPC', 'Upper Layer'],
          ['V-SVZ', 'Deeper Layer'],
          ['V-SVZ'],
          ['V-SVZ', 'Upper Layer']]
```

--------------------------------

### Plot phase portraits for top genes

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_hspc_mix_integration.ipynb

Visualizes phase portraits for the top 20 genes using scVelo scatter plots.

```python
scv.pl.scatter(adata, basis=top_genes[:20], x=f'{key}_uhat', y=f'{key}_chat', color='batch', ncols=5)
```

```python
adata.layers['Mc'] = adata_atac.layers['Mc']
scv.pl.scatter(adata, basis=top_genes[:20], x='Mu', y='Mc', color='batch', ncols=5)
```

--------------------------------

### Plot phase portraits

Source: https://github.com/welch-lab/multivelovae/blob/main/paper-notebooks/mvv_EB_multiome.ipynb

Generates phase portraits for the selected genes.

```python
# Phase portraits of original data
scv.pl.scatter(adata, basis=gene_plot, color='leiden')
```
Tokens:
Raw