Neurogenesis
[1]:
import getopt
import pandas as pd
import numpy as np
import pickle
import os
import sys
import phlower
import scanpy as sc
from datetime import datetime
from matplotlib import pyplot as plt
from collections import Counter
[2]:
## from scvelo.datasets.dentategyrus_lamanno()
print(datetime.now(),"loading...")
adata = phlower.dataset.neurogenesis()
adata.var_names_make_unique()
print(datetime.now(),"mt...")
# annotate the group of mitochondrial genes as "mt"
adata.var["mt"] = adata.var_names.str.startswith("mt-")
sc.pp.calculate_qc_metrics(
adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
)
2025-01-14 15:52:33.429898 loading...
2025-01-14 15:52:34.920234 mt...
[3]:
print(datetime.now(),"preprocessing...")
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
adata = adata[:, adata.var.highly_variable]
sc.pp.regress_out(adata, ["total_counts", "pct_counts_mt"])
print(datetime.now(),"pca...")
sc.tl.pca(adata, svd_solver="arpack")
adata.obs['group'] = adata.obs['clusters']
2025-01-14 15:52:35.494529 preprocessing...
2025-01-14 15:53:05.777980 pca...
[4]:
phlower.ext.ddhodge(adata, basis='X_pca', roots=(adata.obs.group=="nIPC"), k=12, npc=40,ndc=5,s=1, verbose=True)
2025-01-14 15:53:08.446259 distance_matrix
2025-01-14 15:53:41.729824 Diffusionmaps:
done.
2025-01-14 15:55:10.265315 diffusion distance:
2025-01-14 15:55:21.895666 transition matrix:
2025-01-14 15:55:34.563674 graph from A
2025-01-14 15:55:35.451916 Rewiring:
2025-01-14 15:55:35.452000 div(g_o)...
2025-01-14 15:57:24.801157 edge weight...
2025-01-14 15:57:28.130695 cholesky solve ax=b...
+ 1e-06 I is positive-definite
2025-01-14 15:58:57.280354 grad...
2025-01-14 16:17:55.681036 potential.
2025-01-14 16:36:47.933767 ddhodge done.
done.
2025-01-14 16:36:49.143244 calculate layouts
2025-01-14 17:33:56.925509 done
[5]:
phlower.tl.construct_delaunay_persistence(adata, start_n=10, end_n=10, circle_quant=0.7, min_persistence=0.1)
start clusters ['CA', 'Granule', 'CA2-3-4', 'RadialGlia', 'ImmAstro', 'ImmGranule2', 'OPC', 'ImmGranule1', 'Nbl2', 'nIPC', 'RadialGlia2', 'Nbl1', 'CA1-Sub', 'GlialProg']
end clusters ['CA', 'Granule', 'CA2-3-4', 'RadialGlia', 'ImmAstro', 'ImmGranule2', 'OPC', 'ImmGranule1', 'Nbl2', 'nIPC', 'RadialGlia2', 'Nbl1', 'CA1-Sub']
[6]:
fig, axes = plt.subplots(1, 3, figsize=(20, 5))
phlower.pl.persisitence_barcodes(adata, ax=axes[0], include_holes=False, barcodes_dim=[0])
phlower.pl.persisitence_barcodes(adata, ax=axes[1], include_holes=False, barcodes_dim=[0,1])
phlower.pl.persisitence_birth_death(adata, ax=axes[2], include_holes=False)
[6]:
<Axes: title={'center': 'Persistence diagram'}, xlabel='Birth', ylabel='Death'>
[7]:
fig, ax = plt.subplots(1,2, figsize=(12,4))
phlower.pl.nxdraw_group(adata, layout_name='X_pca_ddhodge_g', graph_name='X_pca_ddhodge_g_triangulation', node_size=1, show_edges=True, ax=ax[0], show_legend=False, labelstyle='text', labelsize=8)
phlower.pl.nxdraw_group(adata, layout_name='X_pca_ddhodge_g', graph_name='X_pca_ddhodge_g_triangulation_circle', node_size=2, show_edges=True, ax=ax[1], show_legend=True, labelstyle='text', labelsize=8)
[8]:
phlower.pl.nxdraw_score(adata, layout_name='X_pca_ddhodge_g', graph_name='X_pca_ddhodge_g_triangulation', node_size=1)
Hodge decomposition
[9]:
from datetime import datetime
dt1 = datetime.now()
print('L1 Datetime Start:', dt1)
phlower.tl.L1Norm_decomp(adata, L1_mode='sym', check_symmetric=False, isnorm = True)
dt2 = datetime.now()
print('L1 Datetime End:', dt2)
delta = dt2 - dt1
print('L1 Difference is:', delta)
L1 Datetime Start: 2025-01-14 17:34:46.617858
3599.1217019557953 sec
L1 Datetime End: 2025-01-14 18:34:47.351646
L1 Difference is: 1:00:00.733788
determine harmonic
[10]:
_, ax = plt.subplots(1,2, figsize=(9,3))
phlower.pl.plot_eigen_line(adata, n_eig=10,ax=ax[0])
phlower.pl.plot_eigen_line(adata, n_eig=100, step_size=10, ax=ax[1])
[11]:
phlower.tl.knee_eigen(adata)
knee eigen value is 50
preference random walk
[12]:
phlower.tl.random_climb_knn(adata, n=10000, roots_ratio=0.001, knn_edges_k=16, knn_type='diffusion')
100%|██████████| 10000/10000 [01:38<00:00, 101.81it/s]
Project paths to harmonic space
[13]:
phlower.tl.trajs_matrix(adata)
2025-01-14 18:38:49.683893 projecting trajectories to simplics...
2025-01-14 18:39:40.052710 Embedding trajectory harmonics...
eigen_n < 1, use knee_eigen to find the number of eigen vectors to use: 50
/home/sz753404/miniconda3/envs/R422/lib/python3.10/site-packages/umap/umap_.py:1943: UserWarning: n_jobs value -1 overridden to 1 by setting random_state. Use no seed for parallelism.
warn(f"n_jobs value {self.n_jobs} overridden to 1 by setting random_state. Use no seed for parallelism.")
2025-01-14 18:40:24.916227 done.
Paths clustering
[14]:
phlower.tl.trajs_clustering(adata, eps=0.3, oname_basis='')
[15]:
Counter(adata.uns['trajs_clusters'])
[15]:
Counter({0: 8297, 1: 1044, 2: 484, 4: 97, -1: 49, 3: 19, 5: 10})
[16]:
phlower.tl.unique_trajectory_clusters(adata,verbose=False)
[17]:
phlower.pl.plot_trajs_embedding(adata, clusters="trajs_clusters")
[18]:
fig, ax = phlower.pl.plot_density_grid(adata, figsize=(10, 7), cluster_name='trajs_clusters', bg_alpha=0.05, return_fig=True)
z_min, z_max 4.738726832774441e-09 1.0701098323429395e-06
Tree creating
[19]:
phlower.tl.harmonic_stream_tree(adata, verbose=True, min_bin_number=20, min_kde_quant_rm=0.01, trim_end=False)
2025-01-14 18:41:01.634833 trajectory groups ranking...
eigen_n < 1, use knee_eigen to find the number of eigen vectors to use
eigen_n: 50
{0, 1, 2, 3}
traj bins: 100%|██████████| 4/4 [00:58<00:00, 14.57s/it]
2025-01-14 18:42:49.353899 sync bins to by pseudo time...
2025-01-14 18:42:49.765746 calculating distances between bin pairs...
2025-01-14 18:42:50.011717 fate tree creating...
[(0, 1, 2, 3)]
2025-01-14 18:42:50.012303 tree attributes adding...
2025-01-14 18:42:51.286191 plugin to STREAM...
2025-01-14 18:42:51.290195 projecting trajs to harmonic...
2025-01-14 18:43:36.229265 cumsum...
2025-01-14 18:43:39.484656 edge cumsum dict...
edges projection: 100%|██████████| 9844/9844 [00:42<00:00, 234.34it/s]
2025-01-14 18:44:33.938307 done...
there are 1444 nodes not visited in the tree
will approximate the node cumsum coordinate by 5 nearest accessible neighbors
2025-01-14 18:45:07.457487 done...
stream plot
[20]:
phlower.ext.plot_stream(adata, fig_size=(8,5), dist_scale=1, log_scale=True, factor_min_win=1, factor_num_win=5, show_legend=True)
[21]:
pickle.dump(adata, open("save/hip_phlower_final.pickle", 'wb'))
[22]:
import session_info
session_info.show()
[22]:
Click to view session information
----- anndata 0.9.2 matplotlib 3.9.1.post1 numpy 1.23.5 pandas 2.2.3 phlower 0.1.3 scanpy 1.9.3 session_info 1.0.0 -----
Click to view modules imported as dependencies
Cython 0.29.33 PIL 9.4.0 anyio NA argcomplete NA arrow 1.2.3 asciitree NA asttokens NA attr 22.2.0 attrs 22.2.0 babel 2.12.1 backcall 0.2.0 brotli NA certifi 2024.08.30 cffi 1.15.1 charset_normalizer 3.1.0 cloudpickle 2.2.1 colorama 0.4.6 colorcet 3.0.1 comm 0.2.2 cycler 0.12.1 cython 0.29.33 cython_runtime NA dask 2024.6.2 dateutil 2.8.2 debugpy 1.6.6 decorator 5.1.1 defusedxml 0.7.1 dill 0.3.6 dot_parser NA exceptiongroup 1.1.3 executing 1.2.0 fastjsonschema NA fqdn NA gudhi 3.10.1 h5py 3.10.0 idna 3.4 igraph 0.10.5 importlib_metadata NA ipykernel 6.21.3 ipython_genutils 0.2.0 ipywidgets 8.1.3 isoduration NA jedi 0.19.0 jinja2 3.1.2 joblib 1.3.2 json5 NA jsonpointer 2.3 jsonschema 4.17.3 jsonschema_specifications NA jupyter_events 0.6.3 jupyter_server 2.4.0 jupyterlab_server 2.20.0 kiwisolver 1.4.5 leidenalg 0.9.1 llvmlite 0.39.1 louvain 0.8.0 markupsafe 2.1.3 matplotlib_inline 0.1.6 mpl_toolkits NA mpmath 1.3.0 msgpack 1.0.5 natsort 8.4.0 nbformat 5.9.2 networkx 2.8.8 numba 0.56.4 numcodecs 0.12.1 numexpr 2.8.4 nvfuser NA opt_einsum v3.3.0 packaging 23.1 parso 0.8.3 patsy 0.5.3 pexpect 4.8.0 pickleshare 0.7.5 pkg_resources NA platformdirs 3.8.1 plotly 5.23.0 prometheus_client NA prompt_toolkit 3.0.39 psutil 5.9.4 ptyprocess 0.7.0 pure_eval 0.2.2 pyarrow 16.1.0 pycparser 2.21 pydev_ipython NA pydevconsole NA pydevd 2.9.5 pydevd_file_utils NA pydevd_plugins NA pydevd_tracing NA pydot 1.4.2 pygments 2.15.1 pynndescent 0.5.11 pyparsing 3.1.2 pythonjsonlogger NA pytz 2023.3 referencing NA requests 2.29.0 rfc3339_validator 0.1.4 rfc3986_validator 0.1.1 rpds NA scipy 1.14.0 seaborn 0.13.2 send2trash NA setuptools_scm NA six 1.16.0 sklearn 1.5.1 sksparse 0.4.8 sniffio 1.3.0 socks 1.7.1 sparse 0.14.0 sphinxcontrib NA stack_data 0.6.2 statsmodels 0.14.0 sympy 1.11.1 tblib 3.0.0 texttable 1.6.7 threadpoolctl 3.4.0 tlz 0.12.0 toolz 0.12.0 torch 2.0.1+cu117 tornado 6.3 tqdm 4.65.0 traitlets 5.9.0 typing_extensions NA umap 0.5.5 uri_template NA urllib3 1.26.15 wcwidth 0.2.6 webcolors 1.11.1 websocket 1.5.1 yaml 6.0 zarr 2.18.2 zipp NA zmq 25.0.2
----- IPython 8.15.0 jupyter_client 8.2.0 jupyter_core 5.3.0 jupyterlab 3.6.1 notebook 6.5.3 ----- Python 3.10.8 | packaged by conda-forge | (main, Nov 22 2022, 08:23:14) [GCC 10.4.0] Linux-4.18.0-477.21.1.el8_8.x86_64-x86_64-with-glibc2.28 ----- Session information updated at 2025-01-14 18:45