aristoteleo / dynamo-release

Inclusive model of expression dynamics with conventional or metabolic labeling based scRNA-seq / multiomics, vector field reconstruction and differential geometry analyses
https://dynamo-release.readthedocs.io/en/latest/
BSD 3-Clause "New" or "Revised" License
415 stars 59 forks source link

GATA2-PLEK pair cannot be found in 'Minimal network for basophil lineage commitment' tutorial #375

Closed hyjforesight closed 1 year ago

hyjforesight commented 2 years ago

Hello Dynamo, When I was running this tutorial (https://dynamo-release.readthedocs.io/en/latest/notebooks/tutorial_hsc_dynamo_basophil_lineage/tutorial_hsc_dynamo_basophil_lineage.html), only 4 rows of switch gene-pairs were found. There is no GATA2-PLEK pair like the tutorial said. Do I miss something? Thanks! Best, YJ

dyn.vf.rank_jacobian_genes(adata, groups='cell_type', mode='switch')
adata.uns['switch']
Mon Mon_values  Meg Meg_values  MEP-like    MEP-like_values Ery Ery_values  Bas Bas_values  GMP-like    GMP-like_values HSC HSC_values  Neu Neu_values
0   GATA1 - SPI1    8.633986e-10    GATA1 - GATA1   1.777718e-08    GATA1 - SPI1    1.527943e-09    GATA1 - GATA1   1.985191e-09    GATA1 - GATA1   6.195085e-10    GATA1 - SPI1    2.298747e-09    GATA1 - SPI1    2.451547e-09    SPI1 - SPI1 2.616558e-11
1   SPI1 - GATA1    8.633986e-10    GATA1 - SPI1    1.086596e-10    SPI1 - GATA1    1.527943e-09    GATA1 - SPI1    4.970615e-10    GATA1 - SPI1    3.603877e-10    SPI1 - GATA1    2.298747e-09    SPI1 - GATA1    2.451547e-09    GATA1 - SPI1    8.766959e-12
2   GATA1 - GATA1   9.225353e-11    SPI1 - GATA1    1.086596e-10    GATA1 - GATA1   5.012768e-11    SPI1 - GATA1    4.970615e-10    SPI1 - GATA1    3.603877e-10    GATA1 - GATA1   0.000000e+00    SPI1 - SPI1 3.050378e-14    SPI1 - GATA1    8.766959e-12
3   SPI1 - SPI1 2.043805e-13    SPI1 - SPI1 2.603615e-11    SPI1 - SPI1 6.670578e-13    SPI1 - SPI1 4.424103e-12    SPI1 - SPI1 3.731447e-11    SPI1 - SPI1 0.000000e+00    GATA1 - GATA1   0.000000e+00    GATA1 - GATA1   7.800211e-13
Xiaojieqiu commented 2 years ago

looks like you only calculated Jacobian between SPI1 and GATA1. You should use all top pca genes, etc.

hyjforesight commented 2 years ago

Hello @Xiaojieqiu Thanks for the response. The tutorial did as below

adata = dyn.sample_data.hematopoiesis()
adata
AnnData object with n_obs × n_vars = 1947 × 1956
    obs: 'batch', 'time', 'cell_type', 'nGenes', 'nCounts', 'pMito', 'pass_basic_filter', 'new_Size_Factor', 'initial_new_cell_size', 'total_Size_Factor', 'initial_total_cell_size', 'spliced_Size_Factor', 'initial_spliced_cell_size', 'unspliced_Size_Factor', 'initial_unspliced_cell_size', 'Size_Factor', 'initial_cell_size', 'ntr', 'cell_cycle_phase', 'leiden', 'control_point_pca', 'inlier_prob_pca', 'obs_vf_angle_pca', 'pca_ddhodge_div', 'pca_ddhodge_potential', 'acceleration_pca', 'curvature_pca', 'n_counts', 'mt_frac', 'jacobian_det_pca', 'manual_selection', 'divergence_pca', 'curv_leiden', 'curv_louvain', 'SPI1->GATA1_jacobian', 'jacobian', 'umap_ori_leiden', 'umap_ori_louvain', 'umap_ddhodge_div', 'umap_ddhodge_potential', 'curl_umap', 'divergence_umap', 'acceleration_umap', 'control_point_umap_ori', 'inlier_prob_umap_ori', 'obs_vf_angle_umap_ori', 'curvature_umap_ori'
    var: 'gene_name', 'gene_id', 'nCells', 'nCounts', 'pass_basic_filter', 'use_for_pca', 'frac', 'ntr', 'time_3_alpha', 'time_3_beta', 'time_3_gamma', 'time_3_half_life', 'time_3_alpha_b', 'time_3_alpha_r2', 'time_3_gamma_b', 'time_3_gamma_r2', 'time_3_gamma_logLL', 'time_3_delta_b', 'time_3_delta_r2', 'time_3_bs', 'time_3_bf', 'time_3_uu0', 'time_3_ul0', 'time_3_su0', 'time_3_sl0', 'time_3_U0', 'time_3_S0', 'time_3_total0', 'time_3_beta_k', 'time_3_gamma_k', 'time_5_alpha', 'time_5_beta', 'time_5_gamma', 'time_5_half_life', 'time_5_alpha_b', 'time_5_alpha_r2', 'time_5_gamma_b', 'time_5_gamma_r2', 'time_5_gamma_logLL', 'time_5_bs', 'time_5_bf', 'time_5_uu0', 'time_5_ul0', 'time_5_su0', 'time_5_sl0', 'time_5_U0', 'time_5_S0', 'time_5_total0', 'time_5_beta_k', 'time_5_gamma_k', 'use_for_dynamics', 'gamma', 'gamma_r2', 'use_for_transition', 'gamma_k', 'gamma_b'
    uns: 'PCs', 'VecFld_pca', 'VecFld_umap', 'X_umap_neighbors', 'cell_phase_genes', 'cell_type_colors', 'dynamics', 'explained_variance_ratio_', 'feature_selection', 'grid_velocity_pca', 'grid_velocity_umap', 'grid_velocity_umap_ori_perturbation', 'grid_velocity_umap_test', 'jacobian_pca', 'leiden', 'neighbors', 'pca_mean', 'pp', 'response'
    obsm: 'X', 'X_pca', 'X_pca_SparseVFC', 'X_umap', 'X_umap_SparseVFC', 'X_umap_ori_perturbation', 'X_umap_test', 'acceleration_pca', 'acceleration_umap', 'cell_cycle_scores', 'curvature_pca', 'curvature_umap', 'j_delta_x_perturbation', 'velocity_pca', 'velocity_pca_SparseVFC', 'velocity_umap', 'velocity_umap_SparseVFC', 'velocity_umap_ori_perturbation', 'velocity_umap_test'
    layers: 'M_n', 'M_nn', 'M_t', 'M_tn', 'M_tt', 'X_new', 'X_total', 'velocity_alpha_minus_gamma_s'
    obsp: 'X_umap_connectivities', 'X_umap_distances', 'connectivities', 'cosine_transition_matrix', 'distances', 'fp_transition_rate', 'moments_con', 'pca_ddhodge', 'perturbation_transition_matrix', 'umap_ddhodge'

dyn.vf.rank_jacobian_genes(bif_adata, groups='cell_type', mode='switch')

But there is no bif_adata?? So I did

adata = dyn.sample_data.hematopoiesis()
dyn.vf.rank_jacobian_genes(adata, groups='cell_type', mode='switch')
adata.uns['switch']
Mon | Mon_values | Meg | Meg_values | MEP-like | MEP-like_values | Ery | Ery_values | Bas | Bas_values | GMP-like | GMP-like_values | HSC | HSC_values | Neu | Neu_values
-- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | --
GATA1 - SPI1 | 8.633986e-10 | GATA1 - GATA1 | 1.777718e-08 | GATA1 - SPI1 | 1.527943e-09 | GATA1 - GATA1 | 1.985191e-09 | GATA1 - GATA1 | 6.195085e-10 | GATA1 - SPI1 | 2.298747e-09 | GATA1 - SPI1 | 2.451547e-09 | SPI1 - SPI1 | 2.616558e-11
SPI1 - GATA1 | 8.633986e-10 | GATA1 - SPI1 | 1.086596e-10 | SPI1 - GATA1 | 1.527943e-09 | GATA1 - SPI1 | 4.970615e-10 | GATA1 - SPI1 | 3.603877e-10 | SPI1 - GATA1 | 2.298747e-09 | SPI1 - GATA1 | 2.451547e-09 | GATA1 - SPI1 | 8.766959e-12
GATA1 - GATA1 | 9.225353e-11 | SPI1 - GATA1 | 1.086596e-10 | GATA1 - GATA1 | 5.012768e-11 | SPI1 - GATA1 | 4.970615e-10 | SPI1 - GATA1 | 3.603877e-10 | GATA1 - GATA1 | 0.000000e+00 | SPI1 - SPI1 | 3.050378e-14 | SPI1 - GATA1 | 8.766959e-12
SPI1 - SPI1 | 2.043805e-13 | SPI1 - SPI1 | 2.603615e-11 | SPI1 - SPI1 | 6.670578e-13 | SPI1 - SPI1 | 4.424103e-12 | SPI1 - SPI1 | 3.731447e-11 | SPI1 - SPI1 | 0.000000e+00 | GATA1 - GATA1 | 0.000000e+00 | GATA1 - GATA1 | 7.800211e-13

I didn't limit the genes to GATA1 and SP1. Thanks! Best, YJ

Xiaojieqiu commented 2 years ago

you need to run jacobian with all genes

Xiaojieqiu commented 2 years ago

in more details, you need to do this first: dyn.vf.jacobian(adata_labeling, regulators=adata_labeling.var_names, effectors=adata_labeling.var_names)

hyjforesight commented 2 years ago

Thanks @Xiaojieqiu. If I didn't misunderstand, I think this code line is missing in the tutorial. That's why it made me confused. Thanks for the response!

hyjforesight commented 2 years ago

Hello @Xiaojieqiu Could you please tell me how you get the table showing Bas_Ery_curv? image I can only get the table like below. Note that the first row has nothing to do with Bas_Ery_curv. image

dyn.vf.jacobian(adata, regulators=adata.var_names, effectors=adata.var_names, sample_ncells=1000, basis='pca', Qkey='PCs', method='analytical')
dyn.vf.rank_jacobian_genes(adata, groups='cell_type', mode='switch')
adata.uns['switch'].to_csv('C:/Users/hyjfo/Documents/rank_jacobian_genes.csv')

Thanks! Best, YJ

github-actions[bot] commented 2 years ago

This issue is stale because it has been open 90 days with no activity. Remove stale label or comment or this will be closed in 14 days