Note

This page was generated from 3_morphogenesis_feature.ipynb. Interactive online version: . Some tutorial content may look better in light mode.

3.morphogenesis features#

Packages#

[1]:

import os
import sys

import gseapy as gp
import numpy as np
import dynamo as dyn
import spateo as st

/home/yao/anaconda3/envs/BGIpy38_tf2/lib/python3.8/site-packages/requests/__init__.py:102: RequestsDependencyWarning: urllib3 (1.26.9) or chardet (5.0.0)/charset_normalizer (2.0.12) doesn't match a supported version!
  warnings.warn("urllib3 ({}) or chardet ({})/charset_normalizer ({}) doesn't match a supported "

|-----> setting visualization default mode in dynamo. Your customized matplotlib settings might be overritten.

/home/yao/.local/lib/python3.8/site-packages/nxviz/__init__.py:18: UserWarning:
nxviz has a new API! Version 0.7.3 onwards, the old class-based API is being
deprecated in favour of a new API focused on advancing a grammar of network
graphics. If your plotting code depends on the old API, please consider
pinning nxviz at version 0.7.3, as the new API will break your old code.

To check out the new API, please head over to the docs at
https://ericmjl.github.io/nxviz/ to learn more. We hope you enjoy using it!

(This deprecation message will go away in version 1.0.)

  warnings.warn(
/home/yao/anaconda3/envs/BGIpy38_tf2/lib/python3.8/site-packages/spaghetti/network.py:36: FutureWarning:

The next major release of pysal/spaghetti (2.0.0) will drop support for all ``libpysal.cg`` geometries. This change is a first step in refactoring ``spaghetti`` that is expected to result in dramatically reduced runtimes for network instantiation and operations. Users currently requiring network and point pattern input as ``libpysal.cg`` geometries should prepare for this simply by converting to ``shapely`` geometries.

Data source#

[2]:

# The anndata object is reconstructed by ``2.morphogenesis vector fields and trajectories``
stage_adata = st.sample_data.drosophila(filename="E7-9h_cellbin_tdr_v3_midgut.h5ad")
stage_adata = stage_adata[:, np.sum(stage_adata.layers["counts_X"], axis=0) != 0]
stage_adata.X = stage_adata.layers["counts_X"].copy()
stage_adata.uns["pp"] = {}
dyn.pp.normalize_cell_expr_by_size_factors(adata=stage_adata, layers="X", X_total_layers=True,  skip_log=True)

# The pc model and mesh model are reconstructed by ``2.morphogenesis vector fields and trajectories``
stage_pc = st.tdr.read_model(f"vector_pc_model_midgut.vtk")
stage_mesh = st.tdr.read_model(f"vector_mesh_model_midgut.vtk")

cpo = [(531.4285063139628 / 1.5, 1120.6071331378873 / 1.5, 226.6879332984771 / 1.5),
 (1.9670869138005287, -6.902875264241757, -2.2120172004343885),
 (-0.1349909265808914, -0.13586573883518835, 0.9814876212931068)]

stage_adata

|-----> rounding expression data of layer: X during size factor calculation
|-----> size factor normalize following layers: ['X']
|-----? `total_szfactor` is not `None` and it is not in adata object.
|-----> skipping log transformation as input requires...
|-----> applying None to layer<X>
|-----> set adata <X> to normalized data.
|-----> <insert> pp.norm_method to uns in AnnData Object.

[2]:

AnnData object with n_obs × n_vars = 2326 × 7685
    obs: 'area', 'slices', 'nGenes', 'nCounts', 'pMito', 'nMito', 'pass_basic_filter', 'scc', 'auto_anno', 'anno_cell_type', 'anno_tissue', 'anno_germ_layer', 'actual_stage', 'Size_Factor', 'initial_cell_size'
    uns: 'PCs', 'VecFld_morpho', '__type', 'auto_anno_result', 'dendrogram_anno_cell_type', 'dendrogram_anno_germ_layer', 'dendrogram_anno_tissue', 'explained_variance_ratio_', 'fate_morpho', 'latter_models_align', 'neighbors', 'pca_mean', 'pca_valid_ind', 'pearson_X_neighbors', 'pp', 'rank_genes_groups', 'rank_genes_groups_anno_cell_type', 'rank_genes_groups_anno_germ_layer', 'rank_genes_groups_anno_tissue', 'scc', 'spatial'
    obsm: '3d_align_spatial', 'V_cells_mapping', 'X_cells_mapping', 'align_spatial', 'bbox', 'before_3d_align_spatial', 'contour', 'pearson_X_pca', 'pearson_X_umap', 'spatial', 'tdr_spatial'
    layers: 'counts_X', 'log1p_X', 'pearson_X', 'spliced', 'unspliced'

Generate the trajectory model#

[3]:

trajectory_model = st.tdr.construct_trajectory(
        adata=stage_adata,
        fate_key="fate_morpho",
        key_added="obs_index",
        n_sampling=500,
        sampling_method="trn",
        label=np.asarray(stage_adata.obs.index),
        )

|-----> [Running TRN] in progress: 100.0000%
|-----> [Running TRN] finished [23.7339s]

Morphogenesis features#

[5]:

glm_dict = {}

Velocity#

[6]:

key = "velocity"
st.tdr.morphofield_velocity(adata=stage_adata, vf_key="VecFld_morpho", key_added=key,)
stage_adata.obsm["velocity"]

[6]:

array([[-1.91148383e-03, -1.82671905e-03, -2.71598635e-03],
       [-1.95678207e-03, -1.91036185e-03, -2.87477111e-03],
       [-1.55106979e-03, -1.68945969e-03, -2.41216961e-03],
       ...,
       [-2.10651796e-06,  1.09092371e-04,  1.71498103e-03],
       [ 1.89100976e-06,  8.49589979e-05,  1.72119623e-03],
       [-8.09060396e-05,  1.64108170e-04,  1.67981541e-03]])

Acceleration#

[7]:

key = "acceleration"
glm_key = f"glm_degs_{key}"
st.tdr.morphofield_acceleration(adata=stage_adata, vf_key="VecFld_morpho", key_added=key,)
st.tl.glm_degs(adata=stage_adata, layer=None, fullModelFormulaStr=f'~cr({key}, df=3)', key_added=glm_key, qval_threshold=0.01, llf_threshold=-2500)
glm_data = stage_adata.uns[glm_key]["glm_result"]
glm_dict[key] = glm_data
glm_data

|-----> [Calculating acceleration] in progress: 100.0000%
|-----> [Calculating acceleration] finished [0.0412s]
|-----? Gene expression matrix must be normalized by the size factor, please check if the input gene expression matrix is correct.If you don't have the size factor normalized gene expression matrix, please run `dynamo.pp.normalize_cell_expr_by_size_factors(skip_log = True)`.
|-----> [Detecting genes via Generalized Additive Models (GAMs)] in progress: 100.0000%
|-----> [Detecting genes via Generalized Additive Models (GAMs)] finished [375.0744s]

[7]:

	status	family	log-likelihood	pval	qval
Bacc	ok	NB2	-3671.859375	4.873121e-11	7.489986e-09
28SrRNA-Psi:CR45862	ok	NB2	-4358.506348	6.277172e-10	6.432008e-08
Cys	ok	NB2	-4643.826172	4.387060e-08	2.148794e-06
28SrRNA-Psi:CR45859	ok	NB2	-4632.636230	2.513353e-06	4.769165e-05
eEF1beta	ok	NB2	-3303.432861	1.365682e-05	1.904767e-04
Nph	ok	NB2	-3012.880371	3.456467e-05	3.855290e-04
fax	ok	NB2	-2740.145508	5.348858e-05	5.373330e-04
His3.3A	ok	NB2	-3064.583496	6.052202e-05	5.887490e-04
Fib	ok	NB2	-2566.735840	2.022182e-04	1.532591e-03
awd	ok	NB2	-3193.938477	6.848338e-04	3.859636e-03
mt:lrRNA	ok	NB2	-6083.094238	7.883165e-04	4.318041e-03
tsr	ok	NB2	-3093.676758	9.017585e-04	4.782618e-03
28SrRNA-Psi:CR40596	ok	NB2	-3915.597412	1.051448e-03	5.379748e-03
pAbp	ok	NB2	-3751.640137	1.692311e-03	7.820453e-03

[8]:

st.pl.glm_fit(
    adata=stage_adata,
    gene=glm_data.index.tolist(),
    ncols=5,
    feature_x=key,
    feature_y="expression",
    feature_fit="mu",
    glm_key=glm_key,
    lowess=False,
    frac=0.2,
    save_show_or_return="all",
)
st.pl.glm_fit(
    adata=stage_adata,
    gene=glm_data.index.tolist(),
    ncols=5,
    feature_x=key,
    feature_y="expression",
    feature_fit="mu",
    glm_key=glm_key,
    lowess=True,
    frac=0.2,
    save_show_or_return="all",
)
st.pl.acceleration(
    adata=stage_adata,
    model=st.tdr.collect_models([stage_pc, trajectory_model]),
    acceleration_key=key,
    colormap="default_cmap",
    jupyter="static",
    model_style=["points", "wireframe"],
    model_size=[5, 2],
    cpo=cpo,
    window_size=(1024, 1024),
    legend_kwargs=dict(title="", fmt="%.2e", legend_loc=(0.87, 0.3), label_font_size=20),
)

Saving figure to develop/develop_midgut/genesis_feature/E7-9h_midgut_acceleration/E7-9h_midgut_glm_degs_fit_acceleration.pdf...
Done

../../../_images/tutorials_notebooks_6_morphogenesis_3_morphogenesis_feature_13_1.png

Saving figure to develop/develop_midgut/genesis_feature/E7-9h_midgut_acceleration/E7-9h_midgut_glm_degs_fit_acceleration_lowess.pdf...
Done

../../../_images/tutorials_notebooks_6_morphogenesis_3_morphogenesis_feature_13_3.png

../../../_images/tutorials_notebooks_6_morphogenesis_3_morphogenesis_feature_13_4.png

[9]:

st.tdr.add_model_labels(
    model=trajectory_model,
    key_added=key,
    labels=np.asarray(stage_adata[np.asarray(trajectory_model.point_data["obs_index"])].obs[key]),
    where="point_data",
    inplace=True,
)
cells_linear_models = st.tdr.construct_genesis_X(
    stages_X=[stage_adata.obsm["3d_align_spatial"], stage_adata.obsm["X_cells_mapping"]],
    n_spacing=100-1,
    key_added=key,
    label=[np.asarray(stage_adata.obs[key])] * 100,
)
st.pl.three_d_animate(
     models=cells_linear_models,
    stable_model=trajectory_model,
    key=key,
    stable_kwargs=dict(
        key=key,
        model_style="wireframe",
        model_size=2,
        ambient=0.5,
        opacity=0.3,
        colormap="default_cmap",
        show_legend=False,
    ),
    filename=f"E7-9h_midgut_{key}_linear_movie_with_trajectory.mp4",
    model_style="points",
    model_size=5,
    ambient=0.5,
    colormap="default_cmap",
    show_legend=True,
    legend_kwargs=dict(title="", fmt="%.2e", legend_loc=(0.87, 0.3), label_font_size=20),
    off_screen=True,
    jupyter=False,
    background="black",
    cpo=cpo,
    window_size=(1024, 1024),
)
st.pl.three_d_animate(
    models=cells_linear_models,
    stable_model=None,
    key=key,
    filename=f"E7-9h_midgut_{key}_linear_movie.mp4",
    model_style="points",
    model_size=5,
    ambient=0.5,
    colormap="default_cmap",
    show_legend=True,
    legend_kwargs=dict(title="", fmt="%.2e", legend_loc=(0.87, 0.3), label_font_size=20),
    off_screen=True,
    jupyter=False,
    background="black",
    cpo=cpo,
    window_size=(1024, 1024),
)

Curvature#

[10]:

key = "curvature"
glm_key = f"glm_degs_{key}"

st.tdr.morphofield_curvature(adata=stage_adata, vf_key="VecFld_morpho", key_added=key,)
st.tl.glm_degs(adata=stage_adata, layer=None, fullModelFormulaStr=f'~cr({key}, df=3)', key_added=glm_key, qval_threshold=0.01, llf_threshold=-2500)
glm_data = stage_adata.uns[glm_key]["glm_result"]
glm_dict[key] = glm_data
glm_data

|-----> [Calculating acceleration] in progress: 100.0000%
|-----> [Calculating acceleration] finished [0.0304s]
|-----> [Calculating curvature] in progress: 100.0000%
|-----> [Calculating curvature] finished [0.0381s]
|-----? Gene expression matrix must be normalized by the size factor, please check if the input gene expression matrix is correct.If you don't have the size factor normalized gene expression matrix, please run `dynamo.pp.normalize_cell_expr_by_size_factors(skip_log = True)`.
|-----> [Detecting genes via Generalized Additive Models (GAMs)] in progress: 100.0000%
|-----> [Detecting genes via Generalized Additive Models (GAMs)] finished [287.0102s]

[10]:

	status	family	log-likelihood	pval	qval
Ance	ok	NB2	-3339.842529	5.916981e-13	1.136800e-09
Fib	ok	NB2	-2560.309814	3.273910e-07	8.064102e-06
28SrRNA-Psi:CR40596	ok	NB2	-3911.105713	1.177748e-05	1.544538e-04
mt:ND5	ok	NB2	-2526.593506	2.593650e-05	2.823258e-04
RpL24-like	ok	NB2	-3122.291748	2.842352e-05	3.025412e-04
Cys	ok	NB2	-4650.455078	3.319294e-05	3.410264e-04
28SrRNA-Psi:CR45862	ok	NB2	-4370.237793	7.809990e-05	6.751381e-04
RNASEK	ok	NB2	-2597.170898	1.023775e-04	8.238442e-04
His2Av	ok	NB2	-3063.656494	1.208075e-04	9.377835e-04
sesB	ok	NB2	-2753.985596	1.806184e-04	1.286111e-03
COX7C	ok	NB2	-2909.798096	3.349247e-04	2.114103e-03
Obp99a	ok	NB2	-3229.527588	4.188366e-04	2.513024e-03
Cyp1	ok	NB2	-3645.378418	8.155234e-04	4.231801e-03
chic	ok	NB2	-2536.232910	8.432105e-04	4.340303e-03
RpS16	ok	NB2	-4785.575684	1.652613e-03	7.349729e-03

[11]:

st.pl.glm_fit(
    adata=stage_adata,
    gene=glm_data.index.tolist(),
    ncols=5,
    feature_x=key,
    feature_y="expression",
    feature_fit="mu",
    glm_key=glm_key,
    lowess=False,
    frac=0.2,
    save_show_or_return="both",
)
st.pl.curvature(
    adata=stage_adata,
    model=st.tdr.collect_models([stage_pc, trajectory_model]),
    curvature_key=key,
    colormap="default_cmap",
    jupyter="static",
    model_style=["points", "wireframe"],
    model_size=[5, 2],
    cpo=cpo,
    window_size=(1024, 1024),
    legend_kwargs=dict(title="", fmt="%.2e", legend_loc=(0.87, 0.3), label_font_size=20)
)

Saving figure to develop/develop_midgut/genesis_feature/E7-9h_midgut_curvature/E7-9h_midgut_glm_degs_fit_curvature.pdf...
Done

../../../_images/tutorials_notebooks_6_morphogenesis_3_morphogenesis_feature_17_1.png

../../../_images/tutorials_notebooks_6_morphogenesis_3_morphogenesis_feature_17_2.png

[12]:

st.tdr.add_model_labels(
    model=trajectory_model,
    key_added=key,
    labels=np.asarray(stage_adata[np.asarray(trajectory_model.point_data["obs_index"])].obs[key]),
    where="point_data",
    inplace=True,
)
cells_linear_models = st.tdr.construct_genesis_X(
    stages_X=[stage_adata.obsm["3d_align_spatial"], stage_adata.obsm["X_cells_mapping"]],
    n_spacing=100-1,
    key_added=key,
    label=[np.asarray(stage_adata.obs[key])] * 100,
)
st.pl.three_d_animate(
     models=cells_linear_models,
    stable_model=trajectory_model,
    key=key,
    stable_kwargs=dict(
        key=key,
        model_style="wireframe",
        model_size=2,
        ambient=0.5,
        opacity=0.3,
        colormap="default_cmap",
        show_legend=False,
    ),
    filename=f"E7-9h_midgut_{key}_linear_movie_with_trajectory.mp4",
    model_style="points",
    model_size=5,
    ambient=0.5,
    colormap="default_cmap",
    show_legend=True,
    legend_kwargs=dict(title="", fmt="%.2e", legend_loc=(0.87, 0.3), label_font_size=20),
    off_screen=True,
    jupyter=False,
    background="black",
    cpo=cpo,
    window_size=(1024, 1024),
)
st.pl.three_d_animate(
    models=cells_linear_models,
    stable_model=None,
    key=key,
    filename=f"E7-9h_midgut_{key}_linear_movie.mp4",
    model_style="points",
    model_size=5,
    ambient=0.5,
    colormap="default_cmap",
    show_legend=True,
    legend_kwargs=dict(title="", fmt="%.2e", legend_loc=(0.87, 0.3), label_font_size=20),
    off_screen=True,
    jupyter=False,
    background="black",
    cpo=cpo,
    window_size=(1024, 1024),
)

Curl#

[13]:

key = "curl"
glm_key = f"glm_degs_{key}"

st.tdr.morphofield_curl(adata=stage_adata, vf_key="VecFld_morpho", key_added=key,)
st.tl.glm_degs(adata=stage_adata, layer=None, fullModelFormulaStr=f'~cr({key}, df=3)', key_added=glm_key, qval_threshold=0.01, llf_threshold=-2500)
glm_data = stage_adata.uns[glm_key]["glm_result"]
glm_dict[key] = glm_data
glm_data

Calculating 3-D curl: 100%|██████████| 2326/2326 [00:00<00:00, 57530.08it/s]

|-----? Gene expression matrix must be normalized by the size factor, please check if the input gene expression matrix is correct.If you don't have the size factor normalized gene expression matrix, please run `dynamo.pp.normalize_cell_expr_by_size_factors(skip_log = True)`.

|-----> [Detecting genes via Generalized Additive Models (GAMs)] in progress: 100.0000%
|-----> [Detecting genes via Generalized Additive Models (GAMs)] finished [267.5932s]

[13]:

	status	family	log-likelihood	pval	qval
Fib	ok	NB2	-2560.906738	5.947219e-07	0.000015
pAbp	ok	NB2	-3744.795410	1.802494e-06	0.000038
COX8	ok	NB2	-3233.309570	2.275150e-05	0.000267
28SrRNA-Psi:CR40596	ok	NB2	-3912.000488	2.882065e-05	0.000321
Ance	ok	NB2	-3358.986572	1.219921e-04	0.001011
tsr	ok	NB2	-3092.293213	2.260249e-04	0.001625
Nph	ok	NB2	-3015.566162	5.069769e-04	0.003041
RpL24-like	ok	NB2	-3125.715576	8.723374e-04	0.004005
Bacc	ok	NB2	-3688.385254	7.326030e-04	0.004059
28SrRNA-Psi:CR45859	ok	NB2	-4638.458984	8.491511e-04	0.004522
fax	ok	NB2	-2743.038330	9.650021e-04	0.005024
NHP2	ok	NB2	-2605.645508	1.179141e-03	0.005976
CNBP	ok	NB2	-3570.270264	1.538176e-03	0.006671
ATPsynG	ok	NB2	-3090.621582	1.784211e-03	0.008645
mt:ND1	ok	NB2	-3178.442871	2.424773e-03	0.009808

[14]:

st.pl.glm_fit(
    adata=stage_adata,
    gene=glm_data.index.tolist(),
    ncols=5,
    feature_x=key,
    feature_y="expression",
    feature_fit="mu",
    glm_key=glm_key,
    lowess=False,
    frac=0.2,
    save_show_or_return="both",
)
st.pl.curl(
    adata=stage_adata,
    model=st.tdr.collect_models([stage_pc, trajectory_model]),
    curl_key=key,
    colormap="default_cmap",
    jupyter="static",
    model_style=["points", "wireframe"],
    model_size=[5, 2],
    cpo=cpo,
    window_size=(1024, 1024),
    legend_kwargs=dict(title="", fmt="%.2e", legend_loc=(0.87, 0.3), label_font_size=20),
)

Saving figure to develop/develop_midgut/genesis_feature/E7-9h_midgut_curl/E7-9h_midgut_glm_degs_fit_curl.pdf...
Done

../../../_images/tutorials_notebooks_6_morphogenesis_3_morphogenesis_feature_21_1.png

../../../_images/tutorials_notebooks_6_morphogenesis_3_morphogenesis_feature_21_2.png

[15]:

st.tdr.add_model_labels(
    model=trajectory_model,
    key_added=key,
    labels=np.asarray(stage_adata[np.asarray(trajectory_model.point_data["obs_index"])].obs[key]),
    where="point_data",
    inplace=True,
)
cells_linear_models = st.tdr.construct_genesis_X(
    stages_X=[stage_adata.obsm["3d_align_spatial"], stage_adata.obsm["X_cells_mapping"]],
    n_spacing=100-1,
    key_added=key,
    label=[np.asarray(stage_adata.obs[key])] * 100,
)
st.pl.three_d_animate(
     models=cells_linear_models,
    stable_model=trajectory_model,
    key=key,
    stable_kwargs=dict(
        key=key,
        model_style="wireframe",
        model_size=2,
        ambient=0.5,
        opacity=0.3,
        colormap="default_cmap",
        show_legend=False,
    ),
    filename=f"E7-9h_midgut_{key}_linear_movie_with_trajectory.mp4",
    model_style="points",
    model_size=5,
    ambient=0.5,
    colormap="default_cmap",
    show_legend=True,
    legend_kwargs=dict(title="", fmt="%.2e", legend_loc=(0.87, 0.3), label_font_size=20),
    off_screen=True,
    jupyter=False,
    background="black",
    cpo=cpo,
    window_size=(1024, 1024),
)
st.pl.three_d_animate(
    models=cells_linear_models,
    stable_model=None,
    key=key,
    filename=f"E7-9h_midgut_{key}_linear_movie.mp4",
    model_style="points",
    model_size=5,
    ambient=0.5,
    colormap="default_cmap",
    show_legend=True,
    legend_kwargs=dict(title="", fmt="%.2e", legend_loc=(0.87, 0.3), label_font_size=20),
    off_screen=True,
    jupyter=False,
    background="black",
    cpo=cpo,
    window_size=(1024, 1024),
)

Torsion#

[16]:

key = "torsion"
glm_key = f"glm_degs_{key}"

st.tdr.morphofield_torsion(adata=stage_adata, vf_key="VecFld_morpho", key_added=key,)
st.tl.glm_degs(adata=stage_adata, layer=None, fullModelFormulaStr=f'~cr({key}, df=3)', key_added=glm_key, qval_threshold=0.01, llf_threshold=-2500)
glm_data = stage_adata.uns[glm_key]["glm_result"]
glm_dict[key] = glm_data
glm_data

|-----> [Calculating acceleration] in progress: 100.0000%
|-----> [Calculating acceleration] finished [0.0342s]

Calculating torsion: 100%|██████████| 2326/2326 [00:00<00:00, 123461.80it/s]

|-----? Gene expression matrix must be normalized by the size factor, please check if the input gene expression matrix is correct.If you don't have the size factor normalized gene expression matrix, please run `dynamo.pp.normalize_cell_expr_by_size_factors(skip_log = True)`.

|-----> [Detecting genes via Generalized Additive Models (GAMs)] in progress: 100.0000%
|-----> [Detecting genes via Generalized Additive Models (GAMs)] finished [332.7657s]

[16]:

	status	family	log-likelihood	pval	qval
28SrRNA-Psi:CR45862	ok	NB2	-4346.507812	3.861141e-15	2.472739e-12
awd	ok	NB2	-3186.442383	3.803082e-07	1.154048e-05
COX5B	ok	NB2	-2744.803711	2.241822e-06	4.618874e-05
COX8	ok	NB2	-3231.869873	5.392043e-06	9.460696e-05
28SrRNA-Psi:CR45859	ok	NB2	-4634.142578	1.133787e-05	1.689693e-04
bnb	ok	NB2	-2515.703857	1.134524e-05	1.689693e-04
Cyp1	ok	NB2	-3642.751221	5.894959e-05	6.097518e-04
COX5A	ok	NB2	-2502.031494	6.766557e-05	6.806412e-04
COX7C	ok	NB2	-2909.620850	2.804796e-04	2.050890e-03
mt:lrRNA	ok	NB2	-6082.588867	4.755917e-04	3.090644e-03
Df31	ok	NB2	-3453.745605	5.343023e-04	3.363636e-03
Ance	ok	NB2	-3360.501221	5.547805e-04	3.470894e-03
eEF1beta	ok	NB2	-3307.284668	6.429467e-04	3.881419e-03
CG11267	ok	NB2	-2541.744141	8.072325e-04	4.622639e-03

[17]:

st.pl.glm_fit(
    adata=stage_adata,
    gene=glm_data.index.tolist(),
    ncols=5,
    feature_x=key,
    feature_y="expression",
    feature_fit="mu",
    glm_key=glm_key,
    lowess=False,
    frac=0.2,
    save_show_or_return="both",
)
st.pl.torsion(
    adata=stage_adata,
    model=st.tdr.collect_models([stage_pc, trajectory_model]),
    torsion_key=key,
    colormap="default_cmap",
    jupyter="static",
    model_style=["points", "wireframe"],
    model_size=[5, 2],
    cpo=cpo,
    window_size=(1024, 1024),
    legend_kwargs=dict(title="", fmt="%.2e", legend_loc=(0.87, 0.3), label_font_size=20),
)

Saving figure to develop/develop_midgut/genesis_feature/E7-9h_midgut_torsion/E7-9h_midgut_glm_degs_fit_torsion.pdf...
Done

../../../_images/tutorials_notebooks_6_morphogenesis_3_morphogenesis_feature_25_1.png

../../../_images/tutorials_notebooks_6_morphogenesis_3_morphogenesis_feature_25_2.png

[18]:

st.tdr.add_model_labels(
    model=trajectory_model,
    key_added=key,
    labels=np.asarray(stage_adata[np.asarray(trajectory_model.point_data["obs_index"])].obs[key]),
    where="point_data",
    inplace=True,
)
cells_linear_models = st.tdr.construct_genesis_X(
    stages_X=[stage_adata.obsm["3d_align_spatial"], stage_adata.obsm["X_cells_mapping"]],
    n_spacing=100-1,
    key_added=key,
    label=[np.asarray(stage_adata.obs[key])] * 100,
)
st.pl.three_d_animate(
     models=cells_linear_models,
    stable_model=trajectory_model,
    key=key,
    stable_kwargs=dict(
        key=key,
        model_style="wireframe",
        model_size=2,
        ambient=0.5,
        opacity=0.3,
        colormap="default_cmap",
        show_legend=False,
    ),
    filename=f"E7-9h_midgut_{key}_linear_movie_with_trajectory.mp4",
    model_style="points",
    model_size=5,
    ambient=0.5,
    colormap="default_cmap",
    show_legend=True,
    legend_kwargs=dict(title="", fmt="%.2e", legend_loc=(0.87, 0.3), label_font_size=20),
    off_screen=True,
    jupyter=False,
    background="black",
    cpo=cpo,
    window_size=(1024, 1024),
)
st.pl.three_d_animate(
    models=cells_linear_models,
    stable_model=None,
    key=key,
    filename=f"E7-9h_midgut_{key}_linear_movie.mp4",
    model_style="points",
    model_size=5,
    ambient=0.5,
    colormap="default_cmap",
    show_legend=True,
    legend_kwargs=dict(title="", fmt="%.2e", legend_loc=(0.87, 0.3), label_font_size=20),
    off_screen=True,
    jupyter=False,
    background="black",
    cpo=cpo,
    window_size=(1024, 1024),
)

Jacobian#

[19]:

key = "jacobian"
st.tdr.morphofield_jacobian(adata=stage_adata, vf_key="VecFld_morpho", key_added=key,)
stage_adata.uns["jacobian"]

[19]:

array([[[ 2.08105242e-05,  2.70779418e-05,  3.14283350e-05, ...,
         -8.25586263e-06, -8.49332510e-06, -5.65412876e-06],
        [ 7.23267922e-05,  6.98323637e-05,  6.34322214e-05, ...,
          3.69536432e-07,  2.37314756e-07, -9.78666713e-07],
        [ 4.65323612e-05,  4.93631551e-05,  3.55578438e-05, ...,
         -3.28196112e-06, -3.21871336e-06, -5.28409549e-06]],

       [[ 1.67115495e-06,  7.37013386e-06,  1.69219829e-05, ...,
          7.19850134e-06,  7.21621752e-06,  2.70224705e-06],
        [ 5.42722509e-05,  5.09177864e-05,  5.47200113e-05, ...,
         -7.52587956e-06, -8.50232404e-06, -6.79560366e-06],
        [ 4.87511113e-05,  5.29128280e-05,  4.34329063e-05, ...,
          9.67988261e-07,  3.44231301e-08,  3.61775236e-06]],

       [[ 7.83754913e-06,  1.73498770e-05,  3.04055759e-05, ...,
          7.92778250e-06,  6.45248208e-06, -1.02374499e-05],
        [ 1.02038938e-04,  1.00946200e-04,  9.68161684e-05, ...,
          1.84122214e-05,  1.03350391e-05,  1.36330631e-05],
        [ 6.19007123e-05,  6.78370197e-05,  5.12851558e-05, ...,
          4.45638671e-05,  4.71304287e-05,  5.08665136e-05]]])

[20]:

st.pl.jacobian(
    adata=stage_adata,
    model=st.tdr.collect_models([stage_pc, trajectory_model]),
    jacobian_key="jacobian",
    colormap="default_cmap",
    jupyter="static",
    model_style=["points", "wireframe"],
    model_size=[2.5, 1],
    cpo=[cpo],
    window_size=(512*3, 512*3),
)

../../../_images/tutorials_notebooks_6_morphogenesis_3_morphogenesis_feature_29_0.png

GO analysis#

[21]:

gp.get_library_name(organism="fly")

[21]:

['Allele_LoF_Phenotypes_from_FlyBase_2017',
 'Allele_Phenotypes_from_FlyBase_2017',
 'Anatomy_AutoRIF',
 'Anatomy_AutoRIF_Predicted_zscore',
 'Anatomy_GeneRIF',
 'Anatomy_GeneRIF_Predicted_zscore',
 'Coexpression_Predicted_GO_Biological_Process_2018',
 'Coexpression_Predicted_GO_Cellular_Component_2018',
 'Coexpression_Predicted_GO_Molecular_Function_2018',
 'GO_Biological_Process_2018',
 'GO_Biological_Process_AutoRIF',
 'GO_Biological_Process_AutoRIF_Predicted_zscore',
 'GO_Biological_Process_GeneRIF',
 'GO_Biological_Process_GeneRIF_Predicted_zscore',
 'GO_Cellular_Component_2018',
 'GO_Cellular_Component_AutoRIF',
 'GO_Cellular_Component_AutoRIF_Predicted_zscore',
 'GO_Cellular_Component_GeneRIF',
 'GO_Cellular_Component_GeneRIF_Predicted_zscore',
 'GO_Molecular_Function_2018',
 'GO_Molecular_Function_AutoRIF',
 'GO_Molecular_Function_AutoRIF_Predicted_zscore',
 'GO_Molecular_Function_GeneRIF',
 'GO_Molecular_Function_GeneRIF_Predicted_zscore',
 'Human_Disease_from_FlyBase_2017',
 'InterPro_Domains_2019',
 'KEGG_2019',
 'PPI_Network_Hubs_from_DroID_2017',
 'Pfam_Domains_2019',
 'Phenotype_AutoRIF',
 'Phenotype_AutoRIF_Predicted_zscore',
 'Phenotype_GeneRIF',
 'Phenotype_GeneRIF_Predicted_zscore',
 'Putative_Regulatory_miRNAs_from_DroID_2015',
 'RNAi_Screens_from_GenomeRNAi_2017',
 'TF2DNA_2018',
 'Transcription_Factors_from_DroID_2015',
 'WikiPathways_2018']

[22]:

for key, data in glm_dict.items():
    if len(data.index.tolist()) > 0:
        save_folder = f"E7-9h_midgut_{key}"

        # GO-biologial process
        go_bp = gp.enrichr(
            gene_list=data.index.tolist(),
            gene_sets="GO_Biological_Process_2018",
            organism="fly",
            outdir=save_folder,
            no_plot=True,
            verbose=True
        )
        # GO-cellular component
        gp_cc = gp.enrichr(
            gene_list=data.index.tolist(),
            gene_sets="GO_Cellular_Component_2018",
            organism="fly",
            outdir=save_folder,
            no_plot=True,
            verbose=True
        )
        # GO-molecular function
        go_mf = gp.enrichr(
            gene_list=data.index.tolist(),
            gene_sets="GO_Molecular_Function_2018",
            organism="fly",
            outdir=save_folder,
            no_plot=True,
            verbose=True
        )

2022-10-16 14:11:43,473 Connecting to Enrichr Server to get latest library names
2022-10-16 14:11:44,859 Analysis name: , Enrichr Library: GO_Biological_Process_2018
2022-10-16 14:11:51,166 Save file of enrichment results: Job Id:231b5ebbe467f630d7a0ea9e7a9a39bd
2022-10-16 14:11:51,170 Done.

2022-10-16 14:11:51,171 Connecting to Enrichr Server to get latest library names
2022-10-16 14:11:52,551 Analysis name: , Enrichr Library: GO_Cellular_Component_2018
2022-10-16 14:11:56,542 Save file of enrichment results: Job Id:e5d4d3bb96e241ee70e016b34a31fc08
2022-10-16 14:11:56,546 Done.

2022-10-16 14:11:56,547 Connecting to Enrichr Server to get latest library names
2022-10-16 14:11:57,954 Analysis name: , Enrichr Library: GO_Molecular_Function_2018
2022-10-16 14:12:01,850 Save file of enrichment results: Job Id:44c1a53e9556e824aef5ea46234a3cd6
2022-10-16 14:12:01,852 Done.

2022-10-16 14:12:01,853 Connecting to Enrichr Server to get latest library names
2022-10-16 14:12:03,349 Analysis name: , Enrichr Library: GO_Biological_Process_2018
2022-10-16 14:12:07,927 Save file of enrichment results: Job Id:ed7c2d3f975ee261ffe541326810c97d
2022-10-16 14:12:07,932 Done.

2022-10-16 14:12:07,933 Connecting to Enrichr Server to get latest library names
2022-10-16 14:12:09,320 Analysis name: , Enrichr Library: GO_Cellular_Component_2018
2022-10-16 14:12:13,247 Save file of enrichment results: Job Id:027814f8da96a2305105a3fcc66febcb
2022-10-16 14:12:13,248 Done.

2022-10-16 14:12:13,249 Connecting to Enrichr Server to get latest library names
2022-10-16 14:12:14,632 Analysis name: , Enrichr Library: GO_Molecular_Function_2018
2022-10-16 14:12:18,544 Save file of enrichment results: Job Id:5fc10f71d4ca7b9f89ae182e57f3612f
2022-10-16 14:12:18,551 Done.

2022-10-16 14:12:18,553 Connecting to Enrichr Server to get latest library names
2022-10-16 14:12:19,954 Analysis name: , Enrichr Library: GO_Biological_Process_2018
2022-10-16 14:12:24,336 Save file of enrichment results: Job Id:f81cd3475f2e62680675ad00366355af
2022-10-16 14:12:24,342 Done.

2022-10-16 14:12:24,348 Connecting to Enrichr Server to get latest library names
2022-10-16 14:12:25,726 Analysis name: , Enrichr Library: GO_Cellular_Component_2018
2022-10-16 14:12:30,748 Save file of enrichment results: Job Id:89df4375762b1f95d12d56ee740c75a6
2022-10-16 14:12:30,752 Done.

2022-10-16 14:12:30,754 Connecting to Enrichr Server to get latest library names
2022-10-16 14:12:32,168 Analysis name: , Enrichr Library: GO_Molecular_Function_2018
2022-10-16 14:12:36,127 Save file of enrichment results: Job Id:2600b4c6db1bec44d1f67347d2392e99
2022-10-16 14:12:36,133 Done.

2022-10-16 14:12:36,135 Connecting to Enrichr Server to get latest library names
2022-10-16 14:12:37,611 Analysis name: , Enrichr Library: GO_Biological_Process_2018
2022-10-16 14:12:41,846 Save file of enrichment results: Job Id:c02d8a0119d766863cf00e35ca9e3794
2022-10-16 14:12:41,852 Done.

2022-10-16 14:12:41,854 Connecting to Enrichr Server to get latest library names
2022-10-16 14:12:43,244 Analysis name: , Enrichr Library: GO_Cellular_Component_2018
2022-10-16 14:12:47,747 Save file of enrichment results: Job Id:ee98c8b120087c155ce5116a5dfff9c7
2022-10-16 14:12:47,753 Done.

2022-10-16 14:12:47,757 Connecting to Enrichr Server to get latest library names
2022-10-16 14:12:49,105 Analysis name: , Enrichr Library: GO_Molecular_Function_2018
2022-10-16 14:12:53,686 Save file of enrichment results: Job Id:c4898910b80317f3cbee0ef3c690f369
2022-10-16 14:12:53,690 Done.

Genes expression pattern#

[23]:

gene_list = list(set([gene for data in glm_dict.values() for gene in data.index.tolist()]))

[24]:

for gene in gene_list:
    obs_index = stage_pc.point_data["obs_index"].tolist()
    mkg_exp = stage_adata[obs_index, gene].X.flatten().astype(np.float64)
    st.tdr.add_model_labels(model=stage_pc, labels=mkg_exp, key_added=gene, where="point_data", inplace=True)

    mkg_amap = mkg_exp.copy()
    mkg_amap = mkg_amap / mkg_amap.max()
    mkg_amap[mkg_amap <= 0.5] = 0.5

    rename_gene = str(gene).replace(":", "-") if ":" in gene else gene
    st.pl.three_d_plot(
        model=stage_pc,
        key=gene,
        model_style="points",
        model_size=5,
        # opacity=mkg_amap,
        ambient=0.5,
        colormap="default_cmap",
        text=f"\nGene: {gene}",
        show_legend=True,
        legend_kwargs=dict(title=""),
        off_screen=True,
        jupyter=False,
        background="black",
        cpo=cpo,
        window_size=(1024, 1024),
        filename=f"pc_model_gene_{rename_gene}.pdf"
    )

[ ]: