Gene expression dynamics after perturbation

[anshu957]

Thank you so much for the wonderful package!

I'm trying to calculate an expression of gene (over time) for perturbed vs un-perturbed case and then validate it using the time course data that I've for WT and KO.

Now, in order to get gene dynamics over time, I use kinetics_heatmap function and then use the following command:

gene1_time = adata.uns['kinetics_heatmap'].loc['gene1']

When I perform an in-silico perturbation using dyn.pd.KO(adata_wt,...) it returns the perturbed Vector field class. How do I proceed from there to get the gene expression over time?

[Xiaojieqiu]

Thanks for your interests to our work @anshu957 ! I will give you a more detailed answer for the dyn.pd.KO function. But you may want to look into this function: dyn.pd.perturbation. This function will return an updated adata with a new layer with the key j_delta_x_perturbation that stores the perturbation effects for each gene in each cell

[anshu957]

Thanks @Xiaojieqiu for your reply. To follow up, I did look up dyn.pd.perturbation but before I talk about that, let me briefly explain, how am I calculating gene expression dynamics for un-perturbed system. Since, I have access to time course data (longitudinal study) for my system, I can integrate the learned Vector Field in time and validate the expression dynamics of few key genes with the data just to ensure that I can get correct smooth gene expression dynamics (as a sanity check). Below is the code I use to do that.

# construct vf_pca object from adata
vf_pca = dyn.vf.SvcVectorField()
vf_pca.from_adata(adata_wt, basis='pca', vf_key='VecFld')

# get vector field function from vf_pca
vf_pca_func = vf_pca.func

# create a function to integrate the vector field
def vf_func_scipy(t, x):
    return vf_pca_func(x)

def integrate_traj(initial_state, t_end, t_step):
    

    t= np.linspace(0, t_end, int(t_end/t_step))
    # integrate the vector field
    sol = solve_ivp(vf_func_scipy, [0, t_end], initial_state[0,:], t_eval=t)
    return sol

# for all sampled cells, integrate the vector field

# sampled_cell_pca : (n_cells x pca_expression)

pca_mean = sampled_cells_pca.mean(axis=0).reshape(1, -1)

traj_list = []
for i in range(sampled_cells_pca.shape[0]):
    sol = integrate_traj(sampled_cells_pca[i,:].reshape(1, -1), 500, 10)

    # project the trajectory back to gene expression space
    traj_expr = pca_to_expr(sol.y.T, pcs, pca_mean.reshape(1,-1))

    traj_list.append(traj_expr)

This results in very good agreement for gene dynamics for my system :)

NOTE: I can also use dyn.pd.fate function to do the same but for my own understanding, I integrated myself. Moreover, I'm performing some more dynamical systems analysis on the vector field which is not present in dynamo at present.

Now I want to apply the perturbation and re-calculate the gene dynamics. From peeking into dyn.pd.perturbation function, I see that I can also set velocity of perturbed gene to be zero (which is what I want as the data I have is for a complete gene KO experiment). As per my understanding, this function returns a perturbed velocity vector, Jacobian for the perturbed gene and few other things. However, in order to integrate the perturbed system, I need to have a transformed vector field function. Do I need to re-learn the vector field using perturbed expression and velocities?