MultiVelo

https://github.com/welch-lab/multivelo
Admin
MultiVelo is a mechanistic model that extends RNA velocity by incorporating epigenomic data to infer...
Tokens:29,517
Snippets:234
Trust Score:8.2
License:BSD-3-Clause
Update:3 months ago
Show doc for...
Tokens:
Raw
### Install MultiVelo via pip

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

Use this command to install the package using the Python package manager.

```console
$ pip install multivelo
```

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

### Install MultiVelo via conda

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

Use this command to install the package from the bioconda channel.

```console
$ conda install -c bioconda multivelo
```

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

### Initialize AnnData Objects

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Setup AnnData objects for RNA and ATAC data with necessary layers and metadata.

```python
adata = AnnData(X = X)
adata.obs_names = obs_names
adata.var_names = var_names
adata.layers['unspliced'] = unspliced.copy()
adata.layers['spliced'] = spliced.copy()
adata.obs['true_time'] = true_time / 20

scv.pp.log1p(adata)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
adata.layers['Mu'] = unspliced.copy()
adata.layers['Ms'] = spliced.copy()
adata.var['highly_variable'] = True
```

```python
adata_atac = AnnData(X = atac.copy())
adata_atac.obs_names = obs_names
adata_atac.var_names = var_names

adata_atac.layers['Mc'] = adata_atac.X
adata_atac.obsp['distances'] = adata.obsp['distances']
adata_atac.obsp['connectivities'] = adata.obsp['connectivities']
```

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

### Initialize Environment and Imports

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/agg_fig_data.ipynb

Sets up the Python environment with necessary libraries for data analysis and visualization.

```python
%load_ext autoreload
%autoreload 2

import anndata as ad
import os
import scipy
import numpy as np
import pandas as pd
import sys
import multivelo as mv
import scanpy as sc
import scvelo as scv
import matplotlib.pyplot as plt

sys.path.append("/..")
```

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

### Initialize MultiVelo Environment

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig5.ipynb

Imports necessary libraries and appends the local path to sys.path for module resolution.

```python
import os
import scipy
import numpy as np
import pandas as pd
import scanpy as sc
import scvelo as scv
import sys
import multivelo as mv
import matplotlib.pyplot as plt
import requests
sys.path.append('../')
print(sys.path)
```

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

### Importing required libraries

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/neural_net_example.ipynb

Initializes the environment by importing necessary data science and bioinformatics libraries.

```python
import os
import scipy
import numpy as np
import pandas as pd
import sys
import multivelo as mv
import scanpy as sc
import scvelo as scv
import matplotlib.pyplot as plt
import requests
```

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

### Prepare Data Matrices

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

These snippets initialize observation/variable names and construct data matrices from noise-adjusted lists.

```python
obs_names = [str(i) for i in range(n)]
var_names = [str(i) for i in range(n_gene)]
```

```python
X = np.empty((n,n_gene))
unspliced = np.empty((n,n_gene))
spliced = np.empty((n,n_gene))
atac = np.empty((n,n_gene))
for i in range(n_gene):
    X[:,i] = noise_list[i][:,1] + noise_list[i][:,2]
    unspliced[:,i] = noise_list[i][:,1]
    spliced[:,i] = noise_list[i][:,2]
    atac[:,i] = noise_list[i][:,0]
```

```python
X.shape
```

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

### Directory Structure

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Demo.ipynb

The expected file and folder layout for the MultiVelo demo project.

```text
.
|-- MultiVelo_Demo.ipynb
|-- cell_annotations.tsv
|-- outs
|   |-- analysis
|   |   `-- feature_linkage
|   |       `-- feature_linkage.bedpe
|   |-- filtered_feature_bc_matrix
|   |   |-- barcodes.tsv.gz
|   |   |-- features.tsv.gz
|   |   `-- matrix.mtx.gz
|   `-- peak_annotation.tsv
|-- seurat_wnn
|   |-- nn_cells.txt
|   |-- nn_dist.txt
|   `-- nn_idx.txt
`-- velocyto
    `-- 10X_multiome_mouse_brain.loom
```

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

### Configure Paths and Load Modules

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Sets up the autoreload extension and appends the source directory to the system path for importing custom modules.

```python
%load_ext autoreload
%autoreload 2
sys.path.append("../src/multivelo")
from dynamical_chrom_func import *
import multivelo as mv
```

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

### Setting global parameters

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/neural_net_example.ipynb

Configures verbosity, display options, and printing settings 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)
```

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

### Configure Visualization Data

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/train_neural_net.ipynb

Sets up paths and parameters for loading noiseless data for performance visualization.

```python
read_folder = "./data/simulated_data/" + suffix_name + "_noiseless"
```

```python
genes_to_graph = [1,12]
```

```python
graph_X, graph_t, subset_gene = X_from_file(read_folder, dire, subset_gene=genes_to_graph)
```

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

### Initialize Simulation Data Structures

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Initializes lists and matrices for storing simulation outputs and state variables.

```python
alpha_c = []
alpha = []
beta = []
gamma = []
output_time = []
tau_matrix = []
c0_matrix = []
u0_matrix = []
s0_matrix = []
```

```python
adj_switch_times = []
raw_switch_times = []
```

```python
# np.random.seed(0)
cus_list, cus_sw_list, cus_ss_list, noise_list, para_list = [], [], [], [], []
var_sim_c, var_sim_u, var_sim_s = [], [], []
mean_sim_c, mean_sim_u, mean_sim_s = [], [], []
```

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

### Import Dependencies

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Template.ipynb

Required libraries for the MultiVelo workflow.

```python
import os
import scipy
import numpy as np
import pandas as pd
import scanpy as sc
import scvelo as scv
import multivelo as mv
import matplotlib.pyplot as plt
```

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

### Import Dependencies and Configure Autoreload

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig2.ipynb

Initializes the environment by loading autoreload extensions and importing necessary data science and bioinformatics libraries.

```python
%load_ext autoreload
%autoreload 2

import os
import scipy
import numpy as np
import pandas as pd
import sys
import multivelo as mv
import scanpy as sc
import scvelo as scv
import matplotlib.pyplot as plt

sys.path.append("/..")
```

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

### Directory Structure Visualization

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig3.ipynb

Displays the expected file and folder hierarchy for the MultiVelo analysis project.

```text
.
|-- MultiVelo_Fig3.ipynb
|-- cell_annotations.tsv
|-- outs
|   |-- analysis
|   |   `-- feature_linkage
|   |       `-- feature_linkage.bedpe
|   |-- filtered_feature_bc_matrix
|   |   |-- barcodes.tsv.gz
|   |   |-- features.tsv.gz
|   |   `-- matrix.mtx.gz
|   `-- peak_annotation.tsv
|-- seurat_wnn
|   |-- nn_cells.txt
|   |-- nn_dist.txt
|   `-- nn_idx.txt
`-- velocyto
    `-- 10X_multiome_mouse_brain.loom
```

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

### Run scVelo Workflow

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Standard scVelo pipeline for recovering dynamics, calculating velocity, and generating embeddings.

```python
# scvelo
```

```python
scv.tl.recover_dynamics(adata, n_jobs=12)
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)
```

```python
scv.pl.velocity_embedding_stream(adata, basis='umap')
```

```python
scv.tl.latent_time(adata)
scv.pl.scatter(adata, color='latent_time', color_map='gnuplot', size=80)
```

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

### Load Training Tools

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/train_neural_net.ipynb

Enables autoreload and imports custom training utility functions.

```python
%load_ext autoreload
%autoreload 2

# import function that makes our data
from training_tools import *
```

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

### Initialize Model Weights

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/train_neural_net.ipynb

Applies specific weight initialization strategies based on the model direction.

```python
# Used this: https://stackoverflow.com/questions/49433936/how-do-i-initialize-weights-in-pytorch

def init_weights(m):
    if dire == 0:
        if isinstance(m, nn.Linear):nn.init.xavier_uniform_(m.weight)
    else:
        if isinstance(m, nn.Linear):nn.init.kaiming_normal_(m.weight, nonlinearity='relu')
```

```python
ode_model.apply(init_weights)
```

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

### Directory Structure

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Template.ipynb

Expected file structure for the MultiVelo analysis pipeline.

```text
. 
|-- MultiVelo_Template.ipynb 
|-- outs 
|   |-- analysis 
|   |   `-- feature_linkage 
|   |       `-- feature_linkage.bedpe 
|   |-- filtered_feature_bc_matrix 
|   |   |-- barcodes.tsv.gz 
|   |   |-- features.tsv.gz 
|   |   `-- matrix.mtx.gz 
|   `-- atac_peak_annotation.tsv 
|-- seurat_wnn 
|   |-- nn_cells.txt 
|   |-- nn_dist.txt 
|   `-- nn_idx.txt 
`-- velocyto 
    `-- gex_possorted_bam_XXXXX.loom
```

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

### Initialize Rate Parameter Lists

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Initializes empty lists to store gene rate parameters.

```python
alpha_c = []
alpha = []
beta = []
gamma = []
```

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

### Directory Structure Overview

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig2.ipynb

Visual representation of the required file and folder structure for the MultiVelo Fig2 analysis.

```text
.
|-- MultiVelo_Fig2.ipynb
|-- cell_annotations.tsv
|-- outs
|   |-- analysis
|   |   `-- feature_linkage
|   |       `-- feature_linkage.bedpe
|   |-- filtered_feature_bc_matrix
|   |   |-- barcodes.tsv.gz
|   |   |-- features.tsv.gz
|   |   `-- matrix.mtx.gz
|   `-- peak_annotation.tsv
|-- seurat_wnn
|   |-- nn_cells.txt
|   |-- nn_dist.txt
|   `-- nn_idx.txt
`-- velocyto
    `-- 10X_multiome_mouse_brain.loom
```

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

### Run MultiVelo Dynamics

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Template.ipynb

Execute the main MultiVelo function to recover dynamics.

```python
# Run MultiVelo main function
adata_result = mv.recover_dynamics_chrom(adata_rna,
                                         adata_atac,
                                         max_iter=5, # coordinate-descent like optimization
                                         init_mode="invert", # simple, invert, or grid
                                         verbose=False,
                                         parallel=True,
                                         save_plot=False,
                                         rna_only=False,
                                         fit=True,
                                         n_anchors=500,
                                         extra_color_key='leiden' # used if save_plot=True
                                        )
# Full argument list can be shown with help(mv.recover_dynamics_chrom)
```

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

### Set Simulation Output Directory

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Defines the destination folder where the generated simulated data will be saved.

```python
# name of the folder we save the simulated data to

save_folder = "./data/simulated_data/dir2_m1"
```

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

### Loading datasets

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/neural_net_example.ipynb

Reads the AnnData objects from local paths or downloads them from the provided URLs if missing.

```python
adata_rna = sc.read(rna_path, backup_url=rna_url)
adata_atac = sc.read(atac_path, backup_url=atac_url)
```

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

### Visualize model summaries

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Demo.ipynb

Generates summary plots for the dynamical model results, including pie charts, switch times, and likelihood distributions.

```python
mv.pie_summary(adata_result)
```

```python
mv.switch_time_summary(adata_result)
```

```python
mv.likelihood_plot(adata_result)
```

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

### Save Initial State Matrices

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Saves initial state matrices (s0, u0, c0) to text files after verifying their shapes.

```python
s0_matrix = np.array(s0_matrix)
u0_matrix = np.array(u0_matrix)
c0_matrix = np.array(c0_matrix)

print(s0_matrix.shape)
print(u0_matrix.shape)
print(c0_matrix.shape)

s0_name = save_folder + "/s0_data.txt"
u0_name = save_folder + "/u0_data.txt"
c0_name = save_folder + "/c0_data.txt"

np.savetxt(s0_name, s0_matrix)
np.savetxt(u0_name, u0_matrix)
np.savetxt(c0_name, c0_matrix)
```

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

### Configure Training Parameters

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/train_neural_net.ipynb

Sets the adjustment factor and directionality flags for the model training process.

```python
adj_factor = 21
```

```python
# direction of the data we're fitting
dire = 2

# model 1 repression genes are notoriously hard to train,
# raise this flag if training them and some extra help
# will be added
dir_2_model_1 = True
```

```python
if dir_2_model_1 and not dire == 2:
    raise Exception("Don't set dir_2_model_1 to True if you aren't training repression genes!")
```

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

### Generate Summary Plots

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Template.ipynb

Visualize gene type distributions and switch time statistics.

```python
mv.pie_summary(adata_result) # gene type chart
```

```python
mv.switch_time_summary(adata_result) # switch time statistics
```

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

### Pie Summary

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Generate a pie summary of the results.

```python
pie_summary(adata_result2)
```

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

### Import MultiVelo

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

Standard import statement for the multivelo package.

```default
import multivelo as mv
```

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

### Configure MultiVelo Logging

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig2.ipynb

Sets the log filename for MultiVelo using the current system time to ensure unique log files.

```python
import time

mv.settings.LOG_FILENAME = "Fig2_" + str(time.time()) + ".txt"
```

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

### Generate Switch Time Values

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Calculates switch times by generating random intervals and scaling them to a fixed total time.

```python
# Generate switch time values

np.random.seed(0)
t_sw_array_array = np.empty((n_gene, 3))
for i in range(n_gene):

    # generate 4 intervals drawn from equally wide distributions
    interval1 = np.random.choice(list(range(1,10)))
    interval2 = np.random.choice(list(range(1,10)))
    interval3 = np.random.choice(list(range(1,10)))
    interval4 = np.random.choice(list(range(1,10)))

    # scale the intervals so that the time total is equal to 20
    total = interval1 + interval2 + interval3 + interval4
    scale = 20 / total
    interval1 *= scale
    interval2 *= scale
    interval3 *= scale
    interval4 *= scale

    # calculate the switch times from those intervals
    t_sw_array = np.array([interval1, interval1+interval2, interval1+interval2+interval3])
    t_sw_array_array[i,:] = t_sw_array
        
t_sw_array_array
```

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

### Print AnnData Summary

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig4.ipynb

Display the summary of the resulting AnnData object.

```python
print(adata_result)
```

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

### Configuring logging

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/neural_net_example.ipynb

Sets the log filename using the current timestamp.

```python
import time

mv.settings.LOG_FILENAME = "nn_example_" + str(time.time()) + ".txt"
```

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

### Running standard dynamical model

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/neural_net_example.ipynb

Performs standard RNA velocity recovery using scVelo's dynamical mode.

```python
scv.tl.recover_dynamics(adata_rna)
scv.tl.velocity(adata_rna, mode="dynamical")
scv.tl.velocity_graph(adata_rna, n_jobs=1)
```

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

### Generate pie summary

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig6.ipynb

Creates a pie chart summary of the velocity results.

```python
mv.pie_summary(adata_result)
```

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

### Visualize Rate Parameters

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/agg_fig_data.ipynb

Generates a grid of 2D histograms comparing different rate parameters.

```python
# graph a set of rate parameters
def graph_params(alpha_c, alpha, beta, gamma):

    # the names of the parameters for axis labels
    names = ["alpha_c", "alpha", "beta", "gamma"]

    # create a subplot
    fig, axs = plt.subplots(4, 3, figsize=(12, 12))

    # axis boundaries for each rate parameter
    lims = [0.2, 2, 1.25, 2]

    # assemble the rate parameters into a matrix
    mvln = np.array([alpha_c, alpha, beta, gamma])

    # for each combination of rate parameters...
    for i in range(4):
        col = 0
        for j in range(4):

            # if we're about to plot a rate parameter against
            # itself, then skip this iteration
            if i == j:
                continue

            # graph the rate parameters
            h = axs[i][col].hist2d(np.ravel(mvln[i,:]), np.ravel(mvln[j,:]),
                                   # range=[[0, lims[i]],
                                   #        [0, lims[j]]],
                                   bins=10,
                                   cmap="Greens")
            axs[i][col].set_xlabel(names[i])
            axs[i][col].set_ylabel(names[j])
            fig.colorbar(h[3], ax=axs[i][col])
            col += 1
    
    # # fig.colorbar(axs)
    fig.tight_layout()
```

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

### Calculate Rate Ratios

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Compute and print alpha/beta and alpha/gamma ratios.

```python
print("alpha/beta:")
print(adata_bad.var['fit_alpha'] / adata_bad.var['fit_beta'])
print()
print("alpha/gamma:")
print(adata_bad.var['fit_alpha'] / adata_bad.var['fit_gamma'])
```

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

### Export and Visualize scVelo Metrics

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Exports scVelo results to file and visualizes latent time differences and variance.

```python
a = adata.obs['latent_time'] - adata.obs['true_time']
b = np.var(adata.layers['fit_t']/20, axis=1)
df = pd.DataFrame({'diff':a, 'var':b})
df.to_csv('sim_4-5-1_scvelo.txt', sep='\t', index=False)
```

```python
scv.pl.scatter(adata, color=adata.obs['latent_time'] - adata.obs['true_time'], color_map='gnuplot', size=80)
```

```python
scv.pl.scatter(adata, color=np.var(adata.layers['fit_t']/20, axis=1))
```

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

### Train Neural Network Model

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/train_neural_net.ipynb

Configures hyperparameters, initializes progress bars, and executes the training loop with batch processing and validation checks.

```python
# Number of max epochs for each neural network.
# Each N value was the epoch at which the network ceased
# to improve performance on the developer's validation data
# set for five epochs.
if dire == 0:
    N = 51
elif dire == 1:
    N = 67
elif dire == 2:
    if dir_2_model_1:
        N = 43
    else:
        N = 51

# Use MSE as our loss criterion
criterion = torch.nn.MSELoss(reduction="mean")
val_criterion = torch.nn.MSELoss(reduction="mean")

# learning rate for adam
if dir_2_model_1:
# if dir_2_model_1:
    adam_lr = 1e-3
elif dire == 0:
    adam_lr = 3e-4
else:
    adam_lr = 1e-4

# keep track of the minimum validation loss
min_loss = float('inf')

epoch_loss = []
val_loss_list = []
val_epoch = []

# progress bar keeping track of 
# how many training epochs have
# transpired
epoch_bar = tqdm(total=N)

# progress bar keeping track of
# how many batches within the
# epoch have transpired
batch_bar = tqdm(total=batches)

# mean training loss over each batch
# of the epoch
mean_loss = 0

# set Adam as the training optimizer
optimizer = torch.optim.Adam(ode_model.parameters(), lr=adam_lr)

# print out some hyperparameters
print("Running with lr =", adam_lr, "and base_n =", base_n)

ode_model.train()

# loop thru each epoch
for i in range(N):

    # reset the epoch mean loss
    mean_loss = 0
    mean_context_loss = 0
    
    batch_bar.reset()
    
    sleep(0.001)

    # randomly select data points and assign them to batches
    x_tens, t_tens = make_batches_random(batches, X, t)

    # loop thru each batch
    for j in range(batches):
        
        ode_model.float()
        
        torch.autograd.set_detect_anomaly(True)
        
        # do a forward pass
        t_pred = ode_model(x_tens[j])

        # for loss values, set t_pred and t to 
        # the original time scale
        scaled_t_pred = torch.mul(t_pred, adj_factor)
        scaled_t = torch.mul(t_tens[j], adj_factor)
        
        # check loss
        loss = criterion(scaled_t, scaled_t_pred)
    
        # reset gradients
        optimizer.zero_grad()
    
        # do a backward pass
        loss.backward()

        # step throught the optimizer
        optimizer.step()
        
        batch_bar.update()

        # add the loss of this batch
        # to the running mean batch loss
        mean_loss += loss.item()

    # calculate the mean batch loss
    mean_loss /= batches

    # add loss values to lists for graphing later
    epoch_loss.append(mean_loss)

    # periodically check our progress
    # (can adjust value after modulo in case
    # you want to check progress less regularly)
    if i % 1 == 0:

        ode_model.eval()
        
        # do a forward pass on validation data
        val_t_pred = ode_model(val_X_ten)
        
        # check loss on validation data
        val_loss = val_criterion(torch.mul(val_t_pred, adj_factor), torch.mul(val_t_ten, adj_factor))
        val_loss_list.append(val_loss.item())
        val_epoch.append(i)
        
        print("Epoch", i, "loss:", mean_loss)
        print("Validation loss:", val_loss.item())
    
        # find the relative increase between the two previous
        # loss values for our validation set
        rel_diff = 0
    
        new_val_N = len(val_loss_list)

        # as long as we have at least two items to compare,
        # compare the relative difference in successive
        # validation loss values
        if new_val_N >= 2:
            rel_diff = (val_loss_list[-1] - min_loss) / min_loss
            print(str(rel_diff*100) + "%")

        # # keep track of the lowest validation loss    
        if val_loss_list[-1] < min_loss:
            min_loss = val_loss_list[-1]
            print("New min!")

        print()
        
        ode_model.train()
        
    epoch_bar.update()
    
ode_model.eval()
```

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

### Define Data Directories

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Specifies the folder path for reading original rate parameters used as a basis for simulation.

```python
# the folder from which to read the original
# rate parameters from

# figure data:
read_folder = "./data/rates/fig_rates_train"

# validation data:
# You will need to supply your own validation data in order to run this code!
# read folder =
```

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

### Visualize Distribution Comparisons

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Histograms comparing log-transformed means and variances of observed versus simulated data.

```python
plt.hist(np.log(mean_c), bins=100, alpha=0.5);
plt.hist(np.log(mean_sim_c), bins=100, alpha=0.5);
```

```python
plt.hist(np.log(mean_u), bins=100, alpha=0.5);
plt.hist(np.log(mean_sim_u), bins=100, alpha=0.5);
```

```python
plt.hist(np.log(mean_s), bins=100, alpha=0.5);
plt.hist(np.log(mean_sim_s), bins=100, alpha=0.5);
```

```python
plt.hist(np.log(var_c / mean_c), bins=100, alpha=0.5);
plt.hist(np.log(var_sim_c / mean_sim_c), bins=100, alpha=0.5);
```

```python
plt.hist(np.log(var_u / mean_u), bins=100, alpha=0.5);
plt.hist(np.log(var_sim_u / mean_sim_u), bins=100, alpha=0.5);
```

```python
plt.hist(np.log(var_s / mean_s), bins=100, alpha=0.5);
plt.hist(np.log(var_sim_s / mean_sim_s), bins=100, alpha=0.5);
```

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

### Plot Model Statistics

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Template.ipynb

Visualize likelihood and model parameter statistics.

```python
mv.likelihood_plot(adata_result) # likelihood and model parameter statistics
```

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

### Define Data Paths and Load Data

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/train_neural_net.ipynb

Determines the directory path based on the direction flag and loads the training data.

```python
if dire == 0:
    suffix_name = "dir0"
elif dire == 1:
    suffix_name = "dir1"
elif dire == 2:
    if dir_2_model_1:
        suffix_name = "dir2_m1"
    else:
        suffix_name = "dir2_m2"
```

```python
read_folder = "./data/simulated_data/" + suffix_name
```

```python
X, t = X_from_file(read_folder, dire)
```

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

### Summarize Switch Time

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig3.ipynb

Generates a summary plot for switch times.

```python
mv.switch_time_summary(adata_result)
```

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

### Executing neural network dynamical model

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/neural_net_example.ipynb

Runs the MultiVelo dynamical recovery with neural network time prediction enabled. High parallelization is recommended for this step.

```python
# This will take a while. Parallelization is high recommended.

mv.settings.VERBOSITY = 0

adata_result = mv.recover_dynamics_chrom(adata_rna, 
                                        adata_atac,
                                        max_iter=5, 
                                        init_mode="invert",
                                        parallel=True,
                                        n_jobs = 10,
                                        save_plot=False,
                                        rna_only=False,
                                        fit=True,
                                        n_anchors=500,
                                        neural_net=True,
                                        # device="cuda:0",
                                        extra_color_key='cluster'
                                        )
```

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

### Compare Models and Directions

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Calculate matches between model predictions and directions.

```python
np.sum(model_list[filt] == model_pred[filt])
```

```python
np.sum(dire == dire_pred)
```

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

### Load MultiVelo Results

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig6.ipynb

Loads the MultiVelo result file, with an automatic fallback to download from a Figshare URL if the local file is missing.

```python
h5ad_url = "https://figshare.com/ndownloader/files/40064449"
adata_result = sc.read("multivelo_result_fig6.h5ad", backup_url = h5ad_url)
```

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

### Import Required Libraries

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig3.ipynb

Imports essential data science and bioinformatics libraries, including the development version of MultiVelo.

```python
import os
import scipy
import numpy as np
import pandas as pd
import scanpy as sc
import scvelo as scv
import sys
import multivelo as mv
import matplotlib.pyplot as plt
```

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

### Load WNN Neighbors

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Template.ipynb

Load neighbor indices and distances from Seurat WNN output.

```python
# Back in python, load the neighbors
nn_idx = np.loadtxt("seurat_wnn/nn_idx.txt", delimiter=',')
nn_dist = np.loadtxt("seurat_wnn/nn_dist.txt", delimiter=',')
nn_cells = pd.Index(pd.read_csv("seurat_wnn/nn_cells.txt", header=None)[0])
```

```python
np.all(nn_cells == adata_atac.obs_names) # make sure cell names match
```

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

### Load simulation data from CSV files

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Reads multiple tab-separated simulation result files into pandas DataFrames.

```python
a = pd.read_csv('sim_4-5-1.txt', sep='\t')
b = pd.read_csv('sim_4-4-2.txt', sep='\t')
c = pd.read_csv('sim_4-3-3.txt', sep='\t')
d = pd.read_csv('sim_6-3-1.txt', sep='\t')
e = pd.read_csv('sim_1-8-1.txt', sep='\t')
f = pd.read_csv('sim_4-5-1_low_var.txt', sep='\t')
g = pd.read_csv('sim_4-5-1_high_var.txt', sep='\t')

h = pd.read_csv('sim_4-5-1_scvelo.txt', sep='\t')
i = pd.read_csv('sim_4-4-2_scvelo.txt', sep='\t')
j = pd.read_csv('sim_4-3-3_scvelo.txt', sep='\t')
k = pd.read_csv('sim_6-3-1_scvelo.txt', sep='\t')
l = pd.read_csv('sim_1-8-1_scvelo.txt', sep='\t')
m = pd.read_csv('sim_4-5-1_low_var_scvelo.txt', sep='\t')
n = pd.read_csv('sim_4-5-1_high_var_scvelo.txt', sep='\t')
```

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

### Import Dependencies

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/train_neural_net.ipynb

Standard library and scientific computing imports required for data manipulation and neural network training.

```python
# %load_ext autoreload
# %autoreload 2

import os
import scipy
import numpy as np
import pandas as pd
import sys
import scanpy as sc
import scvelo as scv
import matplotlib.pyplot as plt
import math
import random
import time
```

```python
# Use this for progress bar help: https://stackoverflow.com/questions/3160699/python-progress-bar
from time import sleep
from tqdm.notebook import tqdm
```

```python
#import pytorch
import torch
from torch import nn
```

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

### Configure Outlier Handling

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Defines scaling factors for handling larger parameter values during data generation.

```python
# some gene values are much larger than the majority,
# but aren't large enough to be discarded as outliers.
# When this happens, these two variables help govern how
# to generate extra data around these points

# Governs how much we want the gene multiplier
# to grow for larger parameter values
multiplier = 10
# multiplier = 0

# Governs how much we
# want the radius to increase
# for these larger values
radius_scale = 1
```

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

### Print Likelihood Statistics

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Print length, max, and min of fit likelihood.

```python
print(len(adata_result2.var['fit_likelihood']))
print(max(adata_result2.var['fit_likelihood']))
print(min(adata_result2.var['fit_likelihood']))
```

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

### Define Simulation Parameters

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Sets core simulation variables including cell count, gene directionality, and noise factors.

```python
# number of cells to simulate per gene
n = 750

# the direction of the gene we're simulating
# 0 = complete genes
# 1 = activation genes
# 2 = repression genes
dire = 2

# the list of appropriate gene models to simulate
models = [1]

# the factor to multiply the noise for c/u/s by
factor = 1
# factor = 0
```

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

### Read and Format RNA Data

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig3.ipynb

Load a loom file and ensure observation names are unique and correctly formatted.

```python
adata_rna = scv.read("velocyto/10X_multiome_mouse_brain.loom", cache=True)
adata_rna.obs_names = [x.split(':')[1][:-1] + '-1' for x in adata_rna.obs_names]
adata_rna.var_names_make_unique()
```

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

### Export Rate Parameters

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/agg_fig_data.ipynb

Aggregates rate parameters from multiple AnnData objects and saves them as text files.

```python
# function responsible for saving the relevant data
def write_files(outfile, adatas, graph=False):

        # print the name of the file we're saving to
        print(outfile)

        alpha_c = np.array([])
        alpha = np.array([])
        beta = np.array([])
        gamma = np.array([])

        # assemble the rate parameters of each AnnData object
        for adata in adatas:
        
            alpha_c = np.concatenate((alpha_c, np.array(adata.var['fit_alpha_c'])))
            alpha = np.concatenate((alpha, np.array(adata.var['fit_alpha'])))
            beta = np.concatenate((beta, np.array(adata.var['fit_beta'])))
            gamma = np.concatenate((gamma, np.array(adata.var['fit_gamma'])))

        # graph results if the user specifies it
        if graph:
            graph_results(alpha_c, alpha, beta, gamma)

        # save all of the data
        np.savetxt(outfile + "/alpha_c.txt", alpha_c)
        np.savetxt(outfile + "/alpha.txt", alpha)
        np.savetxt(outfile + "/beta.txt", beta)
        np.savetxt(outfile + "/gamma.txt", gamma)

        # print the total number of genes saved
        print("Number of genes:", alpha_c.shape[0])
```

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

### Load and Sort Rate Parameters

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Loads original rate parameters from text files and sorts them for processing.

```python
ac_orig = np.loadtxt(read_folder+"/alpha_c.txt")
a_orig = np.loadtxt(read_folder+"/alpha.txt")
b_orig = np.loadtxt(read_folder+"/beta.txt")
g_orig = np.loadtxt(read_folder+"/gamma.txt")

big_a_orig = 10
big_ac_orig = 0.1
big_b_orig = 1
big_g_orig = 1.5

print(big_a_orig)

ac_too_big = 0.2
a_too_big = 100
b_too_big = 1.75
g_too_big = 2.2

# m_orig = np.loadtxt(read_folder+"/model_list.txt")

ac_sorted = np.sort(ac_orig)
a_sorted = np.sort(a_orig)
b_sorted = np.sort(b_orig)
g_sorted = np.sort(g_orig)

orig_n_gene = len(ac_orig)
```

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

### Load MultiVelo Result

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig4.ipynb

Load the processed AnnData object, optionally downloading it from a remote URL if the local file is missing.

```python
h5ad_url = "https://figshare.com/ndownloader/files/40064263"
adata_result = sc.read("multivelo_result_fig4.h5ad", backup_url = h5ad_url)
```

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

### Annotate Clusters

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Template.ipynb

Rename and visualize clusters.

```python
# Identify cell types
new_cluster_names = [...]
```

```python
adata_rna.rename_categories('leiden', new_cluster_names) # annotate clusters
```

```python
scv.pl.umap(adata_rna, color='leiden')
```

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

### Visualize Training Data and Model Predictions

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/train_neural_net.ipynb

Plots real versus predicted time and gene expression dynamics using matplotlib. Requires pre-computed graph_X, graph_t, and a trained ode_model.

```python
f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(16, 16))

perf_line = np.arange(0, 20)
ax1.plot(perf_line, perf_line, label="perfect fit", linewidth=1)
ax1.set_xlabel("Real")
ax1.set_ylabel("Pred")

ax2.set_ylabel("Time")
ax2.set_xlabel("C")

ax3.set_ylabel("Time")
ax3.set_xlabel("U")

ax4.set_ylabel("Time")
ax4.set_xlabel("S")

epsilon = 1e-5

for i in range(len(genes_to_graph)):

    subset_X = graph_X[subset_gene[i]:subset_gene[i+1]]
    subset_t = graph_t[subset_gene[i]:subset_gene[i+1]]

    alpha_c = str(np.round(np.exp(subset_X[0,3]) - epsilon, 4))
    alpha = str(np.round(np.exp(subset_X[0,4]) - epsilon, 4))
    beta = str(np.round(np.exp(subset_X[0,5]) - epsilon, 4))
    gamma = str(np.round(np.exp(subset_X[0,6]) - epsilon, 4))

    c = subset_X[:,0]
    u = np.exp(subset_X[:,1] - epsilon)
    s = np.exp(subset_X[:,2] - epsilon)

    t_pred = ode_model(torch.tensor(subset_X).reshape(-1, subset_X.shape[1]))
    t_pred = (t_pred.detach().numpy().reshape(-1) * adj_factor) - 1

    ax1.plot(subset_t, t_pred, linewidth=2, label="alpha_c: " + alpha_c)

    clabel = " t - alpha_c: " + alpha_c
    ax2.plot(c, subset_t, label="real" + clabel, linewidth=1)
    ax2.plot(c, t_pred, label="pred" + clabel, linewidth=1)

    ulabel = " t - alpha: " + alpha + " and beta: " + beta
    ax3.plot(u, subset_t, label="real" + ulabel, linewidth=1)
    ax3.plot(u, t_pred, label="pred" + ulabel, linewidth=1)

    slabel = " t - beta: " + beta + " and gamma: " + gamma
    ax4.plot(s, subset_t, label="real" + slabel, linewidth=1)
    ax4.plot(s, t_pred, label="pred" + slabel, linewidth=1)
    

ax1.legend()
ax2.legend()
ax3.legend()
ax4.legend()
```

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

### Recover Dynamics

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Run MultiVelo dynamic recovery on RNA and ATAC data.

```python
mv.settings.VERBOSITY = 1
mv.settings.LOG_FILENAME = "simdata_test.txt"

adata_result2 = mv.recover_dynamics_chrom(adata, 
                                       adata_atac, 
                                       max_iter=5, 
                                       init_mode="invert",
                                       plot=False, 
                                       parallel=True,
                                       n_jobs=10,
                                       save_plot=False,
                                       rna_only=False,
                                       fit=True
                                       )
```

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

### Save and load MultiVelo results

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig2.ipynb

Persists the resulting AnnData object to disk and reloads it for subsequent analysis.

```python
# Save the result for use later on
adata_result.write("multivelo_result_fig2.h5ad")
```

```python
adata_result = sc.read_h5ad("multivelo_result_fig2.h5ad")
```

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

### Run MultiVelo dynamical model

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig3.ipynb

Executes the MultiVelo dynamical model incorporating chromatin accessibility data. High parallelization is recommended for performance.

```python
# This will take a while. Parallelization is high recommended.

mv.settings.VERBOSITY = 0

adata_result = mv.recover_dynamics_chrom(adata_rna, 
                                        adata_atac,
                                        max_iter=5, 
                                        init_mode="invert",
                                        parallel=True,
                                        n_jobs = 10,
                                        save_plot=False,
                                        rna_only=False,
                                        fit=True,
                                        n_anchors=500, 
                                        extra_color_key='celltype'
                                        )
```

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

### Load MultiVelo Result

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Template.ipynb

Read the processed h5ad file containing MultiVelo results.

```python
# adata_result = sc.read_h5ad('multivelo_result.h5ad')
```

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

### Plot Summary Statistics

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig5.ipynb

Displays summary plots for pie charts and switch times.

```python
mv.pie_summary(adata_result)
```

```python
mv.switch_time_summary(adata_result)
```

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

### Visualize Gene Expression Dynamics

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Generates plots for gene expression over time, including switch time markers.

```python
graph_genes = 500

if graph_genes > n_gene:
    graph_genes = n_gene

fig, axs = plt.subplots(-(-graph_genes // 5), 5, figsize=(14, 2.5 * (-(-graph_genes // 5))), facecolor='white')
fig.patch.set_facecolor('white')
marker_type = ['o', 'X', 'D', '*']
for i in range(graph_genes):
    c = cus_list[i][:,0]
    u = cus_list[i][:,1]
    s = cus_list[i][:,2]
    t = (output_time.T)[i]
    row = i // 5
    col = i % 5
    axs[row, col].scatter(t, c, s=5)
    # only one line may be specified; full height
    #print(t_sw_array_array[i])
    # for line in t_sw_array_array[i]:
    for line in adj_switch_times[i]:
        axs[row, col].axvline(x = line, color = 'r')
    #for j, sw_exp in enumerate(cus_sw_list[i]):
    #    axs[row, col].plot([output_time.T)], [sw_exp[0,0]], marker_type[j]+'r', markersize=8)
    axs[row, col].plot(t[-1], c[-1], marker_type[3]+'r', markersize=8)
    axs[row, col].set_xlabel('t')
    axs[row, col].set_ylabel('c')
    axs[row, col].set_title(f'S{i} M{int(model_list[i])} D{dire}')
fig.tight_layout()
```

```python
graph_genes = 500

if graph_genes > n_gene:
    graph_genes = n_gene

fig, axs = plt.subplots(-(-graph_genes // 5), 5, figsize=(14, 2.5 * (-(-graph_genes // 5))), facecolor='white')
fig.patch.set_facecolor('white')
marker_type = ['o', 'X', 'D', '*']
for i in range(graph_genes):
    c = cus_list[i][:,0]
    u = cus_list[i][:,1]
    s = cus_list[i][:,2]
    t = (output_time.T)[i]
    row = i // 5
    col = i % 5
    axs[row, col].scatter(t, u, s=5)
    # only one line may be specified; full height
    #print(t_sw_array_array[i])
    #for line in t_sw_array_array[i]:
    for line in adj_switch_times[i]:
        axs[row, col].axvline(x = line, color = 'r')
    #for j, sw_exp in enumerate(cus_sw_list[i]):
    #    axs[row, col].plot([output_time.T)], [sw_exp[0,0]], marker_type[j]+'r', markersize=8)
    axs[row, col].plot(t[-1], u[-1], marker_type[3]+'r', markersize=8)
    axs[row, col].set_xlabel('t')
    axs[row, col].set_ylabel('u')
    axs[row, col].set_title(f'S{i} M{int(model_list[i])} D{dire}')
fig.tight_layout()
```

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

### Visualize Scatter Plots

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Generate scatter plots for gene expression or metadata.

```python
scv.pl.scatter(adata, adata.var_names[:10])
```

```python
scv.pl.scatter(adata, color='true_time', color_map='gnuplot', size=80)
```

```python
scv.pl.scatter(adata_result2, color='latent_time', color_map='gnuplot', size=80)
```

```python
scv.pl.scatter(adata_result2, color=['root_cells'])
```

```python
scv.pl.scatter(adata_result2, color=adata_result2.obs['latent_time'] - adata_result2.obs['true_time'], color_map='gnuplot', size=80)
```

```python
scv.pl.scatter(adata_result2, color=np.var(adata_result2.layers['fit_t']/20, axis=1))
```

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

### Perform Leiden Clustering

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Fig3.ipynb

Runs the Leiden clustering algorithm on the AnnData object.

```python
sc.tl.leiden(adata_result)
```

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

### Read H5AD File

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Load an AnnData object from an H5AD file.

```python
adata_result = sc.read_h5ad('../10x_mouse_brain/multivelo.h5ad')
```

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

### Set Random Seeds

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/train_neural_net.ipynb

Ensures reproducibility by setting seeds for both PyTorch and NumPy.

```python
torch.manual_seed(42)
np.random.seed(42)
```

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

### Load and Filter RNA Data

Source: https://github.com/welch-lab/multivelo/blob/main/Examples/MultiVelo_Template.ipynb

Read the loom file and apply basic cell filtering.

```python
# Read in RNA and filter
adata_rna = scv.read('velocyto/gex_possorted_bam_XXXXX.loom', cache=True)
adata_rna.obs_names = [x.split(':')[1][:-1] + '-1' for x in adata_rna.obs_names]
adata_rna.var_names_make_unique()
sc.pp.filter_cells(adata_rna, min_counts=...)
sc.pp.filter_cells(adata_rna, max_counts=...)
```

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

### Load AnnData Objects

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/agg_fig_data.ipynb

Loads pre-computed AnnData objects from H5AD files for training or validation.

```python
# training data

# data used to generate figures in paper
fig3 = sc.read_h5ad("../Examples/multivelo_result_fig3.h5ad")
fig4 = sc.read_h5ad("../Examples/multivelo_result_fig4.h5ad")
fig5 = sc.read_h5ad("../Examples/multivelo_result_fig5.h5ad")
fig6 = sc.read_h5ad("../Examples/multivelo_result_fig6.h5ad")

figs = [fig3, fig4, fig5, fig6]
```

```python
# vali data
# You will need to supply your own validation data in order to run this code!
# val_data = sc.read_h5ad()
# figs = [val_data]
```

```python
N = len(figs)
```

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

### Save Switch Times

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Converts switch times to an array and saves them to a text file.

```python
adj_switch_times = np.array(adj_switch_times)

sw_name = save_folder + "/switch_times.txt"

np.savetxt(sw_name, adj_switch_times)
```

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

### Save Data

Source: https://github.com/welch-lab/multivelo/blob/main/train_neural_network/sim_data_from_rates.ipynb

Write the result object to an H5AD file.

```python
adata_result2.write("sim_data.h5ad")
```
Tokens:
Raw