# Developing-Heart-Local-Improved - Supplementary

**Pages:** 9

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Supplementary_materials/Macrophages/Mo-MP_no8_YSMP-regulon-markers.html

**Contents:**
- Variables¶
- Read in anndata object¶
- Plotting¶

**Examples:**

Example 1 (python):
```python
import scanpy as sc
import pandas as pd
import numpy as np
import anndata
import scipy
import matplotlib.pyplot as plt
import seaborn as sns

import matplotlib as mpl
# enables correct plotting of text
mpl.rcParams['pdf.fonttype'] = 42
mpl.rcParams['ps.fonttype'] = 42
```

Example 2 (python):
```python
import session_info
session_info.show()
```

Example 3 (python):
```python
-----
anndata             0.9.2
matplotlib          3.7.2
numpy               1.24.4
pandas              1.5.3
scanpy              1.9.3
scipy               1.10.1
seaborn             0.12.2
session_info        1.0.0
-----
```

Example 4 (python):
```python
PIL                         10.0.0
anyio                       NA
arrow                       1.2.3
asttokens                   NA
attr                        23.1.0
attrs                       23.1.0
babel                       2.12.1
backcall                    0.2.0
brotli                      1.0.9
certifi                     2023.11.17
cffi                        1.15.1
charset_normalizer          3.2.0
colorama                    0.4.6
comm                        0.1.4
cycler                      0.10.0
cython_runtime              NA
dateutil                    2.8.2
debugpy                     1.6.8
decorator                   5.1.1
defusedxml                  0.7.1
executing                   1.2.0
fastjsonschema              NA
fqdn                        NA
h5py                        3.9.0
idna                        3.4
igraph                      0.10.6
importlib_metadata          NA
importlib_resources         NA
ipykernel                   6.25.1
isoduration                 NA
jedi                        0.19.0
jinja2                      3.1.2
joblib                      1.3.2
json5                       NA
jsonpointer                 2.0
jsonschema                  4.19.0
jsonschema_specifications   NA
jupyter_events              0.7.0
jupyter_server              2.7.1
jupyterlab_server           2.24.0
kiwisolver                  1.4.5
leidenalg                   0.10.1
llvmlite                    0.40.1
markupsafe                  2.1.3
mpl_toolkits                NA
natsort                     8.4.0
nbformat                    5.9.2
numba                       0.57.1
overrides                   NA
packaging                   23.1
parso                       0.8.3
patsy                       0.5.3
pexpect                     4.8.0
pickleshare                 0.7.5
pkg_resources               NA
plotly                      5.16.1
prometheus_client           NA
prompt_toolkit              3.0.39
psutil                      5.9.5
ptyprocess                  0.7.0
pure_eval                   0.2.2
pydev_ipython               NA
pydevconsole                NA
pydevd                      2.9.5
pydevd_file_utils           NA
pydevd_plugins              NA
pydevd_tracing              NA
pygments                    2.16.1
pyparsing                   3.0.9
pythonjsonlogger            NA
pytz                        2023.3
referencing                 NA
requests                    2.31.0
rfc3339_validator           0.1.4
rfc3986_validator           0.1.1
rpds                        NA
ruamel                      NA
send2trash                  NA
six                         1.16.0
sklearn                     1.3.0
sniffio                     1.3.0
socks                       1.7.1
sphinxcontrib               NA
stack_data                  0.6.2
statsmodels                 0.14.0
texttable                   1.6.7
threadpoolctl               3.2.0
tornado                     6.3.3
traitlets                   5.9.0
typing_extensions           NA
uri_template                NA
urllib3                     2.0.4
wcwidth                     0.2.6
webcolors                   1.13
websocket                   1.6.2
yaml                        6.0.1
zipp                        NA
zmq                         25.1.1
zstandard                   0.19.0
```

Example 5 (unknown):
```unknown
-----
IPython             8.12.0
jupyter_client      8.3.0
jupyter_core        4.12.0
jupyterlab          4.0.5
notebook            7.0.2
-----
Python 3.8.17 (default, Jul  5 2023, 21:04:15) [GCC 11.2.0]
Linux-5.4.0-169-generic-x86_64-with-glibc2.17
-----
Session information updated at 2024-02-05 21:15
```

Example 6 (python):
```python
sc.settings.set_figure_params(dpi=120)
```

Example 7 (python):
```python
import importlib.util
import sys
spec = importlib.util.spec_from_file_location("module.name", "/nfs/team205/kk18/function/python/utils.py")
utils = importlib.util.module_from_spec(spec)
sys.modules["module.name"] = utils
spec.loader.exec_module(utils)
```

Example 8 (python):
```python
sc.settings.set_figure_params(dpi=80,dpi_save=300)
```

Example 9 (python):
```python
figdir = '/nfs/team205/heart/fetal_ms_figures/raw_plots_kk'
adata_path_dir = '/nfs/team205/heart/anndata_objects/Foetal/Immune/MoMP_Aug2023'
```

Example 10 (python):
```python
adata = sc.read_h5ad(f'{adata_path_dir}/MultiVI-post_foetal_Mono-MP_withContCov-CCCov_13Mar2023_feature-select-0003.h5ad')
adata
```

Example 11 (python):
```python
AnnData object with n_obs × n_vars = 7674 × 142283
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score', '_indices', '_scvi_batch', '_scvi_labels', 'fine_grain_withcycling', 't', 'seg', 'edge', 't_sd', 'milestones', 'fine_grain_mod'
    var: 'modality', 'gene_name_scRNA-0', 'gene_id', 'n_cells'
    uns: '_scvi_manager_uuid', '_scvi_uuid', 'adata_auc_gene', 'adata_auc_region', 'adata_gex', 'adata_gex_lognorm', 'adata_peak_tfidf', 'draw_graph', 'fine_grain_colors', 'fine_grain_mod_colors', 'fine_grain_withcycling_colors', 'graph', 'milestones_colors', 'neighbors', 'ppt', 'pseudotime_list', 'seg_colors', 'umap', 'week_colors'
    obsm: 'MultiVI_latent', 'X_R', 'X_draw_graph_fa', 'X_palantir', 'X_umap', 'X_umap13', '_scvi_extra_categorical_covs', '_scvi_extra_continuous_covs'
    obsp: 'connectivities', 'distances'
```

Example 12 (python):
```python
# get regulon score data, gene expression based
adata_auc = adata.uns['adata_auc_gene'].copy()
adata_auc.uns['draw_graph'] = adata.uns['draw_graph']
print(adata_auc.X.data[:5])
adata_auc
```

Example 13 (unknown):
```unknown
<memory at 0x151287d85380>
```

Example 14 (unknown):
```unknown
AnnData object with n_obs × n_vars = 1899 × 280
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score', '_indices', '_scvi_batch', '_scvi_labels', 'fine_grain_withcycling', 't', 'seg', 'edge', 't_sd', 'milestones'
    uns: 'draw_graph'
    obsm: 'MultiVI_latent', 'X_R', 'X_draw_graph_fa', 'X_palantir', 'X_umap', 'X_umap13', '_scvi_extra_categorical_covs', '_scvi_extra_continuous_covs'
```

Example 15 (python):
```python
markers = ['FLI1_+','MEF2C_+']
markers = [x for x in markers if x in adata_auc.var_names]

plt.rcParams['figure.figsize'] = 1.7,3.5
sc.pl.draw_graph(adata_auc,color=markers,cmap='inferno',wspace=0.45,frameon=False,show=False)
plt.savefig(f"{figdir}/Force-directed-graph_MoMP_YSMP-marker-regulon-scores.pdf", bbox_inches="tight",dpi=300)
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Supplementary_materials/Macrophages/Mo-MP_no6_MP-LYVE1_fitting-features.html

**Contents:**
- Inputs¶
- Read in trajectory object¶
- Show trajectories¶
- Subset populations¶
- Extract gene features and log-normalise¶
- Add normalised peak counts¶
- Read in scenicplus output and create anndata¶
- Differential test¶
- Save anndata¶
- Fitting¶
- Regulon - auc gene¶
- Regulon - auc region¶
- Gene¶
- Region¶
- Save result of differential and fitting result¶

1960 rows × 6 columns

1960 rows × 6 columns

884527 rows × 6 columns

**Examples:**

Example 1 (python):
```python
import scanpy as sc
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import anndata
import scipy
import os

import matplotlib as mpl
# enables correct plotting of text
mpl.rcParams['pdf.fonttype'] = 42
mpl.rcParams['ps.fonttype'] = 42
```

Example 2 (python):
```python
import session_info
session_info.show()
```

Example 3 (python):
```python
-----
anndata             0.9.2
matplotlib          3.7.2
numpy               1.24.4
pandas              1.5.3
scanpy              1.9.3
scipy               1.10.1
session_info        1.0.0
-----
```

Example 4 (python):
```python
PIL                         10.0.0
anyio                       NA
arrow                       1.2.3
asttokens                   NA
attr                        23.1.0
attrs                       23.1.0
babel                       2.12.1
backcall                    0.2.0
brotli                      1.0.9
certifi                     2023.11.17
cffi                        1.15.1
charset_normalizer          3.2.0
colorama                    0.4.6
comm                        0.1.4
cycler                      0.10.0
cython_runtime              NA
dateutil                    2.8.2
debugpy                     1.6.8
decorator                   5.1.1
defusedxml                  0.7.1
executing                   1.2.0
fastjsonschema              NA
fqdn                        NA
h5py                        3.9.0
idna                        3.4
igraph                      0.10.6
importlib_metadata          NA
importlib_resources         NA
ipykernel                   6.25.1
isoduration                 NA
jedi                        0.19.0
jinja2                      3.1.2
joblib                      1.3.2
json5                       NA
jsonpointer                 2.0
jsonschema                  4.19.0
jsonschema_specifications   NA
jupyter_events              0.7.0
jupyter_server              2.7.1
jupyterlab_server           2.24.0
kiwisolver                  1.4.5
leidenalg                   0.10.1
llvmlite                    0.40.1
markupsafe                  2.1.3
mpl_toolkits                NA
natsort                     8.4.0
nbformat                    5.9.2
numba                       0.57.1
overrides                   NA
packaging                   23.1
parso                       0.8.3
pexpect                     4.8.0
pickleshare                 0.7.5
pkg_resources               NA
plotly                      5.16.1
prometheus_client           NA
prompt_toolkit              3.0.39
psutil                      5.9.5
ptyprocess                  0.7.0
pure_eval                   0.2.2
pydev_ipython               NA
pydevconsole                NA
pydevd                      2.9.5
pydevd_file_utils           NA
pydevd_plugins              NA
pydevd_tracing              NA
pygments                    2.16.1
pyparsing                   3.0.9
pythonjsonlogger            NA
pytz                        2023.3
referencing                 NA
requests                    2.31.0
rfc3339_validator           0.1.4
rfc3986_validator           0.1.1
rpds                        NA
ruamel                      NA
send2trash                  NA
six                         1.16.0
sklearn                     1.3.0
sniffio                     1.3.0
socks                       1.7.1
sphinxcontrib               NA
stack_data                  0.6.2
texttable                   1.6.7
threadpoolctl               3.2.0
tornado                     6.3.3
traitlets                   5.9.0
typing_extensions           NA
uri_template                NA
urllib3                     2.0.4
wcwidth                     0.2.6
webcolors                   1.13
websocket                   1.6.2
yaml                        6.0.1
zipp                        NA
zmq                         25.1.1
zstandard                   0.19.0
```

Example 5 (unknown):
```unknown
-----
IPython             8.12.0
jupyter_client      8.3.0
jupyter_core        4.12.0
jupyterlab          4.0.5
notebook            7.0.2
-----
Python 3.8.17 (default, Jul  5 2023, 21:04:15) [GCC 11.2.0]
Linux-5.4.0-169-generic-x86_64-with-glibc2.17
-----
Session information updated at 2024-02-06 16:16
```

Example 6 (python):
```python
sc.settings.set_figure_params(dpi=120)
```

Example 7 (python):
```python
import importlib.util
import sys
spec = importlib.util.spec_from_file_location("module.name", "/nfs/team205/kk18/function/python/utils.py")
utils = importlib.util.module_from_spec(spec)
sys.modules["module.name"] = utils
spec.loader.exec_module(utils)
```

Example 8 (unknown):
```unknown
figdir = '/nfs/team205/heart/fetal_ms_figures/raw_plots_kk'
```

Example 9 (python):
```python
adata_path = '/nfs/team205/heart/anndata_objects/Foetal/Immune/MoMP_Aug2023/MultiVI-post_foetal_Mono-MP_withContCov-CCCov_13Mar2023_feature-select-0003.h5ad'

# methods to subset populations
key_of_groups = 'fine_grain_mod'
# groups = ['VentricularCardiomyocytesLeftCompact','VentricularCardiomyocytesLeftTrabeculated']
```

Example 10 (python):
```python
adata = sc.read_h5ad(adata_path)

# rename obs_names
adata.obs_names = [x.replace('_expression','').replace('_paired','') for x in adata.obs_names]

adata
```

Example 11 (python):
```python
AnnData object with n_obs × n_vars = 7674 × 142283
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score', '_indices', '_scvi_batch', '_scvi_labels', 'fine_grain_withcycling', 't', 'seg', 'edge', 't_sd', 'milestones', 'fine_grain_mod'
    var: 'modality', 'gene_name_scRNA-0', 'gene_id', 'n_cells'
    uns: '_scvi_manager_uuid', '_scvi_uuid', 'adata_auc_gene', 'adata_auc_region', 'adata_gex', 'adata_gex_lognorm', 'adata_peak_tfidf', 'draw_graph', 'fine_grain_colors', 'fine_grain_mod_colors', 'fine_grain_withcycling_colors', 'graph', 'milestones_colors', 'neighbors', 'ppt', 'pseudotime_list', 'seg_colors', 'umap', 'week_colors'
    obsm: 'MultiVI_latent', 'X_R', 'X_draw_graph_fa', 'X_palantir', 'X_umap', 'X_umap13', '_scvi_extra_categorical_covs', '_scvi_extra_continuous_covs'
    obsp: 'connectivities', 'distances'
```

Example 12 (python):
```python
adata.obs[key_of_groups].value_counts()
```

Example 13 (unknown):
```unknown
MacrophagesLYVE1pos     4757
MacrophagesATF3pos       898
MacrophagesCX3CR1pos     804
Monocytes                496
MacrophagesYSderived     207
MonocytesMPOpos          202
MonocyteDerivedCells     199
MacrophagesTIMD4pos      111
Name: fine_grain_mod, dtype: int64
```

Example 14 (python):
```python
sc.pl.draw_graph(adata,color=['t',key_of_groups],color_map="viridis", wspace=0.2)
```

Example 15 (unknown):
```unknown
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
```

Example 16 (python):
```python
set(adata.obs['seg'])
```

Example 17 (unknown):
```unknown
{'1', '2', '3', '4', '5', '6', '7', '8', '9'}
```

Example 18 (python):
```python
from PIL import Image
image = Image.open('./MoMP_trajectory.png')

plt.figure(figsize=(5, 5))
plt.imshow(image)
plt.axis('off')
plt.grid(False)
plt.show()
```

Example 19 (python):
```python
plt.rcParams['figure.figsize'] = 1.5,4
segments_to_select = ['3', '4', '5', '6', '7', '8', '9']
adata = adata[adata.obs['seg'].isin(segments_to_select)]
sc.pl.draw_graph(adata,color=['t',key_of_groups],color_map="viridis", wspace=0.5)
```

Example 20 (unknown):
```unknown
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
```

Example 21 (python):
```python
# scale t
# min-max scaling
adata.obs['t_scaled'] = (adata.obs['t']-adata.obs['t'].min())/(adata.obs['t'].max()-adata.obs['t'].min())

plt.rcParams['figure.figsize'] = 1.7,3.5
sc.pl.draw_graph(adata,color=['t_scaled'],color_map="viridis",frameon=False,
                 title='MP maturation\n(pseudotime)', show=False
                )
plt.savefig(f"{figdir}/Force-directed-graph_MP_pseudotime.pdf", bbox_inches="tight",dpi=300)
```

Example 22 (python):
```python
/tmp/ipykernel_203412/3235555261.py:3: ImplicitModificationWarning: Trying to modify attribute `.obs` of view, initializing view as actual.
  adata.obs['t_scaled'] = (adata.obs['t']-adata.obs['t'].min())/(adata.obs['t'].max()-adata.obs['t'].min())
```

Example 23 (python):
```python
adata.obs['week']
```

Example 24 (unknown):
```unknown
BHF_F_Hea11933667_BHF_F_Hea11596620_TCATGCCTCGTTATCT-1    13W
BHF_F_Hea11933667_BHF_F_Hea11596620_TCTAATCTCCTCACTA-1    13W
7089STDY13216927_BHF_F_Hea13242534_GTTTCTAGTCCTAGTT-1     12W
BHF_F_Hea11933671_BHF_F_Hea11596624_GCTAGCCAGCACTTGG-1     5W
BHF_F_Hea11933668_BHF_F_Hea11596621_GTTTCTAGTTAGCTGA-1    13W
                                                         ... 
BHF_F_Hea11192328_AAGGTTCTCAGAGGTG-1                       8W
BHF_F_Hea11192328_AAGACCTCACACGCTG-1                       8W
BHF_F_Hea11192328_CATATGGCAGCTCGCA-1                       8W
BHF_F_Hea11192328_CAACCAACACGCGAAA-1                       8W
BHF_F_Hea11192328_GGATGTTGTGTTTGTG-1                       8W
Name: week, Length: 3886, dtype: category
Categories (7, object): ['5W', '7W', '8W', '9W', '10W', '12W', '13W']
```

Example 25 (python):
```python
# convert to continuous week
adata.obs['week_int'] = [int(x.replace('W','')) for x in adata.obs['week']]

plt.rcParams['figure.figsize'] = 1.7,3.5
sc.pl.draw_graph(adata,color=['week_int'],color_map="viridis",frameon=False,
                 title='PCW', show=False
                )
plt.savefig(f"{figdir}/Force-directed-graph_MP_PCW.pdf", bbox_inches="tight",dpi=300)
```

Example 26 (python):
```python
# extract gene features
adata_gex = adata.copy()
gene_vars = adata_gex.var_names[adata_gex.var['modality']=='Gene Expression']
adata_gex = adata_gex[:,gene_vars]

# log normalise 
print(adata_gex.X.data[:10])
sc.pp.normalize_total(adata_gex, target_sum=1e4)
sc.pp.log1p(adata_gex)
print(adata_gex.X.data[:10])

# store to .uns
adata.uns['adata_gex_lognorm']=adata_gex.copy()
adata.uns['adata_gex_lognorm']
```

Example 27 (unknown):
```unknown
[1. 7. 1. 4. 1. 4. 1. 1. 9. 5.]
```

Example 28 (python):
```python
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/scanpy/preprocessing/_normalization.py:170: UserWarning: Received a view of an AnnData. Making a copy.
  view_to_actual(adata)
```

Example 29 (unknown):
```unknown
[1.0841573 2.6877544 1.0841573 2.1779046 1.0841573 2.1779046 1.0841573
 1.0841573 2.9238348 2.3781319]
```

Example 30 (python):
```python
/tmp/ipykernel_150291/1223271433.py:13: ImplicitModificationWarning: Trying to modify attribute `._uns` of view, initializing view as actual.
  adata.uns['adata_gex_lognorm']=adata_gex.copy()
```

Example 31 (python):
```python
AnnData object with n_obs × n_vars = 3886 × 15922
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score', '_indices', '_scvi_batch', '_scvi_labels', 'fine_grain_withcycling', 't', 'seg', 'edge', 't_sd', 'milestones', 'fine_grain_mod'
    var: 'modality', 'gene_name_scRNA-0', 'gene_id', 'n_cells'
    uns: '_scvi_manager_uuid', '_scvi_uuid', 'adata_auc_gene', 'adata_auc_region', 'adata_gex', 'adata_gex_lognorm', 'adata_peak_tfidf', 'draw_graph', 'fine_grain_colors', 'fine_grain_mod_colors', 'fine_grain_withcycling_colors', 'graph', 'milestones_colors', 'neighbors', 'ppt', 'pseudotime_list', 'seg_colors', 'umap', 'week_colors', 'log1p'
    obsm: 'MultiVI_latent', 'X_R', 'X_draw_graph_fa', 'X_palantir', 'X_umap', 'X_umap13', '_scvi_extra_categorical_covs', '_scvi_extra_continuous_covs'
    obsp: 'connectivities', 'distances'
```

Example 32 (python):
```python
adata_peak = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/PeakMatrix/Foetal_Peaks_tfidf.h5ad')

shared_cell_ids = list(set(adata.obs_names).intersection(adata_peak.obs_names))

# subset cells and regions
adata.uns['adata_peak_tfidf']=adata_peak[shared_cell_ids,adata.var_names[adata.var['modality']=='Peaks']]
adata_peak = adata.uns['adata_peak_tfidf'].copy()
adata_peak
```

Example 33 (unknown):
```unknown
AnnData object with n_obs × n_vars = 1162 × 126361
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels'
    var: 'seqnames', 'start', 'end', 'width', 'strand', 'score', 'replicateScoreQuantile', 'groupScoreQuantile', 'Reproducibility', 'GroupReplicate', 'distToGeneStart', 'nearestGene', 'peakType', 'distToTSS', 'nearestTSS', 'GC', 'idx', 'N'
```

Example 34 (python):
```python
# add milestones and segments
columns_to_add = ['t', 'seg', 'edge', 't_sd', 'milestones']
adata.uns['adata_peak_tfidf'].obs[columns_to_add] = adata.obs.loc[adata.uns['adata_peak_tfidf'].obs_names, columns_to_add]
adata_peak = adata.uns['adata_peak_tfidf'].copy()
adata_peak
```

Example 35 (python):
```python
/tmp/ipykernel_150291/281467255.py:3: ImplicitModificationWarning: Trying to modify attribute `.obs` of view, initializing view as actual.
  adata.uns['adata_peak_tfidf'].obs[columns_to_add] = adata.obs.loc[adata.uns['adata_peak_tfidf'].obs_names, columns_to_add]
```

Example 36 (unknown):
```unknown
AnnData object with n_obs × n_vars = 1162 × 126361
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 't', 'seg', 'edge', 't_sd', 'milestones'
    var: 'seqnames', 'start', 'end', 'width', 'strand', 'score', 'replicateScoreQuantile', 'groupScoreQuantile', 'Reproducibility', 'GroupReplicate', 'distToGeneStart', 'nearestGene', 'peakType', 'distToTSS', 'nearestTSS', 'GC', 'idx', 'N'
```

Example 37 (python):
```python
# read in regulon score
auc_gene = pd.read_csv('/nfs/team205/heart/pyscenic/plus/AUCscore/MoMP_eRegulon_AUC_gene.csv',index_col=0)
print(auc_gene.shape)
auc_region = pd.read_csv('/nfs/team205/heart/pyscenic/plus/AUCscore/MoMP_eRegulon_AUC_region.csv',index_col=0)
print(auc_region.shape)

# select shared barcodes (ie only multiome data, only monocyte/macrophage)
shared_cell_ids = list(set(auc_gene.index).intersection(auc_region.index).intersection(adata.obs_names))
print(len(shared_cell_ids))
```

Example 38 (unknown):
```unknown
(1899, 280)
(1899, 280)
1162
```

Example 39 (python):
```python
# gene
# select cells in anndata (monocyte and macrophage)
auc_gene = auc_gene.loc[shared_cell_ids]

# create anndata
adata_aucg = anndata.AnnData(
    X=auc_gene,
    obs=adata[shared_cell_ids].obs,
    var=pd.DataFrame(auc_gene.columns).set_index(0),
    obsm=adata[shared_cell_ids].obsm
    # obsp=adata.obsp.copy(),
    # uns=adata.uns.copy(),
)
# modify varnames
adata_aucg.var_names = [x.split('_(')[0] for x in adata_aucg.var_names]

# store to .uns
adata.uns['adata_auc_gene']=adata_aucg.copy()
adata.uns['adata_auc_gene']
```

Example 40 (unknown):
```unknown
AnnData object with n_obs × n_vars = 1162 × 280
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score', '_indices', '_scvi_batch', '_scvi_labels', 'fine_grain_withcycling', 't', 'seg', 'edge', 't_sd', 'milestones', 'fine_grain_mod'
    obsm: 'MultiVI_latent', 'X_R', 'X_draw_graph_fa', 'X_palantir', 'X_umap', 'X_umap13', '_scvi_extra_categorical_covs', '_scvi_extra_continuous_covs'
```

Example 41 (python):
```python
# region
# select cells in anndata (monocyte and macrophage)
auc_region = auc_region.loc[shared_cell_ids]

# create anndata
adata_aucr = anndata.AnnData(
    X=auc_region,
    obs=adata[shared_cell_ids].obs,
    var=pd.DataFrame(auc_region.columns).set_index(0),
    obsm=adata[shared_cell_ids].obsm
    # obsp=adata.obsp.copy(),
    # uns=adata.uns.copy(),
)
# modify varnames
adata_aucr.var_names = [x.split('_(')[0] for x in adata_aucr.var_names]

# store to .uns
adata.uns['adata_auc_region']=adata_aucr.copy()
adata.uns['adata_auc_region']
```

Example 42 (unknown):
```unknown
AnnData object with n_obs × n_vars = 1162 × 280
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score', '_indices', '_scvi_batch', '_scvi_labels', 'fine_grain_withcycling', 't', 'seg', 'edge', 't_sd', 'milestones', 'fine_grain_mod'
    obsm: 'MultiVI_latent', 'X_R', 'X_draw_graph_fa', 'X_palantir', 'X_umap', 'X_umap13', '_scvi_extra_categorical_covs', '_scvi_extra_continuous_covs'
```

Example 43 (unknown):
```unknown
root_milestone = '29'
```

Example 44 (python):
```python
# prepare a dictionary to put result
diff_dict = {}
diff_dict['auc_gene']=np.nan
diff_dict['auc_region']=np.nan
diff_dict['gene']=np.nan
diff_dict['region']=np.nan
diff_dict
```

Example 45 (unknown):
```unknown
{'auc_gene': nan, 'auc_region': nan, 'gene': nan, 'region': nan}
```

Example 46 (python):
```python
# Differential test
# genes
sc.tl.rank_genes_groups(adata_gex,groupby='milestones',reference=root_milestone,method='t-test')
## Get result and store in .uns
DE_df=sc.get.rank_genes_groups_df(adata_gex, group=None)
diff_dict['gene']=DE_df.copy()

# regions
sc.tl.rank_genes_groups(adata_peak,groupby='milestones',reference=root_milestone,method='t-test')
## Get result and store in .uns
DE_df=sc.get.rank_genes_groups_df(adata_peak, group=None)
diff_dict['region']=DE_df.copy()

# auc-gene
sc.tl.rank_genes_groups(adata_aucg,groupby='milestones',reference=root_milestone,method='t-test')
## Get result and store in .uns
DE_df=sc.get.rank_genes_groups_df(adata_aucg, group=None)
diff_dict['auc_gene']=DE_df.copy()

# auc-region
sc.tl.rank_genes_groups(adata_aucr,groupby='milestones',reference=root_milestone,method='t-test')
## Get result and store in .uns
DE_df=sc.get.rank_genes_groups_df(adata_aucr, group=None)
diff_dict['auc_region']=DE_df.copy()
```

Example 47 (unknown):
```unknown
diff_dict['auc_gene']
```

Example 48 (unknown):
```unknown
diff_dict['auc_region']
```

Example 49 (unknown):
```unknown
diff_dict['region']
```

Example 50 (unknown):
```unknown
segments_to_select.sort()
outdir = f'/nfs/team205/heart/anndata_objects/Foetal/Immune/MoMP_Aug2023/seg_{"-".join(segments_to_select)}'
outdir
```

Example 51 (unknown):
```unknown
'/nfs/team205/heart/anndata_objects/Foetal/Immune/MoMP_Aug2023/seg_3-4-5-6-7-8-9'
```

Example 52 (python):
```python
# save anndata
if os.path.exists(outdir):
    print('directory already exists')
else:
    os.makedirs(outdir)
adata.write(f'{outdir}/adata.h5ad')
adata
```

Example 53 (python):
```python
AnnData object with n_obs × n_vars = 3886 × 142283
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score', '_indices', '_scvi_batch', '_scvi_labels', 'fine_grain_withcycling', 't', 'seg', 'edge', 't_sd', 'milestones', 'fine_grain_mod'
    var: 'modality', 'gene_name_scRNA-0', 'gene_id', 'n_cells'
    uns: '_scvi_manager_uuid', '_scvi_uuid', 'adata_auc_gene', 'adata_auc_region', 'adata_gex', 'adata_gex_lognorm', 'adata_peak_tfidf', 'draw_graph', 'fine_grain_colors', 'fine_grain_mod_colors', 'fine_grain_withcycling_colors', 'graph', 'milestones_colors', 'neighbors', 'ppt', 'pseudotime_list', 'seg_colors', 'umap', 'week_colors'
    obsm: 'MultiVI_latent', 'X_R', 'X_draw_graph_fa', 'X_palantir', 'X_umap', 'X_umap13', '_scvi_extra_categorical_covs', '_scvi_extra_continuous_covs'
    obsp: 'connectivities', 'distances'
```

Example 54 (python):
```python
# read in anndata
adata = sc.read_h5ad(f'{outdir}/adata.h5ad')
adata
```

Example 55 (python):
```python
import scFates as scf
```

Example 56 (python):
```python
# prepare a dictionary to put result
fit_dict = {}
fit_dict['auc_gene']=np.nan
fit_dict['auc_region']=np.nan
fit_dict['gene']=np.nan
fit_dict['region']=np.nan
fit_dict
```

Example 57 (unknown):
```unknown
{'auc_gene': nan, 'auc_region': nan, 'gene': nan, 'region': nan}
```

Example 58 (python):
```python
%%time
# fit
scf.tl.test_association(adata.uns['adata_auc_gene'],n_jobs=20)

# extract results and store to the dictionary
df = adata.uns['adata_auc_gene'].var.copy()
df['activator_repressor'] = [x.split('_')[-1] for x in df.index]
df = df.sort_values('A',ascending=False)
fit_dict['auc_gene'] = df.copy()
```

Example 59 (unknown):
```unknown
test features for association with the trajectory
    single mapping :   0%|          | 0/280 [00:00<?, ?it/s]
```

Example 60 (unknown):
```unknown
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
```

Example 61 (unknown):
```unknown
single mapping :   1%|          | 3/280 [01:01<1:34:08, 20.39s/it]
```

Example 62 (unknown):
```unknown
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
```

Example 63 (unknown):
```unknown
single mapping :   2%|▎         | 7/280 [02:18<1:30:18, 19.85s/it]
```

Example 64 (unknown):
```unknown
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/joblib/externals/loky/process_executor.py:752: UserWarning: A worker stopped while some jobs were given to the executor. This can be caused by a too short worker timeout or by a memory leak.
  warnings.warn(
```

Example 65 (unknown):
```unknown
single mapping :  69%|██████▉   | 194/280 [06:17<02:47,  1.94s/it]
```

Example 66 (unknown):
```unknown
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
```

Example 67 (python):
```python
single mapping : 100%|██████████| 280/280 [06:45<00:00,  1.45s/it]
    found 0 significant features (0:06:45) --> added
    .var['p_val'] values from statistical test.
    .var['fdr'] corrected values from multiple testing.
    .var['st'] proportion of mapping in which feature is significant.
    .var['A'] amplitue of change of tested feature.
    .var['signi'] feature is significantly changing along pseudotime.
    .uns['stat_assoc_list'] list of fitted features on the graph for all mappings.
CPU times: user 1.89 s, sys: 2.75 s, total: 4.63 s
Wall time: 6min 45s
```

Example 68 (python):
```python
%%time
# fit
scf.tl.test_association(adata.uns['adata_auc_region'],n_jobs=20)

# extract results and store to the dictionary
df = adata.uns['adata_auc_region'].var.copy()
df['activator_repressor'] = [x.split('_')[-1] for x in df.index]
df = df.sort_values('A',ascending=False)
fit_dict['auc_region'] = df.copy()
```

Example 69 (python):
```python
test features for association with the trajectory
    single mapping : 100%|██████████| 280/280 [01:41<00:00,  2.76it/s]
    found 0 significant features (0:01:41) --> added
    .var['p_val'] values from statistical test.
    .var['fdr'] corrected values from multiple testing.
    .var['st'] proportion of mapping in which feature is significant.
    .var['A'] amplitue of change of tested feature.
    .var['signi'] feature is significantly changing along pseudotime.
    .uns['stat_assoc_list'] list of fitted features on the graph for all mappings.
CPU times: user 1.87 s, sys: 620 ms, total: 2.49 s
Wall time: 1min 41s
```

Example 70 (python):
```python
%%time
# fit
scf.tl.test_association(adata.uns['adata_gex_lognorm'],n_jobs=20)

# extract results and store to the dictionary
df = adata.uns['adata_gex_lognorm'].var.copy()
df = df.sort_values('A',ascending=False)
fit_dict['gene'] = df.copy()
```

Example 71 (python):
```python
test features for association with the trajectory
    single mapping : 100%|██████████| 15922/15922 [1:56:39<00:00,  2.27it/s]
    found 1045 significant features (1:56:39) --> added
    .var['p_val'] values from statistical test.
    .var['fdr'] corrected values from multiple testing.
    .var['st'] proportion of mapping in which feature is significant.
    .var['A'] amplitue of change of tested feature.
    .var['signi'] feature is significantly changing along pseudotime.
    .uns['stat_assoc_list'] list of fitted features on the graph for all mappings.
CPU times: user 1min 31s, sys: 28.4 s, total: 2min
Wall time: 1h 56min 45s
```

Example 72 (python):
```python
%%time
# fit
scf.tl.test_association(adata.uns['adata_peak_tfidf'],n_jobs=20)

# extract results and store to the dictionary
df = adata.uns['adata_peak_tfidf'].var.copy()
df = df.sort_values('A',ascending=False)
fit_dict['region'] = df.copy()
```

Example 73 (unknown):
```unknown
test features for association with the trajectory
    single mapping :   0%|          | 0/126361 [00:00<?, ?it/s]
```

Example 74 (unknown):
```unknown
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
```

Example 75 (unknown):
```unknown
single mapping :   0%|          | 3/126361 [00:48<571:45:08, 16.29s/it]
```

Example 76 (unknown):
```unknown
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
```

Example 77 (unknown):
```unknown
single mapping :   0%|          | 127/126361 [03:54<64:45:58,  1.85s/it]
```

Example 78 (unknown):
```unknown
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
Unable to determine R library path: Command '('/software/isg/private/R-4.1.2/bin/Rscript', '-e', 'cat(Sys.getenv("LD_LIBRARY_PATH"))')' returned non-zero exit status 2.
```

Example 79 (python):
```python
single mapping : 100%|██████████| 126361/126361 [5:08:20<00:00,  6.83it/s]  
    found 700 significant features (5:08:20) --> added
    .var['p_val'] values from statistical test.
    .var['fdr'] corrected values from multiple testing.
    .var['st'] proportion of mapping in which feature is significant.
    .var['A'] amplitue of change of tested feature.
    .var['signi'] feature is significantly changing along pseudotime.
    .uns['stat_assoc_list'] list of fitted features on the graph for all mappings.
CPU times: user 6min 36s, sys: 1min 2s, total: 7min 38s
Wall time: 5h 8min 33s
```

Example 80 (unknown):
```unknown
utils.save_pkl(diff_dict, 
         file=f'{outdir}/differential-test_result.pkl')
```

Example 81 (unknown):
```unknown
utils.save_pkl(fit_dict, 
         file=f'{outdir}/fitting_result.pkl')
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Supplementary_materials/Macrophages/Mo-MP_no2_MultiVI_run-GPU_Mo-MP.html

**Contents:**
- Import modules¶
- Read in data¶
- Concatenate anndata¶
- Setup and training MultiVI (with continuous covariates)¶

1899 rows × 40 columns

7674 rows × 1 columns

544641 rows × 3 columns

MultiVI requires the main batch annotation to correspond to the modality of the samples. Other batch annotation, such as in the case of multiple RNA-only batches, can be specified using categorical_covariate_keys.

**Examples:**

Example 1 (python):
```python
import scvi
import anndata
import numpy as np
import pandas as pd
import scanpy as sc
import matplotlib.pyplot as plt

scvi.settings.seed = 420
```

Example 2 (python):
```python
Global seed set to 0
/nfs/team205/kk18/miniconda3/envs/scvi0190_env/lib/python3.9/site-packages/flax/core/frozen_dict.py:169: FutureWarning: jax.tree_util.register_keypaths is deprecated, and will be removed in a future release. Please use `register_pytree_with_keys()` instead.
  jax.tree_util.register_keypaths(
/nfs/team205/kk18/miniconda3/envs/scvi0190_env/lib/python3.9/site-packages/pytorch_lightning/utilities/warnings.py:53: LightningDeprecationWarning: pytorch_lightning.utilities.warnings.rank_zero_deprecation has been deprecated in v1.6 and will be removed in v1.8. Use the equivalent function from the pytorch_lightning.utilities.rank_zero module instead.
  new_rank_zero_deprecation(
/nfs/team205/kk18/miniconda3/envs/scvi0190_env/lib/python3.9/site-packages/pytorch_lightning/utilities/warnings.py:58: LightningDeprecationWarning: The `pytorch_lightning.loggers.base.rank_zero_experiment` is deprecated in v1.7 and will be removed in v1.9. Please use `pytorch_lightning.loggers.logger.rank_zero_experiment` instead.
  return new_rank_zero_deprecation(*args, **kwargs)
Global seed set to 420
```

Example 3 (python):
```python
import session_info
session_info.show()
```

Example 4 (python):
```python
-----
anndata             0.8.0
matplotlib          3.7.1
numpy               1.23.5
pandas              1.5.3
scanpy              1.9.3
scvi                0.19.0
session_info        1.0.0
-----
```

Example 5 (python):
```python
PIL                 9.4.0
absl                NA
asttokens           NA
attr                22.2.0
backcall            0.2.0
certifi             2022.12.07
cffi                1.15.1
charset_normalizer  3.1.0
chex                0.1.6
comm                0.1.2
contextlib2         NA
cycler              0.10.0
cython_runtime      NA
dateutil            2.8.2
debugpy             1.6.6
decorator           5.1.1
defusedxml          0.7.1
deprecate           0.3.2
docrep              0.3.2
etils               1.1.0
executing           1.2.0
flax                0.6.7
fsspec              2023.3.0
google              NA
h5py                3.8.0
idna                3.4
igraph              0.10.4
importlib_metadata  NA
importlib_resources NA
ipykernel           6.21.3
ipython_genutils    0.2.0
ipywidgets          8.0.4
jax                 0.4.6
jaxlib              0.4.6
jedi                0.18.2
joblib              1.2.0
kiwisolver          1.4.4
leidenalg           0.9.1
llvmlite            0.39.1
ml_collections      NA
mpl_toolkits        NA
msgpack             1.0.5
mudata              0.2.1
multipledispatch    0.6.0
natsort             8.3.1
numba               0.56.4
numpyro             0.11.0
nvidia              NA
opt_einsum          v3.3.0
optax               0.1.4
packaging           23.0
parso               0.8.3
pexpect             4.8.0
pickleshare         0.7.5
pkg_resources       NA
platformdirs        3.1.1
prompt_toolkit      3.0.38
psutil              5.9.4
ptyprocess          0.7.0
pure_eval           0.2.2
pydev_ipython       NA
pydevconsole        NA
pydevd              2.9.5
pydevd_file_utils   NA
pydevd_plugins      NA
pydevd_tracing      NA
pygments            2.14.0
pyparsing           3.0.9
pyro                1.8.4
pytorch_lightning   1.7.7
pytz                2022.7.1
requests            2.28.2
rich                NA
scipy               1.10.1
setuptools          65.6.3
six                 1.16.0
sklearn             1.2.2
stack_data          0.6.2
tensorboard         2.12.0
texttable           1.6.7
threadpoolctl       3.1.0
toolz               0.12.0
torch               1.13.1+cu117
torchmetrics        0.11.4
torchvision         0.14.1+cu117
tornado             6.2
tqdm                4.65.0
traitlets           5.9.0
tree                0.1.8
typing_extensions   NA
urllib3             1.26.15
wcwidth             0.2.6
yaml                6.0
zipp                NA
zmq                 25.0.1
zoneinfo            NA
```

Example 6 (unknown):
```unknown
-----
IPython             8.11.0
jupyter_client      8.0.3
jupyter_core        5.2.0
notebook            6.5.3
-----
Python 3.9.16 (main, Mar  8 2023, 14:00:05) [GCC 11.2.0]
Linux-4.15.0-200-generic-x86_64-with-glibc2.27
-----
Session information updated at 2023-03-13 21:45
```

Example 7 (python):
```python
# paired
adata_paired = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/Immune/MultiVI-pre_paired_foetal_MoMP_raw.h5ad')
adata_paired
```

Example 8 (unknown):
```unknown
AnnData object with n_obs × n_vars = 1899 × 544641
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score'
    var: 'modality', 'gene_name_scRNA-0', 'gene_id'
```

Example 9 (python):
```python
adata_paired.obs
```

Example 10 (python):
```python
sum(adata_paired.var['modality']=='Gene Expression')
```

Example 11 (python):
```python
# scrna
adata_scrna = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/Immune/MultiVI-pre_scRNA_foetal_MoMP_raw.h5ad')
adata_scrna
```

Example 12 (unknown):
```unknown
AnnData object with n_obs × n_vars = 5775 × 36601
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score'
    var: 'modality', 'gene_name_scRNA-0', 'gene_id'
```

Example 13 (python):
```python
adata_mvi = scvi.data.organize_multiome_anndatas(adata_paired, adata_scrna)
adata_mvi
```

Example 14 (python):
```python
/nfs/team205/kk18/miniconda3/envs/scvi0190_env/lib/python3.9/site-packages/anndata/_core/anndata.py:1785: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  [AnnData(sparse.csr_matrix(a.shape), obs=a.obs) for a in all_adatas],
/nfs/team205/kk18/miniconda3/envs/scvi0190_env/lib/python3.9/site-packages/anndata/_core/anndata.py:1785: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.
  [AnnData(sparse.csr_matrix(a.shape), obs=a.obs) for a in all_adatas],
```

Example 15 (unknown):
```unknown
AnnData object with n_obs × n_vars = 7674 × 544641
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score'
    var: 'modality', 'gene_name_scRNA-0', 'gene_id'
```

Example 16 (python):
```python
adata_mvi.obs[['modality']]
```

Example 17 (python):
```python
# MultiVI requires the features to be ordered, 
# such that genes appear before genomic regions. In this case this is already the case, but it’s always good to verify:
adata_mvi = adata_mvi[:, adata_mvi.var["modality"].argsort()].copy()
adata_mvi.var
```

Example 18 (python):
```python
# We also filter features to remove those that appear in fewer than 1% of the cells
# modified to 0.3% to retain around 100000 features
print(adata_mvi.shape)
sc.pp.filter_genes(adata_mvi, min_cells=int(adata_mvi.shape[0] * 0.003))
print(adata_mvi.shape)
```

Example 19 (unknown):
```unknown
(7674, 544641)
(7674, 142283)
```

Example 20 (python):
```python
adata_mvi.var['modality'].value_counts()
```

Example 21 (unknown):
```unknown
Peaks              126361
Gene Expression     15922
Name: modality, dtype: int64
```

Example 22 (python):
```python
adata_mvi.obs['modality'].value_counts()
```

Example 23 (unknown):
```unknown
expression    5775
paired        1899
Name: modality, dtype: int64
```

Example 24 (python):
```python
adata_mvi.obs['batch_key'].value_counts()
```

Example 25 (unknown):
```unknown
C99_Cell_5prime-v2          4129
C94_Cell_5prime-v2          1080
Hst33_Nuclei_Multiome-v1     838
C86_Cell_5prime-v2           566
C92_Nuclei_Multiome-v1       417
C98_Nuclei_Multiome-v1       319
Hst36_Nuclei_Multiome-v1     120
C85_Nuclei_Multiome-v1        81
C87_Nuclei_Multiome-v1        74
C82_Nuclei_Multiome-v1        50
Name: batch_key, dtype: int64
```

Example 26 (python):
```python
scvi.model.MULTIVI.setup_anndata(adata_mvi, batch_key='modality',
                                categorical_covariate_keys=['batch_key'],
                                continuous_covariate_keys=['total_counts','pct_counts_mt','pct_counts_ribo',
                                                          'S_score','G2M_score'])
```

Example 27 (unknown):
```unknown
No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)
```

Example 28 (python):
```python
# When creating the object, we need to specify how many of the features are genes, and how many are genomic regions. 
# This is so MultiVI can determine the exact architecture for each modality.
mvi = scvi.model.MULTIVI(
    adata_mvi,
    n_genes=(adata_mvi.var['modality']=='Gene Expression').sum(),
    n_regions=(adata_mvi.var['modality']=='Peaks').sum(),
)
mvi.view_anndata_setup()
```

Example 29 (unknown):
```unknown
/nfs/team205/kk18/miniconda3/envs/scvi0190_env/lib/python3.9/site-packages/torch/nn/init.py:405: UserWarning: Initializing zero-element tensors is a no-op
  warnings.warn("Initializing zero-element tensors is a no-op")
```

Example 30 (unknown):
```unknown
Anndata setup with scvi-tools version 0.19.0.
```

Example 31 (unknown):
```unknown
Setup via `MULTIVI.setup_anndata` with arguments:
```

Example 32 (unknown):
```unknown
{
│   'layer': None,
│   'batch_key': 'modality',
│   'size_factor_key': None,
│   'categorical_covariate_keys': ['batch_key'],
│   'continuous_covariate_keys': ['total_counts', 'pct_counts_mt', 'pct_counts_ribo', 'S_score', 'G2M_score'],
│   'protein_expression_obsm_key': None,
│   'protein_names_uns_key': None
}
```

Example 33 (unknown):
```unknown
Summary Statistics          
┏━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━┓
┃     Summary Stat Key     ┃ Value  ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━┩
│         n_batch          │   2    │
│         n_cells          │  7674  │
│ n_extra_categorical_covs │   1    │
│ n_extra_continuous_covs  │   5    │
│         n_labels         │   1    │
│          n_vars          │ 142283 │
└──────────────────────────┴────────┘
```

Example 34 (python):
```python
Data Registry                             
┏━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃      Registry Key      ┃            scvi-tools Location             ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│           X            │                  adata.X                   │
│         batch          │          adata.obs['_scvi_batch']          │
│ extra_categorical_covs │ adata.obsm['_scvi_extra_categorical_covs'] │
│ extra_continuous_covs  │ adata.obsm['_scvi_extra_continuous_covs']  │
│         ind_x          │           adata.obs['_indices']            │
│         labels         │         adata.obs['_scvi_labels']          │
└────────────────────────┴────────────────────────────────────────────┘
```

Example 35 (python):
```python
batch State Registry                    
┏━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓
┃    Source Location    ┃ Categories ┃ scvi-tools Encoding ┃
┡━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩
│ adata.obs['modality'] │ expression │          0          │
│                       │   paired   │          1          │
└───────────────────────┴────────────┴─────────────────────┘
```

Example 36 (python):
```python
labels State Registry                      
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓
┃      Source Location      ┃ Categories ┃ scvi-tools Encoding ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩
│ adata.obs['_scvi_labels'] │     0      │          0          │
└───────────────────────────┴────────────┴─────────────────────┘
```

Example 37 (python):
```python
batch State Registry                    
┏━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓
┃    Source Location    ┃ Categories ┃ scvi-tools Encoding ┃
┡━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩
│ adata.obs['modality'] │ expression │          0          │
│                       │   paired   │          1          │
└───────────────────────┴────────────┴─────────────────────┘
```

Example 38 (python):
```python
extra_categorical_covs State Registry                   
┏━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓
┃    Source Location     ┃        Categories        ┃ scvi-tools Encoding ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩
│ adata.obs['batch_key'] │  C82_Nuclei_Multiome-v1  │          0          │
│                        │  C85_Nuclei_Multiome-v1  │          1          │
│                        │    C86_Cell_5prime-v2    │          2          │
│                        │  C87_Nuclei_Multiome-v1  │          3          │
│                        │  C92_Nuclei_Multiome-v1  │          4          │
│                        │    C94_Cell_5prime-v2    │          5          │
│                        │  C98_Nuclei_Multiome-v1  │          6          │
│                        │    C99_Cell_5prime-v2    │          7          │
│                        │ Hst33_Nuclei_Multiome-v1 │          8          │
│                        │ Hst36_Nuclei_Multiome-v1 │          9          │
│                        │                          │                     │
└────────────────────────┴──────────────────────────┴─────────────────────┘
```

Example 39 (python):
```python
extra_continuous_covs State   
            Registry            
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃       Source Location        ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│  adata.obs['total_counts']   │
│  adata.obs['pct_counts_mt']  │
│ adata.obs['pct_counts_ribo'] │
│     adata.obs['S_score']     │
│    adata.obs['G2M_score']    │
└──────────────────────────────┘
```

Example 40 (unknown):
```unknown
mvi.train()
```

Example 41 (unknown):
```unknown
GPU available: True (cuda), used: True
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs
/nfs/team205/kk18/miniconda3/envs/scvi0190_env/lib/python3.9/site-packages/pytorch_lightning/trainer/configuration_validator.py:267: LightningDeprecationWarning: The `Callback.on_epoch_end` hook was deprecated in v1.6 and will be removed in v1.8. Please use `Callback.on_<train/validation/test>_epoch_end` instead.
  rank_zero_deprecation(
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]
```

Example 42 (unknown):
```unknown
Epoch 404/500:  81%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████                              | 404/500 [52:59<12:35,  7.87s/it, loss=6.21e+03, v_num=1]
Monitored metric reconstruction_loss_validation did not improve in the last 50 records. Best score: 12617.060. Signaling Trainer to stop.
```

Example 43 (unknown):
```unknown
mvi.save("/nfs/team205/heart/anndata_objects/Foetal/Immune/MultiVI-model_foetal_Mono-MP_withContCov-CCCov_13Mar2023")
```

Example 44 (python):
```python
adata_mvi.obsm["MultiVI_latent"] = mvi.get_latent_representation()
```

Example 45 (python):
```python
adata_mvi.write('/nfs/team205/heart/anndata_objects/Foetal/Immune/MultiVI-post_foetal_Mono-MP_withContCov-CCCov_13Mar2023.h5ad')
```

Example 46 (python):
```python
# Convergency plot
mvi.history['elbo_train']
x = np.linspace(1, (len(mvi.history['elbo_train'])), (len(mvi.history['elbo_train'])))
plt.plot(x, mvi.history['elbo_train'], label="train")
plt.xlabel('epoch')
plt.ylabel('elbo_train')
plt.show()
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Supplementary_materials/Macrophages/Mo-MP_no3_scFate.html

**Contents:**
- Read in¶
- scFates - Tree analysis¶
- Preprocessing pipeline from Palantir¶
  - Create projection dataframe¶
  - Run Palantir to obtain multiscale diffusion space¶
  - Generate embedding from the multiscale diffusion space¶
  - Plotting¶
- Tree learning with SimplePPT (or ElPiGraph)¶
  - Projecting results onto ForceAtlas2 embedding¶
  - Slecting roots and computing pseudotime¶
- Representing the trajectory and tree (on top of existing embedding)¶
- Save the object¶
- Plots for figure¶

ref: https://scfates.readthedocs.io/en/latest/Tree_Analysis_Bone_marrow_fates.html

142283 rows × 4 columns

7674 rows × 18 columns

**Examples:**

Example 1 (python):
```python
import scanpy as sc
import pandas as pd
import numpy as np
import anndata
import scipy
```

Example 2 (python):
```python
import session_info
session_info.show()
```

Example 3 (unknown):
```unknown
-----
anndata             0.8.0
numpy               1.21.5
pandas              1.3.5
scanpy              1.8.2
scipy               1.7.3
session_info        1.0.0
-----
```

Example 4 (python):
```python
PIL                 9.1.1
asciitree           NA
backcall            0.2.0
beta_ufunc          NA
binom_ufunc         NA
cached_property     1.5.1
cairo               1.20.0
cffi                1.14.4
cloudpickle         1.6.0
colorama            0.4.4
cycler              0.10.0
cython_runtime      NA
cytoolz             0.11.2
dask                2022.01.0
dateutil            2.8.1
decorator           4.4.2
defusedxml          0.6.0
fasteners           0.17.3
fsspec              2022.01.0
google              NA
h5py                3.1.0
igraph              0.10.2
importlib_metadata  1.7.0
ipykernel           5.4.2
ipython_genutils    0.2.0
jedi                0.17.2
jinja2              3.1.2
joblib              1.0.0
kiwisolver          1.3.1
leidenalg           0.9.0
llvmlite            0.35.0
louvain             0.7.0
markupsafe          2.1.1
matplotlib          3.5.2
mpl_toolkits        NA
msgpack             1.0.3
natsort             7.1.0
nbinom_ufunc        NA
numba               0.52.0
numcodecs           0.9.1
numexpr             2.7.2
packaging           21.3
parso               0.7.1
pexpect             4.8.0
pickleshare         0.7.5
pkg_resources       NA
plotly              5.10.0
prompt_toolkit      3.0.10
psutil              5.9.0
ptyprocess          0.7.0
pyarrow             9.0.0
pygments            2.13.0
pyparsing           2.4.7
pytz                2020.5
setuptools_scm      NA
sinfo               0.3.1
six                 1.14.0
sklearn             1.0.2
sphinxcontrib       NA
storemagic          NA
tables              3.6.1
tblib               1.7.0
texttable           1.6.3
threadpoolctl       2.1.0
tlz                 0.11.2
toolz               0.11.2
tornado             6.1
traitlets           5.0.5
typing_extensions   NA
wcwidth             0.2.5
yaml                5.4.1
zarr                2.11.3
zipp                NA
zmq                 21.0.1
zope                NA
```

Example 5 (unknown):
```unknown
-----
IPython             7.19.0
jupyter_client      6.1.11
jupyter_core        4.7.0
notebook            6.2.0
-----
Python 3.7.9 | packaged by conda-forge | (default, Dec  9 2020, 21:08:20) [GCC 9.3.0]
Linux-4.15.0-158-generic-x86_64-with-debian-bullseye-sid
-----
Session information updated at 2023-03-26 09:38
```

Example 6 (python):
```python
sc.settings.set_figure_params(dpi=300)
```

Example 7 (unknown):
```unknown
!ls /nfs/team205/heart/anndata_objects/Foetal/Immune
```

Example 8 (unknown):
```unknown
archive
C97_index.csv
MultiVI-model_foetal_Mono-MP_withContCov-CCCov_13Mar2023_feature-select-0003
MultiVI-model_foetal_Mono-MP_withContCov-CCCov_13Mar2023_feature-select-001
MultiVI-post_foetal_Mono-MP_withContCov-CCCov_13Mar2023_feature-select-0003.h5ad
MultiVI-post_foetal_Mono-MP_withContCov-CCCov_13Mar2023_feature-select-001.h5ad
MultiVI-pre_paired_foetal_MoMP_raw.h5ad
MultiVI-pre_scRNA_foetal_MoMP_raw.h5ad
RNA_foetal_raw_immune_16Feb2023.h5ad
RNA_foetal_raw_immune_16Feb2023_Multiome.h5ad
RNA_foetal_raw_immune_16Feb2023_Multiome.rds
RNA_foetal_raw_immune_16Feb2023.obs.csv
RNA_foetal_raw_Lymphoid_16Feb2023_CCregressed-scaled.h5ad
RNA_foetal_raw_Lymphoid_16Feb2023.h5ad
RNA_foetal_raw_Lymphoid_16Feb2023.obs.csv
RNA_foetal_raw_Myeloid_16Feb2023.h5ad
RNA_foetal_raw_Myeloid_16Feb2023.obs.csv
Tcycling
```

Example 9 (python):
```python
adata = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/Immune/MultiVI-post_foetal_Mono-MP_withContCov-CCCov_13Mar2023_feature-select-0003.h5ad')
adata
```

Example 10 (unknown):
```unknown
AnnData object with n_obs × n_vars = 7674 × 142283
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score', '_indices', '_scvi_batch', '_scvi_labels', 'fine_grain_withcycling', 't', 'seg', 'edge', 't_sd'
    var: 'modality', 'gene_name_scRNA-0', 'gene_id', 'n_cells'
    uns: '_scvi_manager_uuid', '_scvi_uuid', 'draw_graph', 'fine_grain_colors', 'fine_grain_withcycling_colors', 'graph', 'neighbors', 'ppt', 'pseudotime_list', 'seg_colors', 'umap', 'week_colors'
    obsm: 'MultiVI_latent', 'X_R', 'X_draw_graph_fa', 'X_palantir', 'X_umap', 'X_umap13', '_scvi_extra_categorical_covs', '_scvi_extra_continuous_covs'
    obsp: 'connectivities', 'distances'
```

Example 11 (python):
```python
sc.pl.umap(adata, color=['fine_grain_withcycling','week'], size=20, wspace=0.7, components='all', ncols=3)
```

Example 12 (python):
```python
sc.pl.umap(adata, color=['week'], groups=['5W','7W','9W','13W'], size=20, wspace=0.3, components='all', ncols=3)
```

Example 13 (python):
```python
import scFates as scf
import palantir
```

Example 14 (unknown):
```unknown
findfont: Font family ['Raleway'] not found. Falling back to DejaVu Sans.
findfont: Font family ['Lato'] not found. Falling back to DejaVu Sans.
```

Example 15 (python):
```python
projections = pd.DataFrame(adata.obsm["MultiVI_latent"],index=adata.obs_names)
projections
```

Example 16 (unknown):
```unknown
# fate probability mapping of the cells
dm_res = palantir.utils.run_diffusion_maps(projections,knn=50)
ms_data = palantir.utils.determine_multiscale_space(dm_res,n_eigs=5)
```

Example 17 (unknown):
```unknown
Determing nearest neighbor graph...
```

Example 18 (python):
```python
# Generate embedding from the multiscale diffusion space
adata.obsm["X_palantir"]=ms_data.values
sc.pp.neighbors(adata,n_neighbors=50,use_rep="X_palantir")
```

Example 19 (python):
```python
sc.pl.embedding(adata, basis='X_palantir',color=['fine_grain_withcycling'], components='all',ncols=3,wspace=0.6,size=10)
```

Example 20 (unknown):
```unknown
findfont: Font family ['Bitstream Vera Sans'] not found. Falling back to DejaVu Sans.
findfont: Font family ['Bitstream Vera Sans'] not found. Falling back to DejaVu Sans.
```

Example 21 (python):
```python
# draw ForceAtlas2 embedding using UMAP dim1 and dim2 as initial position
# make obsm with UMAP1 and 3
# since it looks well capturing the transition of Mono -> MoMP -> TRM and YS-MP -> TRM
adata.obsm['X_umap13'] = adata.obsm['X_umap'][:,[0,2]]
sc.tl.draw_graph(adata,init_pos='X_umap13',random_state=0)
```

Example 22 (python):
```python
adata.uns["week_colors"][list(adata.obs['week'].cat.categories).index('13W')]='purple'
sc.pl.draw_graph(adata,color="week",groups=['5W','7W', '13W'],color_map="RdBu_r")
```

Example 23 (python):
```python
sc.pl.draw_graph(adata,color="fine_grain_withcycling",color_map="RdBu_r")
```

Example 24 (python):
```python
scf.tl.tree(adata,method="ppt",Nodes=50,use_rep="palantir",
            device="cpu",seed=5)
```

Example 25 (unknown):
```unknown
inferring a principal tree --> parameters used 
    50 principal points, sigma = 0.1, lambda = 1, metric = euclidean
    fitting:  84%|████████▍ | 42/50 [00:04<00:00,  8.86it/s]
    converged
    finished (0:00:04) --> added 
    .uns['ppt'], dictionnary containing inferred tree.
    .obsm['X_R'] soft assignment of cells to principal points.
    .uns['graph']['B'] adjacency matrix of the principal points.
    .uns['graph']['F'] coordinates of principal points in representation space.
```

Example 26 (python):
```python
scf.pl.graph(adata,tips=False,forks=False,nodes=[2,5,29,11,47,28])
```

Example 27 (python):
```python
adata.uns['graph']['pp_seg']
```

Example 28 (python):
```python
# Selecting a root and computing pseudotime
scf.tl.roots(adata,[29,2],meeting=47)
# generate 100 mappings of pseudotime, to account for cell asssignment uncertainty.
# .obs will be saved the mean of all calculated pseudotimes.
scf.tl.pseudotime(adata,n_jobs=20,n_map=1,seed=40)
```

Example 29 (python):
```python
root selected --> added
    2 is the farthest root.
    .uns['graph']['root'] farthest root selected.
    .uns['graph']['root2'] 2nd root selected.
    .uns['graph']['meeting'] meeting point on the tree.
    .uns['graph']['pp_info'] for each PP, its distance vs root and segment assignment.
    .uns['graph']['pp_seg'] segments network information.
projecting cells onto the principal graph
    finished (0:00:05) --> added
    .obs['edge'] assigned edge.
    .obs['t'] pseudotime value.
    .obs['seg'] segment of the tree assigned.
    .obs['milestones'] milestone assigned.
    .uns['pseudotime_list'] list of cell projection from all mappings.
```

Example 30 (python):
```python
scf.pl.trajectory(adata,arrows=True,arrow_offset=3,scale_path=0.5,
                  # root_milestone=29,milestones=27
                 )
```

Example 31 (python):
```python
sc.pl.draw_graph(adata,color=["seg","milestones"],legend_loc='on data',legend_fontsize=10)
# scf.tl.rename_milestones(adata,["bifA","Ery","DC","Root","BifB","Mono"])
# we change the color of the root milestone for better visualisations
# adata.uns["milestones_colors"][3]="#17bece"
```

Example 32 (python):
```python
adata.uns['graph']['pp_seg']
```

Example 33 (python):
```python
# save
adata.write('/nfs/team205/heart/anndata_objects/Foetal/Immune/MultiVI-post_foetal_Mono-MP_withContCov-CCCov_13Mar2023_feature-select-0003.h5ad')
```

Example 34 (python):
```python
# read in
adata = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/Immune/MultiVI-post_foetal_Mono-MP_withContCov-CCCov_13Mar2023_feature-select-0003.h5ad')
adata
```

Example 35 (unknown):
```unknown
AnnData object with n_obs × n_vars = 7674 × 142283
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score', '_indices', '_scvi_batch', '_scvi_labels', 'fine_grain_withcycling', 't', 'seg', 'edge', 't_sd', 'milestones'
    var: 'modality', 'gene_name_scRNA-0', 'gene_id', 'n_cells'
    uns: '_scvi_manager_uuid', '_scvi_uuid', 'draw_graph', 'fine_grain_colors', 'fine_grain_withcycling_colors', 'graph', 'milestones_colors', 'neighbors', 'ppt', 'pseudotime_list', 'seg_colors', 'umap', 'week_colors'
    obsm: 'MultiVI_latent', 'X_R', 'X_draw_graph_fa', 'X_palantir', 'X_umap', 'X_umap13', '_scvi_extra_categorical_covs', '_scvi_extra_continuous_covs'
    obsp: 'connectivities', 'distances'
```

Example 36 (python):
```python
import scanpy as sc
import pandas as pd
import numpy as np
import anndata
import scipy
import scFates as scf
import matplotlib as mpl
import matplotlib.pyplot as plt
# enables correct plotting of text
mpl.rcParams['pdf.fonttype'] = 42
mpl.rcParams['ps.fonttype'] = 42
```

Example 37 (python):
```python
sc.settings.set_figure_params(dpi=80,dpi_save=300)
```

Example 38 (python):
```python
%matplotlib inline
```

Example 39 (unknown):
```unknown
figdir = '/nfs/team205/heart/fetal_ms_figures/raw_plots_kk'
```

Example 40 (python):
```python
# read in adata
adata = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/Immune/MultiVI-post_foetal_Mono-MP_withContCov-CCCov_13Mar2023_feature-select-0003.h5ad')
adata
```

Example 41 (python):
```python
AnnData object with n_obs × n_vars = 7674 × 142283
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score', '_indices', '_scvi_batch', '_scvi_labels', 'fine_grain_withcycling', 't', 'seg', 'edge', 't_sd', 'milestones'
    var: 'modality', 'gene_name_scRNA-0', 'gene_id', 'n_cells'
    uns: '_scvi_manager_uuid', '_scvi_uuid', 'adata_auc_gene', 'adata_auc_region', 'adata_gex', 'adata_gex_lognorm', 'draw_graph', 'fine_grain_colors', 'fine_grain_withcycling_colors', 'graph', 'milestones_colors', 'neighbors', 'ppt', 'pseudotime_list', 'seg_colors', 'umap', 'week_colors'
    obsm: 'MultiVI_latent', 'X_R', 'X_draw_graph_fa', 'X_palantir', 'X_umap', 'X_umap13', '_scvi_extra_categorical_covs', '_scvi_extra_continuous_covs'
    obsp: 'connectivities', 'distances'
```

Example 42 (python):
```python
scf.pl.trajectory(adata,arrows=True,arrow_offset=3,scale_path=0.5, frameon=False,show=False)
plt.savefig(f"{figdir}/Force-directed-graph_MoMP_pseudotime.pdf", bbox_inches="tight",dpi=300)
```

Example 43 (unknown):
```unknown
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored
  cax = scatter(
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Supplementary_materials/Epigenetic_stability/no3_entropy-score.html

**Contents:**
- Functions¶
- Parameters¶
- Read in RNA Milo nhoods data¶
- Read in ATAC Milo nhoods data¶
- Entropy of the diffusion of each RNA nhood in ATAC space¶
- Model entropy¶
- Assign fine-grained celltype to each nhood¶
- Entropy of the diffusion of each RNA nhood in ATAC space (selected cells)¶
- Save¶
- Plot scores on UMAP¶
- Plot entropy scores per cell type¶
- All mid-grained cell types¶
- Kruskal-Wallis Test for two groups¶
  - Atrial cardiomyocytes¶
  - Ventricular cardiomyocytes¶
  - Myeloid cells¶
- For figure plots¶
  - UMAPs¶
  - Bar plots¶

**Examples:**

Example 1 (python):
```python
import warnings
warnings.simplefilter("ignore", UserWarning)

import scanpy as sc
import pandas as pd
import numpy as np
import anndata
import os
from progressbar import progressbar

import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import seaborn as sns

import matplotlib as mpl
# enables correct plotting of text
mpl.rcParams['pdf.fonttype'] = 42
mpl.rcParams['ps.fonttype'] = 42 

from scipy.stats import pearsonr
```

Example 2 (python):
```python
import milopy.core as milo
import milopy.plot as milopl
import milopy.utils
```

Example 3 (python):
```python
import session_info
session_info.show()
```

Example 4 (python):
```python
-----
anndata             0.9.2
matplotlib          3.7.2
milopy              0.1.1
numpy               1.24.4
pandas              1.5.3
progressbar         4.3.2
scanpy              1.9.3
scipy               1.10.1
seaborn             0.12.2
session_info        1.0.0
-----
```

Example 5 (python):
```python
PIL                         10.0.0
anyio                       NA
arrow                       1.2.3
asttokens                   NA
attr                        23.1.0
attrs                       23.1.0
babel                       2.12.1
backcall                    0.2.0
backports                   NA
brotli                      1.0.9
certifi                     2023.11.17
cffi                        1.15.1
charset_normalizer          3.2.0
colorama                    0.4.6
comm                        0.1.4
cycler                      0.10.0
cython_runtime              NA
dateutil                    2.8.2
debugpy                     1.6.8
decorator                   5.1.1
defusedxml                  0.7.1
executing                   1.2.0
fastjsonschema              NA
fqdn                        NA
h5py                        3.9.0
idna                        3.4
igraph                      0.10.6
importlib_metadata          NA
importlib_resources         NA
ipykernel                   6.25.1
isoduration                 NA
jedi                        0.19.0
jinja2                      3.1.2
joblib                      1.3.2
json5                       NA
jsonpointer                 2.0
jsonschema                  4.19.0
jsonschema_specifications   NA
jupyter_events              0.7.0
jupyter_server              2.7.1
jupyterlab_server           2.24.0
kiwisolver                  1.4.5
leidenalg                   0.10.1
llvmlite                    0.40.1
markupsafe                  2.1.3
mpl_toolkits                NA
natsort                     8.4.0
nbformat                    5.9.2
numba                       0.57.1
overrides                   NA
packaging                   23.1
parso                       0.8.3
patsy                       0.5.3
pexpect                     4.8.0
pickleshare                 0.7.5
pkg_resources               NA
plotly                      5.16.1
prometheus_client           NA
prompt_toolkit              3.0.39
psutil                      5.9.5
ptyprocess                  0.7.0
pure_eval                   0.2.2
pydev_ipython               NA
pydevconsole                NA
pydevd                      2.9.5
pydevd_file_utils           NA
pydevd_plugins              NA
pydevd_tracing              NA
pygments                    2.16.1
pyparsing                   3.0.9
python_utils                NA
pythonjsonlogger            NA
pytz                        2023.3
referencing                 NA
requests                    2.31.0
rfc3339_validator           0.1.4
rfc3986_validator           0.1.1
rpds                        NA
rpy2                        3.5.14
ruamel                      NA
send2trash                  NA
six                         1.16.0
sklearn                     1.3.0
sniffio                     1.3.0
socks                       1.7.1
sphinxcontrib               NA
stack_data                  0.6.2
statsmodels                 0.14.0
texttable                   1.6.7
threadpoolctl               3.2.0
tornado                     6.3.3
traitlets                   5.9.0
typing_extensions           NA
tzlocal                     NA
uri_template                NA
urllib3                     2.0.4
wcwidth                     0.2.6
webcolors                   1.13
websocket                   1.6.2
yaml                        6.0.1
zipp                        NA
zmq                         25.1.1
zstandard                   0.19.0
```

Example 6 (unknown):
```unknown
-----
IPython             8.12.0
jupyter_client      8.3.0
jupyter_core        4.12.0
jupyterlab          4.0.5
notebook            7.0.2
-----
Python 3.8.17 (default, Jul  5 2023, 21:04:15) [GCC 11.2.0]
Linux-4.15.0-220-generic-x86_64-with-glibc2.17
-----
Session information updated at 2024-02-15 14:08
```

Example 7 (python):
```python
import importlib.util
import sys
spec = importlib.util.spec_from_file_location("module.name", "/nfs/team205/kk18/function/python/utils.py")
utils = importlib.util.module_from_spec(spec)
sys.modules["module.name"] = utils
spec.loader.exec_module(utils)
```

Example 8 (python):
```python
sc.settings.set_figure_params(dpi=80)
```

Example 9 (python):
```python
from scipy.stats import entropy
import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline

def calculate_entropy(df1, df2, relative_to_max=True):
    """
    Calculate the entropy-based diffusion scores for each group defined in df1 and apply it in df2.

    Parameters:
    - df1: DataFrame representing cell-by-cluster matrix for defining groups.
    - df2: DataFrame representing cell-by-cluster matrix for calculating entropy.

    Returns:
    A dictionary with group names as keys and their entropy scores as values.
    """
    scores = {}
    for group_name in progressbar(df1.columns):
        # Identify cells belonging to the current group
        group_indices = df1.index[df1[group_name] == 1].tolist()
        # Subset df2 for the current group
        group_df2 = df2.loc[group_indices]
        # Calculate the distribution of cells across clusters in df2
        cluster_distribution = group_df2.sum(axis=0) / group_df2.sum(axis=0).sum()
        # Calculate entropy for the group
        if relative_to_max:
            # Prepare maximum distribution
            # Create a DataFrame filled with zeros
            max_distr_df = pd.DataFrame(0, index=group_df2.index, columns=group_df2.columns)
            # Replace elements with 1
            col_indices = [i % max_distr_df.shape[1] for i in range(len(group_indices))]
            max_distr_df.values[range(len(group_indices)), col_indices] = 1
            # display(max_distr_df)
            # Calculate disctribution
            max_distribution = max_distr_df.sum(axis=0) / max_distr_df.sum(axis=0).sum()
            # print(max_distribution)
            scores[group_name] = entropy(pk=cluster_distribution,
                                         qk=max_distribution,
                                        )
        else:
            scores[group_name] = entropy(pk=cluster_distribution)

    return scores



def calculate_kl_divergence(df1, df2):
    """
    Calculate the Kullback-Leibler divergence for each group defined in df1 and apply it in df2.

    Parameters:
    - df1: DataFrame representing cell-by-cluster matrix for defining groups.
    - df2: DataFrame representing cell-by-cluster matrix for calculating KL divergence.

    Returns:
    A dictionary with group names as keys and their KL divergence as values.
    """
    scores = {}
    num_clusters = df2.shape[1]
    uniform_distribution = np.ones(num_clusters) / num_clusters

    for group_name in df1.columns:
        # Identify cells belonging to the current group
        group_indices = df1.index[df1[group_name] == 1].tolist()
        # Subset df2 for the current group
        group_df2 = df2.loc[group_indices]
        # Calculate the distribution of cells across clusters in df2
        cluster_distribution = group_df2.sum(axis=0) / group_df2.sum(axis=0).sum()
        # Calculate KL divergence for the group
        scores[group_name] = entropy(cluster_distribution, uniform_distribution)

    return scores

def model_entropy(y, x:list, degree=2):
    """
    Fit a trend to the nhood-wise entropy against the data size.

    Parameters:
    y (array-like): Numeric vector of dependent variable 
    x (list): list of numeric vector of independent variables (regressor)
    degree (int, optional): The degree of polynomials to be used in the loess fit.

    Returns:
    dict: A dictionary containing the fitted values and the model.
    """

    # Preparing the data for fitting
    X_list = []
    for xx in x:
        X_list.append(np.array(xx))
    Y = np.array(y)
    
    # Combine X1 and X2 into a single feature matrix
    if len(X_list)>1:
        X = np.column_stack(X_list)
    else:
        X = X_list[0]

    # Creating a polynomial regression model
    model = make_pipeline(PolynomialFeatures(degree), LinearRegression())
    model.fit(X, Y)

    # Predicting the fitted values
    fitted_values = model.predict(X)

    return {"fitted_values": fitted_values, "model": model}
```

Example 10 (python):
```python
# dir_to_adata_peakvi = '/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/PeakMatrix/mid_grain_mod'
path_to_gex_milo_nhoods = '/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/Milo_nhoods_gex/scVI_mid_grain_mod'
path_to_atac_milo_nhoods = '/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/Milo_nhoods/mid_grain_mod'
path_to_adata_gex = '/nfs/team205/heart/anndata_objects/Foetal/Feb28ObjectRaw_Multiome.h5ad'
celltype_col = 'mid_grain_mod'
```

Example 11 (unknown):
```unknown
figdir = '/nfs/team205/heart/fetal_ms_figures/raw_plots_kk'
```

Example 12 (unknown):
```unknown
# order of mid grain (mod) cell types
celltypes = [
    'AtrialCardiomyocytes_plusCCS',
    'VentricularCardiomyocytes_plusCCS',
    'Fibroblasts',
    'MuralCells',
    'PericardialCells', 
    'EpicardialCells',
    'BloodVesselEndothelialCells',
    'EndocardialCells', 
     'Neurons',
    'Glia', 
     'MyeloidCells',
    'LymphoidCells',   
             ]
```

Example 13 (python):
```python
# nhood results
prefix = 'one-layer_added-cov-region_subsampled-n1000_'
suffix = '_milo_nhoods_result.csv'

rna_nhoods_res = {}
for file in os.listdir(path_to_gex_milo_nhoods):
    if (prefix in file)&(suffix in file):
        celltype = file.replace(prefix,'').replace(suffix,'')
        print(celltype)
        rna_nhoods_res[celltype] = pd.read_csv(f'{path_to_gex_milo_nhoods}/{prefix}{celltype}{suffix}',
                                           index_col=0
                                          )
        # rename nhood names
        rna_nhoods_res[celltype].index = celltype+'_'+rna_nhoods_res[celltype].index.astype('str')
rna_nhoods_res[celltype].head(n=3)
```

Example 14 (unknown):
```unknown
LymphoidCells
LymphaticEndothelialCells
EndocardialCells
AtrialCardiomyocytes_plusCCS
MuralCells
Neurons
VentricularCardiomyocytes_plusCCS
Fibroblasts
EpicardialCells
MyeloidCells
BloodVesselEndothelialCells
Glia
PericardialCells
```

Example 15 (python):
```python
prefix = 'one-layer_added-cov-region_subsampled-n1000_'
suffix = '_milo_cell-by-nhood.csv'

rna_cell_by_nhood_dfs = {}
for file in os.listdir(path_to_gex_milo_nhoods):
    if (prefix in file)&(suffix in file):
        celltype = file.replace(prefix,'').replace(suffix,'')
        print(celltype)
        rna_cell_by_nhood_dfs[celltype] = pd.read_csv(f'{path_to_gex_milo_nhoods}/{prefix}{celltype}{suffix}',
                                           index_col=0
                                          )
        # rename nhood names
        rna_cell_by_nhood_dfs[celltype].columns = celltype+'_'+rna_cell_by_nhood_dfs[celltype].columns.astype('str')
rna_cell_by_nhood_dfs[celltype].head(n=3)
```

Example 16 (unknown):
```unknown
VentricularCardiomyocytes_plusCCS
AtrialCardiomyocytes_plusCCS
MyeloidCells
PericardialCells
LymphoidCells
BloodVesselEndothelialCells
LymphaticEndothelialCells
Neurons
MuralCells
EpicardialCells
Glia
EndocardialCells
Fibroblasts
```

Example 17 (python):
```python
midfix = '_added-cov-region_subsampled-n1000'
suffix = '_milo_cell-by-nhood.csv'

atac_cell_by_nhood_dfs = {}
for file in os.listdir(path_to_atac_milo_nhoods):
    if (midfix in file)&(suffix in file):
        celltype = file.replace(midfix,'').replace(suffix,'')
        print(celltype)
        atac_cell_by_nhood_dfs[celltype] = pd.read_csv(f'{path_to_atac_milo_nhoods}/{celltype}{midfix}{suffix}',
                                           index_col=0
                                          )
        # rename nhood names
        atac_cell_by_nhood_dfs[celltype].columns = 'ATAC_'+celltype+'_'+atac_cell_by_nhood_dfs[celltype].columns.astype('str')

atac_cell_by_nhood_dfs[celltype].head(n=3)
```

Example 18 (unknown):
```unknown
PericardialCells
Glia
BloodVesselEndothelialCells
AtrialCardiomyocytes_plusCCS
MuralCells
VentricularCardiomyocytes_plusCCS
Neurons
Fibroblasts
EpicardialCells
MyeloidCells
LymphoidCells
EndocardialCells
LymphaticEndothelialCells
```

Example 19 (python):
```python
# for each celltype
for i,celltype in enumerate(celltypes):
    # entropy score
    res = calculate_entropy(df1 = rna_cell_by_nhood_dfs[celltype], 
                            df2 = atac_cell_by_nhood_dfs[celltype],
                            relative_to_max=False
                                                             )
    # add to nhood result dataframe
    rna_nhoods_res[celltype].loc[res.keys(),'entropy score'] = list(res.values())
    # checkign whether there's nan
    if rna_nhoods_res[celltype]['entropy score'].isna().sum()>0:
        print(f'{celltype} has nan in the result')

rna_nhoods_res[celltype].head(n=3)
```

Example 20 (unknown):
```unknown
100% (195 of 195) |######################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (300 of 300) |######################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (300 of 300) |######################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (118 of 118) |######################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (83 of 83) |########################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (65 of 65) |########################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (138 of 138) |######################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (65 of 65) |########################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (50 of 50) |########################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (80 of 80) |########################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (161 of 161) |######################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (284 of 284) |######################| Elapsed Time: 0:00:00 Time:  0:00:00
```

Example 21 (python):
```python
for celltype in celltypes:
    nhood_sizes = np.array(rna_nhoods_res[celltype]['Nhood_size'])
    variances = np.array(rna_nhoods_res[celltype]['variance_scVI-latent'])
    entropy_scores = np.array(rna_nhoods_res[celltype]['entropy score'])
    fitted_values = model_entropy(y=entropy_scores, x=[nhood_sizes,variances], degree=1)['fitted_values']
    residuals = entropy_scores - fitted_values
    # add to the dataframe
    rna_nhoods_res[celltype]['entropy score, fitted'] = fitted_values.copy()
    rna_nhoods_res[celltype]['entropy score, residual'] = residuals.copy()
```

Example 22 (python):
```python
# Plotting original value with fitted line
# prepare figure
nrows = 2
ncols = 6
fig,axs=plt.subplots(nrows=nrows,ncols=ncols,figsize=(22,7)) 
for i,celltype in enumerate(celltypes):
    df = rna_nhoods_res[celltype]
    corr_res = pearsonr(df['Nhood_size'],df['entropy score'])
    axs[i//ncols,i%ncols].scatter(x=df['Nhood_size'],y=df['entropy score'])
    # Sort the sizes array and reorder the fitted values accordingly
    sorted_indices = np.argsort(df['Nhood_size'])
    sorted_sizes = np.array(df['Nhood_size'])[sorted_indices]
    sorted_fitted_values = df['entropy score, fitted'][sorted_indices]
    axs[i//ncols,i%ncols].plot(sorted_sizes,sorted_fitted_values,color='red')
    axs[i//ncols,i%ncols].set_title(f'{celltype}\nPearson corr coeff: {round(corr_res[0],2)}, pval: {round(corr_res[1],3)}')
    axs[i//ncols,i%ncols].set_xlabel('Nhood_size')
    axs[i//ncols,i%ncols].set_ylabel('entropy score')
fig.suptitle('before regressing out')
fig.tight_layout()

fig,axs=plt.subplots(nrows=nrows,ncols=ncols,figsize=(22,7)) 
for i,celltype in enumerate(celltypes):
    df = rna_nhoods_res[celltype]
    corr_res = pearsonr(df['variance_scVI-latent'],df['entropy score'])
    axs[i//ncols,i%ncols].scatter(x=df['variance_scVI-latent'],y=df['entropy score'])
    # Sort the sizes array and reorder the fitted values accordingly
    sorted_indices = np.argsort(df['variance_scVI-latent'])
    sorted_sizes = np.array(df['variance_scVI-latent'])[sorted_indices]
    sorted_fitted_values = df['entropy score, fitted'][sorted_indices]
    axs[i//ncols,i%ncols].plot(sorted_sizes,sorted_fitted_values,color='red')
    axs[i//ncols,i%ncols].set_title(f'{celltype}\nPearson corr coeff: {round(corr_res[0],2)}, pval: {round(corr_res[1],3)}')
    axs[i//ncols,i%ncols].set_xlabel('variance_scVI-latent')
    axs[i//ncols,i%ncols].set_ylabel('entropy score')
fig.suptitle('before regressing out')
fig.tight_layout()
```

Example 23 (python):
```python
# Plot after correction

# prepare figure
nrows = 2
ncols = 6
fig,axs=plt.subplots(nrows=nrows,ncols=ncols,figsize=(22,7)) 
for i,celltype in enumerate(celltypes):
    df = rna_nhoods_res[celltype]
    corr_res = pearsonr(df['Nhood_size'],df['entropy score, residual'])
    axs[i//ncols,i%ncols].scatter(x=df['Nhood_size'],y=df['entropy score, residual'])
    axs[i//ncols,i%ncols].set_title(f'{celltype}\nPearson corr coeff: {round(corr_res[0],2)}, pval: {round(corr_res[1],4)}')
    axs[i//ncols,i%ncols].set_xlabel('Nhood_size')
    axs[i//ncols,i%ncols].set_ylabel('entropy score, residual')
fig.suptitle('regressed out')
fig.tight_layout()

fig,axs=plt.subplots(nrows=nrows,ncols=ncols,figsize=(22,7)) 
for i,celltype in enumerate(celltypes):
    df = rna_nhoods_res[celltype]
    corr_res = pearsonr(df['variance_scVI-latent'],df['entropy score, residual'])
    axs[i//ncols,i%ncols].scatter(x=df['variance_scVI-latent'],y=df['entropy score, residual'])
    axs[i//ncols,i%ncols].set_title(f'{celltype}\nPearson corr coeff: {round(corr_res[0],2)}, pval: {round(corr_res[1],4)}')
    axs[i//ncols,i%ncols].set_xlabel('variance_scVI-latent')
    axs[i//ncols,i%ncols].set_ylabel('entropy score, residual')
fig.suptitle('regressed out')
fig.tight_layout()
```

Example 24 (python):
```python
purity_thresh = 0.6

# read in celltype info
ad = sc.read_h5ad(path_to_adata_gex,backed='r')
obs = ad.obs.copy()

# prepare figure
nrows = 2
ncols = 6
fig,axs=plt.subplots(nrows=nrows,ncols=ncols,figsize=(22,5)) 

rna_cell_by_nhood_dfs_sel = {}
# assign
for i,celltype in enumerate(celltypes):
    df = rna_cell_by_nhood_dfs[celltype].copy()
    # get cell number of each cell type in each nhood
    df.index = obs['fine_grain'].astype('str').reindex(df.index)
    df = df.reset_index().groupby('fine_grain').sum()
    # purity (proportion of the largest cell type)
    purity_arr = df.max(0)/df.sum(0)
    # get pure nhood
    pure_nhoods = list(purity_arr.index[purity_arr>purity_thresh])
    axs[i//ncols,i%ncols].hist(purity_arr,bins=50)
    axs[i//ncols,i%ncols].axvline(x=purity_thresh, color='red', linestyle='--', linewidth=1.5)
    axs[i//ncols,i%ncols].set_title(f'{celltype}\nprop of pure nhoods: {round(len(pure_nhoods)/len(purity_arr),3)}')
    # assign celltype
    df_pure = df[pure_nhoods]
    mask = df_pure == df_pure.max(0)
    assigned_pre = mask.reset_index().melt(id_vars=['fine_grain']).rename(columns={'value':'assigned',
                                                                               'variable':'nhood',
                                                                               'fine_grain':'assigned_fine_grain'
                                                                              })
    assigned_pre['purity_fine_grain'] = list(purity_arr[assigned_pre['nhood']])
    assigned = assigned_pre[assigned_pre['assigned']]
    # get cells of the assigned_celltype
    rna_cell_by_nhood_dfs_sel[celltype] = pd.DataFrame(0,
                                                       columns=rna_cell_by_nhood_dfs[celltype].columns,
                                                       index=rna_cell_by_nhood_dfs[celltype].index)
    for ind in assigned.index:
        celltype_assigned = assigned.loc[ind,'assigned_fine_grain']
        nh = assigned.loc[ind,'nhood']
        nh_cells = list(rna_cell_by_nhood_dfs[celltype].index[rna_cell_by_nhood_dfs[celltype][nh]==1])
        celltype_cells = list(obs.index[obs['fine_grain']==celltype_assigned])
        cells = list(set(nh_cells).intersection(celltype_cells))
        rna_cell_by_nhood_dfs_sel[celltype].loc[cells,nh] = 1
        rna_nhoods_res[celltype].loc[nh,'cells_of_assigned_celltype'] = ','.join(cells)
        rna_nhoods_res[celltype].loc[nh,'Nhood_size_sel'] = len(cells)
    
    # add to nhood result dataframe
    rna_nhoods_res[celltype].loc[assigned['nhood'],'assigned_fine_grain'] = list(assigned['assigned_fine_grain'])
    rna_nhoods_res[celltype].loc[assigned['nhood'],'purity_fine_grain'] = list(assigned['purity_fine_grain'])
    # sanity check
    if sum(rna_nhoods_res[celltype]['assigned_fine_grain'].isna()==False)/len(rna_nhoods_res[celltype])!=len(pure_nhoods)/len(purity_arr):
        print(f'{celltype}: something went wrong when adding to result dataframe')

fig.tight_layout()
rna_nhoods_res[celltype].head(n=3)
```

Example 25 (unknown):
```unknown
# just checking the calculations are consistent
sum(rna_nhoods_res[celltype]['Nhood_size_sel']/rna_nhoods_res[celltype]['Nhood_size'] == rna_nhoods_res[celltype]['purity_fine_grain'])
```

Example 26 (python):
```python
# for each celltype
for i,celltype in enumerate(celltypes):
    # entropy score
    res = calculate_entropy(df1 = rna_cell_by_nhood_dfs_sel[celltype], 
                            df2 = atac_cell_by_nhood_dfs[celltype],
                            relative_to_max=False
                                                             )
    # add to nhood result dataframe
    rna_nhoods_res[celltype].loc[res.keys(),'entropy score, selected cells'] = list(res.values())
    # checkign whether there's more nan than expected
    if rna_nhoods_res[celltype]['entropy score, selected cells'].isna().sum()>rna_nhoods_res[celltype]['assigned_fine_grain'].isna().sum():
        print(f'{celltype} has larger nan than expected')

rna_nhoods_res[celltype].head(n=3)
```

Example 27 (unknown):
```unknown
100% (195 of 195) |######################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (300 of 300) |######################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (300 of 300) |######################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (118 of 118) |######################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (83 of 83) |########################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (65 of 65) |########################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (138 of 138) |######################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (65 of 65) |########################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (50 of 50) |########################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (80 of 80) |########################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (161 of 161) |######################| Elapsed Time: 0:00:00 Time:  0:00:00
100% (284 of 284) |######################| Elapsed Time: 0:00:00 Time:  0:00:00
```

Example 28 (python):
```python
# regress out the effect of Nhood_size
for celltype in celltypes:
    df_sel = rna_nhoods_res[celltype].dropna(subset='assigned_fine_grain')
    nhood_sizes = np.array(df_sel['Nhood_size'])
    variances = np.array(df_sel['variance_scVI-latent'])
    entropy_scores = np.array(df_sel['entropy score'])
    fitted_values = model_entropy(y=entropy_scores, x=[nhood_sizes,variances], degree=1)['fitted_values']
    residuals = entropy_scores - fitted_values
    df_out = pd.DataFrame({'nhood':df_sel.index,'fitted':fitted_values,'residual':residuals})
    df_out.set_index('nhood',inplace=True)
    # add to the dataframe
    rna_nhoods_res[celltype]['entropy score, selected cells, fitted'] = df_out['fitted'].reindex(rna_nhoods_res[celltype].index).copy()
    rna_nhoods_res[celltype]['entropy score, selected cells, residual'] = df_out['residual'].reindex(rna_nhoods_res[celltype].index).copy()
```

Example 29 (unknown):
```unknown
utils.save_pkl(rna_nhoods_res,
               '/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/Milo_nhoods_gex/scVI_mid_grain_mod/rna_nhoods_res.pkl'
              )
```

Example 30 (unknown):
```unknown
# read in
rna_nhoods_res = utils.read_pkl('/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/Milo_nhoods_gex/scVI_mid_grain_mod/rna_nhoods_res.pkl')
rna_nhoods_res.keys()
```

Example 31 (unknown):
```unknown
dict_keys(['LymphoidCells', 'LymphaticEndothelialCells', 'EndocardialCells', 'AtrialCardiomyocytes_plusCCS', 'MuralCells', 'Neurons', 'VentricularCardiomyocytes_plusCCS', 'Fibroblasts', 'EpicardialCells', 'MyeloidCells', 'BloodVesselEndothelialCells', 'Glia', 'PericardialCells'])
```

Example 32 (python):
```python
# read in anndata include UMAP embedding
gex_adata_dict = utils.read_pkl(f'{path_to_gex_milo_nhoods}/one-layer_added-cov-region_subsampled-n1000_post_milo_adata_dictionary.pkl')
# atac_adata_dict = utils.read_pkl(f'{path_to_atac_milo_nhoods}/added-cov-region_subsampled-n1000_post_milo_adata_dictionary.pkl')
```

Example 33 (python):
```python
?sc.pl.embedding
```

Example 34 (python):
```python
Signature:
sc.pl.embedding(
    adata: anndata._core.anndata.AnnData,
    basis: str,
    *,
    color: Union[str, Sequence[str], NoneType] = None,
    gene_symbols: Union[str, NoneType] = None,
    use_raw: Union[bool, NoneType] = None,
    sort_order: bool = True,
    edges: bool = False,
    edges_width: float = 0.1,
    edges_color: Union[str, Sequence[float], Sequence[str]] = 'grey',
    neighbors_key: Union[str, NoneType] = None,
    arrows: bool = False,
    arrows_kwds: Union[Mapping[str, Any], NoneType] = None,
    groups: Union[str, NoneType] = None,
    components: Union[str, Sequence[str]] = None,
    dimensions: Union[Tuple[int, int], Sequence[Tuple[int, int]], NoneType] = None,
    layer: Union[str, NoneType] = None,
    projection: Literal['2d', '3d'] = '2d',
    scale_factor: Union[float, NoneType] = None,
    color_map: Union[matplotlib.colors.Colormap, str, NoneType] = None,
    cmap: Union[matplotlib.colors.Colormap, str, NoneType] = None,
    palette: Union[str, Sequence[str], cycler.Cycler, NoneType] = None,
    na_color: Union[str, Tuple[float, ...]] = 'lightgray',
    na_in_legend: bool = True,
    size: Union[float, Sequence[float], NoneType] = None,
    frameon: Union[bool, NoneType] = None,
    legend_fontsize: Union[int, float, Literal['xx-small', 'x-small', 'small', 'medium', 'large', 'x-large', 'xx-large'], NoneType] = None,
    legend_fontweight: Union[int, Literal['light', 'normal', 'medium', 'semibold', 'bold', 'heavy', 'black']] = 'bold',
    legend_loc: str = 'right margin',
    legend_fontoutline: Union[int, NoneType] = None,
    colorbar_loc: Union[str, NoneType] = 'right',
    vmax: Union[str, float, Callable[[Sequence[float]], float], Sequence[Union[str, float, Callable[[Sequence[float]], float]]], NoneType] = None,
    vmin: Union[str, float, Callable[[Sequence[float]], float], Sequence[Union[str, float, Callable[[Sequence[float]], float]]], NoneType] = None,
    vcenter: Union[str, float, Callable[[Sequence[float]], float], Sequence[Union[str, float, Callable[[Sequence[float]], float]]], NoneType] = None,
    norm: Union[matplotlib.colors.Normalize, Sequence[matplotlib.colors.Normalize], NoneType] = None,
    add_outline: Union[bool, NoneType] = False,
    outline_width: Tuple[float, float] = (0.3, 0.05),
    outline_color: Tuple[str, str] = ('black', 'white'),
    ncols: int = 4,
    hspace: float = 0.25,
    wspace: Union[float, NoneType] = None,
    title: Union[str, Sequence[str], NoneType] = None,
    show: Union[bool, NoneType] = None,
    save: Union[bool, str, NoneType] = None,
    ax: Union[matplotlib.axes._axes.Axes, NoneType] = None,
    return_fig: Union[bool, NoneType] = None,
    **kwargs,
) -> Union[matplotlib.figure.Figure, matplotlib.axes._axes.Axes, NoneType]
Docstring:
Scatter plot for user specified embedding basis (e.g. umap, pca, etc)

Parameters
----------
basis : str
    Name of the `obsm` basis to use.
adata : AnnData
    Annotated data matrix.
color : typing.Union[str, typing.Sequence[str], NoneType], optional (default: None)
    Keys for annotations of observations/cells or variables/genes, e.g.,
    `'ann1'` or `['ann1', 'ann2']`.
gene_symbols : typing.Union[str, NoneType], optional (default: None)
    Column name in `.var` DataFrame that stores gene symbols. By default `var_names`
    refer to the index column of the `.var` DataFrame. Setting this option allows
    alternative names to be used.
use_raw : typing.Union[bool, NoneType], optional (default: None)
    Use `.raw` attribute of `adata` for coloring with gene expression. If `None`,
    defaults to `True` if `layer` isn't provided and `adata.raw` is present.
layer : typing.Union[str, NoneType], optional (default: None)
    Name of the AnnData object layer that wants to be plotted. By default
    adata.raw.X is plotted. If `use_raw=False` is set, then `adata.X` is plotted.
    If `layer` is set to a valid layer name, then the layer is plotted. `layer`
    takes precedence over `use_raw`.
edges : bool, optional (default: False)
    Show edges.
edges_width : float, optional (default: 0.1)
    Width of edges.
edges_color : typing.Union[str, typing.Sequence[float], typing.Sequence[str]], optional (default: 'grey')
    Color of edges. See :func:`~networkx.drawing.nx_pylab.draw_networkx_edges`.
neighbors_key : typing.Union[str, NoneType], optional (default: None)
    Where to look for neighbors connectivities.
    If not specified, this looks .obsp['connectivities'] for connectivities
    (default storage place for pp.neighbors).
    If specified, this looks
    .obsp[.uns[neighbors_key]['connectivities_key']] for connectivities.
arrows : bool, optional (default: False)
    Show arrows (deprecated in favour of `scvelo.pl.velocity_embedding`).
arrows_kwds : typing.Union[typing.Mapping[str, typing.Any], NoneType], optional (default: None)
    Passed to :meth:`~matplotlib.axes.Axes.quiver`
sort_order : bool, optional (default: True)
    For continuous annotations used as color parameter, plot data points
    with higher values on top of others.
groups : typing.Union[str, NoneType], optional (default: None)
    Restrict to a few categories in categorical observation annotation.
    The default is not to restrict to any groups.
dimensions : typing.Union[typing.Tuple[int, int], typing.Sequence[typing.Tuple[int, int]], NoneType], optional (default: None)
    0-indexed dimensions of the embedding to plot as integers. E.g. [(0, 1), (1, 2)].
    Unlike `components`, this argument is used in the same way as `colors`, e.g. is
    used to specify a single plot at a time. Will eventually replace the components
    argument.
components : typing.Union[str, typing.Sequence[str]], optional (default: None)
    For instance, `['1,2', '2,3']`. To plot all available components use
    `components='all'`.
projection : typing.Literal['2d', '3d'], optional (default: '2d')
    Projection of plot (default: `'2d'`).
legend_loc : str, optional (default: 'right margin')
    Location of legend, either `'on data'`, `'right margin'` or a valid keyword
    for the `loc` parameter of :class:`~matplotlib.legend.Legend`.
legend_fontsize : typing.Union[int, float, typing.Literal['xx-small', 'x-small', 'small', 'medium', 'large', 'x-large', 'xx-large'], NoneType], optional (default: None)
    Numeric size in pt or string describing the size.
    See :meth:`~matplotlib.text.Text.set_fontsize`.
legend_fontweight : typing.Union[int, typing.Literal['light', 'normal', 'medium', 'semibold', 'bold', 'heavy', 'black']], optional (default: 'bold')
    Legend font weight. A numeric value in range 0-1000 or a string.
    Defaults to `'bold'` if `legend_loc == 'on data'`, otherwise to `'normal'`.
    See :meth:`~matplotlib.text.Text.set_fontweight`.
legend_fontoutline : typing.Union[int, NoneType], optional (default: None)
    Line width of the legend font outline in pt. Draws a white outline using
    the path effect :class:`~matplotlib.patheffects.withStroke`.
colorbar_loc : typing.Union[str, NoneType], optional (default: 'right')
    Where to place the colorbar for continous variables. If `None`, no colorbar
    is added.
size : typing.Union[float, typing.Sequence[float], NoneType], optional (default: None)
    Point size. If `None`, is automatically computed as 120000 / n_cells.
    Can be a sequence containing the size for each cell. The order should be
    the same as in adata.obs.
color_map : typing.Union[matplotlib.colors.Colormap, str, NoneType], optional (default: None)
    Color map to use for continous variables. Can be a name or a
    :class:`~matplotlib.colors.Colormap` instance (e.g. `"magma`", `"viridis"`
    or `mpl.cm.cividis`), see :func:`~matplotlib.cm.get_cmap`.
    If `None`, the value of `mpl.rcParams["image.cmap"]` is used.
    The default `color_map` can be set using :func:`~scanpy.set_figure_params`.
palette : typing.Union[str, typing.Sequence[str], cycler.Cycler, NoneType], optional (default: None)
    Colors to use for plotting categorical annotation groups.
    The palette can be a valid :class:`~matplotlib.colors.ListedColormap` name
    (`'Set2'`, `'tab20'`, …), a :class:`~cycler.Cycler` object, a dict mapping
    categories to colors, or a sequence of colors. Colors must be valid to
    matplotlib. (see :func:`~matplotlib.colors.is_color_like`).
    If `None`, `mpl.rcParams["axes.prop_cycle"]` is used unless the categorical
    variable already has colors stored in `adata.uns["{var}_colors"]`.
    If provided, values of `adata.uns["{var}_colors"]` will be set.
na_color : typing.Union[str, typing.Tuple[float, ...]], optional (default: 'lightgray')
    Color to use for null or masked values. Can be anything matplotlib accepts as a
    color. Used for all points if `color=None`.
na_in_legend : bool, optional (default: True)
    If there are missing values, whether they get an entry in the legend. Currently
    only implemented for categorical legends.
frameon : typing.Union[bool, NoneType], optional (default: None)
    Draw a frame around the scatter plot. Defaults to value set in
    :func:`~scanpy.set_figure_params`, defaults to `True`.
title : typing.Union[str, typing.Sequence[str], NoneType], optional (default: None)
    Provide title for panels either as string or list of strings,
    e.g. `['title1', 'title2', ...]`.

vmin : typing.Union[str, float, typing.Callable[[typing.Sequence[float]], float], typing.Sequence[typing.Union[str, float, typing.Callable[[typing.Sequence[float]], float]]], NoneType], optional (default: None)
    The value representing the lower limit of the color scale. Values smaller than vmin are plotted
    with the same color as vmin. vmin can be a number, a string, a function or `None`. If
    vmin is a string and has the format `pN`, this is interpreted as a vmin=percentile(N).
    For example vmin='p1.5' is interpreted as the 1.5 percentile. If vmin is function, then
    vmin is interpreted as the return value of the function over the list of values to plot.
    For example to set vmin tp the mean of the values to plot, `def my_vmin(values): return
    np.mean(values)` and then set `vmin=my_vmin`. If vmin is None (default) an automatic
    minimum value is used as defined by matplotlib `scatter` function. When making multiple
    plots, vmin can be a list of values, one for each plot. For example `vmin=[0.1, 'p1', None, my_vmin]`
vmax : typing.Union[str, float, typing.Callable[[typing.Sequence[float]], float], typing.Sequence[typing.Union[str, float, typing.Callable[[typing.Sequence[float]], float]]], NoneType], optional (default: None)
    The value representing the upper limit of the color scale. The format is the same as for `vmin`.
vcenter : typing.Union[str, float, typing.Callable[[typing.Sequence[float]], float], typing.Sequence[typing.Union[str, float, typing.Callable[[typing.Sequence[float]], float]]], NoneType], optional (default: None)
    The value representing the center of the color scale. Useful for diverging colormaps.
    The format is the same as for `vmin`.
    Example: sc.pl.umap(adata, color='TREM2', vcenter='p50', cmap='RdBu_r')
add_outline : typing.Union[bool, NoneType], optional (default: False)
    If set to True, this will add a thin border around groups of dots. In some situations
    this can enhance the aesthetics of the resulting image
outline_color : typing.Tuple[str, str], optional (default: ('black', 'white'))
    Tuple with two valid color names used to adjust the add_outline. The first color is the
    border color (default: black), while the second color is a gap color between the
    border color and the scatter dot (default: white).
outline_width : typing.Tuple[float, float], optional (default: (0.3, 0.05))
    Tuple with two width numbers used to adjust the outline. The first value is the width
    of the border color as a fraction of the scatter dot size (default: 0.3). The second value is
    width of the gap color (default: 0.05).
ncols : int, optional (default: 4)
    Number of panels per row.
wspace : typing.Union[float, NoneType], optional (default: None)
    Adjust the width of the space between multiple panels.
hspace : float, optional (default: 0.25)
    Adjust the height of the space between multiple panels.
return_fig : typing.Union[bool, NoneType], optional (default: None)
    Return the matplotlib figure.
kwargs : _empty
    Arguments to pass to :func:`matplotlib.pyplot.scatter`,
    for instance: the maximum and minimum values (e.g. `vmin=-2, vmax=5`).
show : typing.Union[bool, NoneType], optional (default: None)
     Show the plot, do not return axis.
save : typing.Union[bool, str, NoneType], optional (default: None)
    If `True` or a `str`, save the figure.
    A string is appended to the default filename.
    Infer the filetype if ending on {`'.pdf'`, `'.png'`, `'.svg'`}.
ax : typing.Union[matplotlib.axes._axes.Axes, NoneType], optional (default: None)
    A matplotlib axes object. Only works if plotting a single component.

Returns
-------
If `show==False` a :class:`~matplotlib.axes.Axes` or a list of it.
File:      /nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/scanpy/plotting/_tools/scatterplots.py
Type:      function
```

Example 35 (unknown):
```unknown
len(celltypes)
```

Example 36 (unknown):
```unknown
i=0, row:0, col:0,1
i=1, row:0, col:2,3
i=2, row:1, col:0,1
i=3, row:1, col:2,3
i=4, row:2, col:0,1
i=5, row:2, col:2,3
```

Example 37 (python):
```python
i = 3
if i%2==0:
    col_pos = i%4
else:
    col_pos = (i%4)+1
print(col_pos)
print(col_pos+1)
```

Example 38 (python):
```python
# prepare dictionary to store labels and colors
colors = {}
fig,axs=plt.subplots(nrows=6,ncols=4,figsize=(20,30))
for i,celltype in enumerate(celltypes):
    ad = gex_adata_dict[celltype].uns['nhood_adata']
    ad.obs['nhood'] = [f'{celltype}_nhood_{str(int(x)+1)}' for x in ad.obs.index]
    # add values
    ad.obs[f'assigned_fine_grain'] = list(rna_nhoods_res[celltype]['assigned_fine_grain'].reindex(ad.obs['nhood']))
    ad.obs[f'entropy score'] = list(rna_nhoods_res[celltype]['entropy score, residual'].reindex(ad.obs['nhood']))
    # print(celltype)
    row_pos = i//2
    if i%2==0:
        col_pos_1 = 0
        col_pos_2 = 1
    else:
        col_pos_1 = 2
        col_pos_2 = 3
    sc.pl.embedding(ad,basis='X_milo_graph',size=ad.obs['Nhood_size']*7,
                    color=['entropy score'],
                    cmap='inferno', frameon=False, wspace=0.15, 
                    title=f'{celltype}\nentropy score', show=False,
                    ax=axs[row_pos,col_pos_1]
                   )
    sc.pl.embedding(ad,basis='X_milo_graph',size=ad.obs['Nhood_size']*7,
                    color=['assigned_fine_grain'],
                    cmap='inferno', frameon=False, wspace=0.15, 
                    title=f'{celltype}\nassigned_fine_grain', show=False,
                    ax=axs[row_pos,col_pos_2]
                   )
    # store colors
    colors[celltype]={}
    colors[celltype]['labels'] = list(ad.obs[f'assigned_fine_grain'].cat.categories)
    colors[celltype]['colors'] = ad.uns[f'assigned_fine_grain_colors']
fig.savefig(f"{figdir}/entropyscore_fine-grain-assignment.pdf", bbox_inches="tight",dpi=300)
```

Example 39 (python):
```python
nrows = 3
ncols = 4
fig,axs=plt.subplots(nrows=nrows,ncols=ncols,figsize=(30,12))
for i,celltype in enumerate(celltypes):
    ad = gex_adata_dict[celltype].uns['nhood_adata']
    ad.obs['nhood'] = [f'{celltype}_nhood_{str(int(x)+1)}' for x in ad.obs.index]
    # add fine-grained labels
    ad.obs[celltype] = list(rna_nhoods_res[celltype]['assigned_fine_grain'].reindex(ad.obs['nhood']))
    sc.pl.embedding(ad,basis='X_milo_graph',size=ad.obs['Nhood_size']*7,
                    color=[celltype], # legend_loc='on data', 
                    legend_fontsize=12,
                    cmap='inferno', frameon=False,
                    show=False, ax=axs[i//ncols,i%ncols]
                   )
fig.tight_layout(w_pad=0.2)
```

Example 40 (python):
```python
nrows = 3
ncols = 4
fig,axs=plt.subplots(nrows=nrows,ncols=ncols,figsize=(20,12))
for i,celltype in enumerate(celltypes):
    ad = gex_adata_dict[celltype].uns['nhood_adata']
    ad.obs['nhood'] = [f'{celltype}_nhood_{str(int(x)+1)}' for x in ad.obs.index]
    # add entropy scores
    ad.obs[celltype] = list(rna_nhoods_res[celltype]['entropy score, residual'].reindex(ad.obs['nhood']))
    sc.pl.embedding(ad,basis='X_milo_graph',size=ad.obs['Nhood_size']*7,
                    color=[celltype],
                    cmap='inferno', frameon=False,
                    show=False, ax=axs[i//ncols,i%ncols]
                   )
fig.tight_layout()
```

Example 41 (python):
```python
# Prepare concatenated dataframe
for i,celltype in enumerate(celltypes):
    subset = rna_nhoods_res[celltype].copy()
    subset['mid_grain_mod'] = celltype
    if i==0:
        nhoods_res_concat = subset.copy()
    else:
        nhoods_res_concat = pd.concat([nhoods_res_concat,subset])
nhoods_res_concat['mid_grain_mod']=nhoods_res_concat['mid_grain_mod'].astype('category')
nhoods_res_concat['mid_grain_mod'].cat.reorder_categories(celltypes,inplace=True)
nhoods_res_concat=nhoods_res_concat.dropna(subset='assigned_fine_grain')

# Create a catplot with boxplots for each gene
g = sns.catplot(x='assigned_fine_grain', y='entropy score, selected cells, residual', col='mid_grain_mod', data=nhoods_res_concat, kind='box',
                col_wrap=3, sharey=False, sharex=False, height=6, aspect=0.8, palette="Blues")

# Overlay swarmplot on each Axes
# We'll use the 'hue' parameter to match the subsets and turn off dodge to align the points with the boxes.
for celltype, ax in zip(celltypes, g.axes.flat):
    subset = nhoods_res_concat[nhoods_res_concat['mid_grain_mod'] == celltype]
    sns.swarmplot(x='assigned_fine_grain', y='entropy score, selected cells, residual', data=subset, color='k', size=5, ax=ax, dodge=False)
    
# Update titles to remove " ="
for ax, celltype in zip(g.axes.flat, celltypes):
    ax.set_title(celltype)

# Rotate x-tick labels and remove x-axis labels
for i,ax in enumerate(g.axes.flat):
    ax.set_xticklabels(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor")  # Rotate x-tick labels
    # ax.set_xticklabels('')
    ax.set_xlabel('')  # Remove x-axis labels
    if i!=0:
        ax.set_ylabel('')

# Adjust the spacing of the plots
g.fig.subplots_adjust(wspace=0.2, hspace=1.2)

# Show the plot
plt.show()
```

Example 42 (unknown):
```unknown
/tmp/ipykernel_60954/2274242884.py:10: FutureWarning: The `inplace` parameter in pandas.Categorical.reorder_categories is deprecated and will be removed in a future version. Reordering categories will always return a new Categorical object.
  nhoods_res_concat['mid_grain_mod'].cat.reorder_categories(celltypes,inplace=True)
```

Example 43 (python):
```python
# Prepare concatenated dataframe
for i,celltype in enumerate(celltypes):
    subset = rna_nhoods_res[celltype].copy()
    subset['mid_grain_mod'] = celltype
    if i==0:
        nhoods_res_concat = subset.copy()
    else:
        nhoods_res_concat = pd.concat([nhoods_res_concat,subset])
nhoods_res_concat['mid_grain_mod']=nhoods_res_concat['mid_grain_mod'].astype('category')
nhoods_res_concat['mid_grain_mod'].cat.reorder_categories(celltypes,inplace=True)
nhoods_res_concat=nhoods_res_concat.dropna(subset='assigned_fine_grain')

# Create a catplot with boxplots for each gene
g = sns.catplot(x='assigned_fine_grain', y='purity_fine_grain', col='mid_grain_mod', data=nhoods_res_concat, kind='box',
                col_wrap=6, sharey=False, sharex=False, height=6, aspect=0.8, palette="Blues")

# Overlay swarmplot on each Axes
# We'll use the 'hue' parameter to match the subsets and turn off dodge to align the points with the boxes.
for celltype, ax in zip(celltypes, g.axes.flat):
    subset = nhoods_res_concat[nhoods_res_concat['mid_grain_mod'] == celltype]
    sns.swarmplot(x='assigned_fine_grain', y='purity_fine_grain', data=subset, color='k', size=5, ax=ax, dodge=False)
    
# Update titles to remove " ="
for ax, celltype in zip(g.axes.flat, celltypes):
    ax.set_title(celltype)

# Rotate x-tick labels and remove x-axis labels
for i,ax in enumerate(g.axes.flat):
    ax.set_xticklabels(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor")  # Rotate x-tick labels
    # ax.set_xticklabels('')
    ax.set_xlabel('')  # Remove x-axis labels
    if i!=0:
        ax.set_ylabel('')

# Adjust the spacing of the plots
g.fig.subplots_adjust(wspace=0.2, hspace=1)

# Show the plot
plt.show()
```

Example 44 (unknown):
```unknown
/tmp/ipykernel_60954/3286563489.py:10: FutureWarning: The `inplace` parameter in pandas.Categorical.reorder_categories is deprecated and will be removed in a future version. Reordering categories will always return a new Categorical object.
  nhoods_res_concat['mid_grain_mod'].cat.reorder_categories(celltypes,inplace=True)
```

Example 45 (python):
```python
import scipy.stats as stats
import scikit_posthocs as spost
```

Example 46 (unknown):
```unknown
celltype = 'AtrialCardiomyocytes_plusCCS'
assigned = rna_nhoods_res[celltype].copy()
assigned = assigned.dropna(subset='assigned_fine_grain')
set(assigned['assigned_fine_grain'])
```

Example 47 (unknown):
```unknown
{'AtrialCardiomyocytesCycling',
 'AtrialCardiomyocytesLeft',
 'AtrialCardiomyocytesRight',
 'AtrioventricularNodeCardiomyocytes',
 'SinoatrialNodeCardiomyocytes'}
```

Example 48 (python):
```python
celltype = 'AtrialCardiomyocytes_plusCCS'
map_celltypes = {
    'AtrialCardiomyocytesCycling':'AtrialWorkingCM',
     'AtrialCardiomyocytesLeft':'AtrialWorkingCM',
     'AtrialCardiomyocytesRight':'AtrialWorkingCM',
     'AtrioventricularNodeCardiomyocytes':'AtrialCCS',
     'SinoatrialNodeCardiomyocytes':'AtrialCCS'
}

assigned = rna_nhoods_res[celltype].copy()
assigned = assigned.dropna(subset='assigned_fine_grain')
assigned['mapped_celltype'] = assigned['assigned_fine_grain'].copy()
assigned.replace({'mapped_celltype':map_celltypes},inplace=True)

# Grouping data by SampleLabel
grouped_data = assigned.groupby('mapped_celltype')['entropy score, selected cells, residual'].apply(list)

# Kruskal-Wallis Test
stat, p_value = stats.kruskal(*grouped_data)
print("Kruskal-Wallis Test:")
print(f"Statistic: {stat}, P-value: {p_value}")
print("")

# also confirm with mann-whitney-u test (scipy)
stat, p_value = stats.mannwhitneyu(*grouped_data)
print("Mann-Whitney U rank Test:")
print(f"Statistic: {stat}, P-value: {p_value}")
```

Example 49 (unknown):
```unknown
Kruskal-Wallis Test:
Statistic: 10.054975332593926, P-value: 0.0015193705181662676

Mann-Whitney U rank Test:
Statistic: 1791.0, P-value: 0.0015293858581563156
```

Example 50 (unknown):
```unknown
celltype = 'VentricularCardiomyocytes_plusCCS'
assigned = rna_nhoods_res[celltype].copy()
assigned = assigned.dropna(subset='assigned_fine_grain')
set(assigned['assigned_fine_grain'])
```

Example 51 (unknown):
```unknown
{'VentricularCardiomyocytesCycling',
 'VentricularCardiomyocytesLeftCompact',
 'VentricularCardiomyocytesLeftTrabeculated',
 'VentricularCardiomyocytesRightCompact',
 'VentricularCardiomyocytesRightTrabeculated',
 'VentricularConductionSystemDistal',
 'VentricularConductionSystemProximal'}
```

Example 52 (python):
```python
celltype = 'VentricularCardiomyocytes_plusCCS'
map_celltypes = {
    'VentricularCardiomyocytesCycling':'VentricularWorkingCM',
     'VentricularCardiomyocytesLeftCompact':'VentricularWorkingCM',
     'VentricularCardiomyocytesLeftTrabeculated':'VentricularWorkingCM',
     'VentricularCardiomyocytesRightCompact':'VentricularWorkingCM',
     'VentricularCardiomyocytesRightTrabeculated':'VentricularWorkingCM',
     'VentricularConductionSystemDistal':'VentricularCCS',
     'VentricularConductionSystemProximal':'VentricularCCS'
}

assigned = rna_nhoods_res[celltype].copy()
assigned = assigned.dropna(subset='assigned_fine_grain')
assigned['mapped_celltype'] = assigned['assigned_fine_grain'].copy()
assigned.replace({'mapped_celltype':map_celltypes},inplace=True)

# Grouping data by SampleLabel
grouped_data = assigned.groupby('mapped_celltype')['entropy score, selected cells, residual'].apply(list)

# Kruskal-Wallis Test
stat, p_value = stats.kruskal(*grouped_data)
print("Kruskal-Wallis Test:")
print(f"Statistic: {stat}, P-value: {p_value}")
print("")

# also confirm with mann-whitney-u test (scipy)
stat, p_value = stats.mannwhitneyu(*grouped_data)
print("Mann-Whitney U rank Test:")
print(f"Statistic: {stat}, P-value: {p_value}")
```

Example 53 (unknown):
```unknown
Kruskal-Wallis Test:
Statistic: 39.17342310868321, P-value: 3.877830158632655e-10

Mann-Whitney U rank Test:
Statistic: 974.0, P-value: 3.922386003787531e-10
```

Example 54 (unknown):
```unknown
celltype = 'MyeloidCells'
assigned = rna_nhoods_res[celltype].copy()
assigned = assigned.dropna(subset='assigned_fine_grain')
set(assigned['assigned_fine_grain'])
```

Example 55 (unknown):
```unknown
{'DendriticCellsMature',
 'DendriticCellsType1',
 'MacrophagesATF3pos',
 'MacrophagesCX3CR1pos',
 'MacrophagesLYVE1pos',
 'MacrophagesTIMD4pos',
 'MastCells',
 'Megakaryocytes',
 'MonocyteDerivedCells',
 'Monocytes',
 'PlasmacytoidDendriticCells'}
```

Example 56 (python):
```python
celltype = 'MyeloidCells'
map_celltypes = {
    'Monocytes':'Mo,MP',
    'MonocyteDerivedCells':'Mo,MP',
    'MacrophagesCX3CR1pos':'Mo,MP',
    'MacrophagesLYVE1pos':'Mo,MP',
    'MacrophagesTIMD4pos':'Mo,MP',
    'MacrophagesATF3pos':'Mo,MP',
    'DendriticCellsType1':'DC',
    'DendriticCellsMature':'DC',
    'PlasmacytoidDendriticCells':'DC',
     'MastCells':'Mast',
     'Megakaryocytes':'Mega',
}

assigned = rna_nhoods_res[celltype].copy()
assigned = assigned.dropna(subset='assigned_fine_grain')
assigned['mapped_celltype'] = assigned['assigned_fine_grain'].copy()
assigned.replace({'mapped_celltype':map_celltypes},inplace=True)

# Grouping data by SampleLabel
grouped_data = assigned.groupby('mapped_celltype')['entropy score, selected cells, residual'].apply(list)

# Kruskal-Wallis Test
stat, p_value = stats.kruskal(*grouped_data)
print("Kruskal-Wallis Test:")
print(f"Statistic: {stat}, P-value: {p_value}")

# Perform Post Hoc analysis if Kruskal-Wallis test is significant
if p_value < 0.05:
    print("\nPerforming Post Hoc Analysis with Dunn's test...")

    # Post Hoc analysis using scikit-posthocs
    posthoc_results = spost.posthoc_dunn(assigned, 
                                         val_col='entropy score, selected cells, residual', 
                                         group_col='mapped_celltype', 
                                         p_adjust='bonferroni')

    # Print results
    print("Post Hoc Analysis (Dunn's test with Bonferroni correction):")
    display(pd.DataFrame(posthoc_results))
else:
    print("No significant differences found. No need for Post Hoc analysis.")
```

Example 57 (unknown):
```unknown
Kruskal-Wallis Test:
Statistic: 89.22165969298436, P-value: 3.219039117656211e-19

Performing Post Hoc Analysis with Dunn's test...
Post Hoc Analysis (Dunn's test with Bonferroni correction):
```

Example 58 (unknown):
```unknown
['AtrialCardiomyocytes_plusCCS',
 'VentricularCardiomyocytes_plusCCS',
 'Fibroblasts',
 'MuralCells',
 'PericardialCells',
 'EpicardialCells',
 'BloodVesselEndothelialCells',
 'EndocardialCells',
 'Neurons',
 'Glia',
 'MyeloidCells',
 'LymphoidCells']
```

Example 59 (python):
```python
cells_to_plot = ['AtrialCardiomyocytes_plusCCS','VentricularCardiomyocytes_plusCCS','MyeloidCells']

nrows = 1
ncols = 3
fig,axs=plt.subplots(nrows=nrows,ncols=ncols,figsize=(17,5))
for i,celltype in enumerate(cells_to_plot):
    ad = gex_adata_dict[celltype].uns['nhood_adata']
    ad.obs['nhood'] = [f'{celltype}_nhood_{str(int(x)+1)}' for x in ad.obs.index]
    # add entropy scores
    ad.obs[celltype] = list(rna_nhoods_res[celltype]['entropy score, residual'].reindex(ad.obs['nhood']))
    sc.pl.embedding(ad,basis='X_milo_graph',size=ad.obs['Nhood_size']*7,
                    color=[celltype],
                    cmap='inferno', frameon=False,
                    show=False, ax=axs[i]
                   )
fig.tight_layout()
fig.savefig(f"{figdir}/entropyscore_UMAPs_for-main.pdf",bbox_inches="tight",dpi=300)
```

Example 60 (python):
```python
celltype = 'AtrialCardiomyocytes_plusCCS'
relabel = {
    'AtrialCardiomyocytesRight':'Right',
    'AtrialCardiomyocytesLeft':'Left',
    'AtrialCardiomyocytesCycling':'Cycling',
    'SinoatrialNodeCardiomyocytes':'SAN_P',
     'AtrioventricularNodeCardiomyocytes':'AVN_P',
}
map_celltypes = {
    'AtrialCardiomyocytesRight':'WorkingCM',
    'AtrialCardiomyocytesLeft':'WorkingCM',
    'AtrialCardiomyocytesCycling':'WorkingCM',
    'SinoatrialNodeCardiomyocytes':'CCS',
     'AtrioventricularNodeCardiomyocytes':'CCS',
}
map_colors = {
    'WorkingCM':'lightsalmon',
     'CCS':'lightskyblue',
}

assigned = rna_nhoods_res[celltype].copy()
assigned = assigned.dropna(subset='assigned_fine_grain')
assigned['mapped_celltype'] = assigned['assigned_fine_grain'].copy()
assigned.replace({'mapped_celltype':map_celltypes},inplace=True)
assigned['relabeled'] = assigned['assigned_fine_grain'].copy()
assigned.replace({'relabeled':relabel},inplace=True)
assigned['relabeled'] = pd.Categorical(assigned['relabeled'],
                                              categories=relabel.values(),
                                              ordered=True)

# Create a catplot with boxplots
g = sns.catplot(x='relabeled', y='entropy score, selected cells, residual', 
                col=None, data=assigned, kind='box',
                col_wrap=None, sharey=False, sharex=False, height=4, aspect=0.9, 
                palette=[map_colors[x] for x in map_celltypes.values()],
               )
sns.swarmplot(x='relabeled', y='entropy score, selected cells, residual', 
              data=assigned, color='k', size=5,dodge=False)

# Create a list of patches for the legend
patches = [mpatches.Patch(color=color, label=label, linewidth=3) for label, color in map_colors.items()]
# Add legend to the plot
plt.legend(handles=patches, frameon=False)
plt.title('Atrial Cardiomyocytes', pad=50)
plt.xticks(rotation=45, ha="right", rotation_mode="anchor")  # Rotate x-tick labels
plt.xlabel('') # Remove x-axis labels
plt.ylabel('Entropy score')
plt.savefig(f"{figdir}/entropyscore_barplot_{celltype}.pdf", bbox_inches="tight",dpi=300)
plt.show()
```

Example 61 (python):
```python
celltype = 'VentricularCardiomyocytes_plusCCS'
relabel = {
    'VentricularCardiomyocytesRightTrabeculated':'RightTrab',
    'VentricularCardiomyocytesRightCompact':'RightComp',
    'VentricularCardiomyocytesLeftTrabeculated':'LeftTrab',
    'VentricularCardiomyocytesLeftCompact':'LeftComp',
    'VentricularCardiomyocytesCycling':'Cycling',
   'VentricularConductionSystemProximal':'VCSproximal',
     'VentricularConductionSystemDistal':'VCSdistal',
}
map_celltypes = {
    'VentricularCardiomyocytesRightTrabeculated':'WorkingCM',
    'VentricularCardiomyocytesRightCompact':'WorkingCM',
    'VentricularCardiomyocytesLeftTrabeculated':'WorkingCM',
    'VentricularCardiomyocytesLeftCompact':'WorkingCM',
    'VentricularCardiomyocytesCycling':'WorkingCM',
   'VentricularConductionSystemProximal':'CCS',
     'VentricularConductionSystemDistal':'CCS',
}
map_colors = {
    'WorkingCM':'lightsalmon',
     'CCS':'lightskyblue',
}

assigned = rna_nhoods_res[celltype].copy()
assigned = assigned.dropna(subset='assigned_fine_grain')
assigned['mapped_celltype'] = assigned['assigned_fine_grain'].copy()
assigned.replace({'mapped_celltype':map_celltypes},inplace=True)
assigned['relabeled'] = assigned['assigned_fine_grain'].copy()
assigned.replace({'relabeled':relabel},inplace=True)
assigned['relabeled'] = pd.Categorical(assigned['relabeled'],
                                              categories=relabel.values(),
                                              ordered=True)

# Create a catplot with boxplots
g = sns.catplot(x='relabeled', y='entropy score, selected cells, residual', 
                col=None, data=assigned, kind='box',
                col_wrap=None, sharey=False, sharex=False, height=4, aspect=1.1, 
                palette=[map_colors[x] for x in map_celltypes.values()],
               )
sns.swarmplot(x='relabeled', y='entropy score, selected cells, residual', 
              data=assigned, color='k', size=5,dodge=False)

# Create a list of patches for the legend
patches = [mpatches.Patch(color=color, label=label, linewidth=3) for label, color in map_colors.items()]
# Add legend to the plot
plt.legend(handles=patches, frameon=False)
plt.title('Ventricular Cardiomyocytes', pad=50)
plt.xticks(rotation=45, ha="right", rotation_mode="anchor")  # Rotate x-tick labels
plt.xlabel('') # Remove x-axis labels
plt.ylabel('Entropy score')
plt.savefig(f"{figdir}/entropyscore_barplot_{celltype}.pdf", bbox_inches="tight",dpi=300)
plt.show()
```

Example 62 (python):
```python
celltype = 'MyeloidCells'
relabel = {
    'Monocytes':'Mono',
    'MonocyteDerivedCells':'MoDerived',
    'MacrophagesCX3CR1pos':'MP_CX3CR1',
    'MacrophagesLYVE1pos':'MP_LYVE1',
    'MacrophagesTIMD4pos':'MP_TIMD4',
    'MacrophagesATF3pos':'MP_ATF3',
    'DendriticCellsType1':'DC1',
    'DendriticCellsMature':'DCmature',
    'PlasmacytoidDendriticCells':'pDC',
     'MastCells':'Mast',
     'Megakaryocytes':'Mega',
}
map_celltypes = {
    'Monocytes':'Mo,MP',
    'MonocyteDerivedCells':'Mo,MP',
    'MacrophagesCX3CR1pos':'Mo,MP',
    'MacrophagesLYVE1pos':'Mo,MP',
    'MacrophagesTIMD4pos':'Mo,MP',
    'MacrophagesATF3pos':'Mo,MP',
    'DendriticCellsType1':'DC',
    'DendriticCellsMature':'DC',
    'PlasmacytoidDendriticCells':'DC',
     'MastCells':'Mast',
     'Megakaryocytes':'Mega',
}
map_colors = {
    'Mo,MP':'lightsalmon',
     'DC':'lightskyblue',
    'Mast':'limegreen',
    'Mega':'plum',
}

assigned = rna_nhoods_res[celltype].copy()
assigned = assigned.dropna(subset='assigned_fine_grain')
assigned['mapped_celltype'] = assigned['assigned_fine_grain'].copy()
assigned.replace({'mapped_celltype':map_celltypes},inplace=True)
assigned['relabeled'] = assigned['assigned_fine_grain'].copy()
assigned.replace({'relabeled':relabel},inplace=True)
assigned['relabeled'] = pd.Categorical(assigned['relabeled'],
                                              categories=relabel.values(),
                                              ordered=True)

# Create a catplot with boxplots
g = sns.catplot(x='relabeled', y='entropy score, selected cells, residual', 
                col=None, data=assigned, kind='box',
                col_wrap=None, sharey=False, sharex=False, height=4, aspect=1.4, 
                palette=[map_colors[x] for x in map_celltypes.values()],
               )
sns.swarmplot(x='relabeled', y='entropy score, selected cells, residual', 
              data=assigned, color='k', size=5,dodge=False)

# Create a list of patches for the legend
patches = [mpatches.Patch(color=color, label=label, linewidth=3) for label, color in map_colors.items()]
# Add legend to the plot
plt.legend(handles=patches, frameon=False, bbox_to_anchor=(1.25, 1))
plt.title('Myeloid cells', pad=90)
plt.xticks(rotation=45, ha="right", rotation_mode="anchor")  # Rotate x-tick labels
plt.xlabel('') # Remove x-axis labels
plt.ylabel('Entropy score')
plt.savefig(f"{figdir}/entropyscore_barplot_{celltype}.pdf", bbox_inches="tight",dpi=300)
plt.show()
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Supplementary_materials/Macrophages/Mo-MP_no5_milestone_pseudotime_plot.html

**Contents:**
- Read in trajectory object¶

**Examples:**

Example 1 (python):
```python
import scanpy as sc
import pandas as pd
import numpy as np
import anndata
import scipy
import os

import matplotlib as mpl
# enables correct plotting of text
mpl.rcParams['pdf.fonttype'] = 42
mpl.rcParams['ps.fonttype'] = 42
```

Example 2 (python):
```python
import scFates as scf
```

Example 3 (python):
```python
import session_info
session_info.show()
```

Example 4 (unknown):
```unknown
-----
anndata             0.9.2
numpy               1.24.4
pandas              1.5.3
scFates             1.0.6
scanpy              1.9.3
scipy               1.10.1
session_info        1.0.0
-----
```

Example 5 (python):
```python
PIL                         10.0.0
adjustText                  0.8
anyio                       NA
arrow                       1.2.3
asttokens                   NA
attr                        23.1.0
attrs                       23.1.0
babel                       2.12.1
backcall                    0.2.0
backports                   NA
brotli                      1.0.9
certifi                     2023.11.17
cffi                        1.15.1
charset_normalizer          3.2.0
colorama                    0.4.6
comm                        0.1.4
cycler                      0.10.0
cython_runtime              NA
dateutil                    2.8.2
debugpy                     1.6.8
decorator                   5.1.1
defusedxml                  0.7.1
executing                   1.2.0
fastjsonschema              NA
fqdn                        NA
h5py                        3.9.0
idna                        3.4
igraph                      0.10.6
importlib_metadata          NA
importlib_resources         NA
ipykernel                   6.25.1
isoduration                 NA
jedi                        0.19.0
jinja2                      3.1.2
joblib                      1.3.2
json5                       NA
jsonpointer                 2.0
jsonschema                  4.19.0
jsonschema_specifications   NA
jupyter_events              0.7.0
jupyter_server              2.7.1
jupyterlab_server           2.24.0
kiwisolver                  1.4.5
leidenalg                   0.10.1
llvmlite                    0.40.1
markupsafe                  2.1.3
matplotlib                  3.7.2
mpl_toolkits                NA
natsort                     8.4.0
nbformat                    5.9.2
networkx                    2.8.8
numba                       0.57.1
overrides                   NA
packaging                   23.1
parso                       0.8.3
patsy                       0.5.3
pexpect                     4.8.0
pickleshare                 0.7.5
pkg_resources               NA
plotly                      5.16.1
prometheus_client           NA
prompt_toolkit              3.0.39
psutil                      5.9.5
ptyprocess                  0.7.0
pure_eval                   0.2.2
pydev_ipython               NA
pydevconsole                NA
pydevd                      2.9.5
pydevd_file_utils           NA
pydevd_plugins              NA
pydevd_tracing              NA
pygments                    2.16.1
pyparsing                   3.0.9
pythonjsonlogger            NA
pytz                        2023.3
referencing                 NA
requests                    2.31.0
rfc3339_validator           0.1.4
rfc3986_validator           0.1.1
rpds                        NA
rpy2                        3.5.14
ruamel                      NA
seaborn                     0.12.2
send2trash                  NA
simpleppt                   1.1.4
six                         1.16.0
sklearn                     1.3.0
skmisc                      0.2.0
sniffio                     1.3.0
socks                       1.7.1
sphinxcontrib               NA
stack_data                  0.6.2
statsmodels                 0.14.0
texttable                   1.6.7
threadpoolctl               3.2.0
tornado                     6.3.3
tqdm                        4.66.1
traitlets                   5.9.0
typing_extensions           NA
tzlocal                     NA
uri_template                NA
urllib3                     2.0.4
wcwidth                     0.2.6
webcolors                   1.13
websocket                   1.6.2
yaml                        6.0.1
zipp                        NA
zmq                         25.1.1
zstandard                   0.19.0
```

Example 6 (unknown):
```unknown
-----
IPython             8.12.0
jupyter_client      8.3.0
jupyter_core        4.12.0
jupyterlab          4.0.5
notebook            7.0.2
-----
Python 3.8.17 (default, Jul  5 2023, 21:04:15) [GCC 11.2.0]
Linux-5.4.0-169-generic-x86_64-with-glibc2.17
-----
Session information updated at 2024-02-06 15:22
```

Example 7 (python):
```python
import importlib.util
import sys
spec = importlib.util.spec_from_file_location("module.name", "/nfs/team205/kk18/function/python/utils.py")
utils = importlib.util.module_from_spec(spec)
sys.modules["module.name"] = utils
spec.loader.exec_module(utils)
```

Example 8 (python):
```python
sc.settings.set_figure_params(dpi=120, dpi_save=300, transparent=True)
```

Example 9 (python):
```python
%matplotlib inline
```

Example 10 (unknown):
```unknown
figdir = '/nfs/team205/heart/fetal_ms_figures/raw_plots_kk'
```

Example 11 (python):
```python
adata = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/Immune/MoMP_Aug2023/MultiVI-post_foetal_Mono-MP_withContCov-CCCov_13Mar2023_feature-select-0003.h5ad')

# rename obs_names
adata.obs_names = [x.replace('_expression','').replace('_paired','') for x in adata.obs_names]
adata
```

Example 12 (python):
```python
AnnData object with n_obs × n_vars = 7674 × 142283
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score', '_indices', '_scvi_batch', '_scvi_labels', 'fine_grain_withcycling', 't', 'seg', 'edge', 't_sd', 'milestones', 'fine_grain_mod'
    var: 'modality', 'gene_name_scRNA-0', 'gene_id', 'n_cells'
    uns: '_scvi_manager_uuid', '_scvi_uuid', 'adata_auc_gene', 'adata_auc_region', 'adata_gex', 'adata_gex_lognorm', 'adata_peak_tfidf', 'draw_graph', 'fine_grain_colors', 'fine_grain_mod_colors', 'fine_grain_withcycling_colors', 'graph', 'milestones_colors', 'neighbors', 'ppt', 'pseudotime_list', 'seg_colors', 'umap', 'week_colors'
    obsm: 'MultiVI_latent', 'X_R', 'X_draw_graph_fa', 'X_palantir', 'X_umap', 'X_umap13', '_scvi_extra_categorical_covs', '_scvi_extra_continuous_covs'
    obsp: 'connectivities', 'distances'
```

Example 13 (python):
```python
sc.pl.draw_graph(adata,color=["fine_grain_mod"],color_map="RdBu_r",wspace=0.7,show=False)
```

Example 14 (unknown):
```unknown
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
```

Example 15 (unknown):
```unknown
<Axes: title={'center': 'fine_grain_mod'}, xlabel='FA1', ylabel='FA2'>
```

Example 16 (python):
```python
scf.pl.trajectory(adata,arrows=True,arrow_offset=2,scale_path=0.5)
```

Example 17 (unknown):
```unknown
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored
  cax = scatter(
```

Example 18 (python):
```python
scf.pl.graph(adata,tips=True,forks=True)
```

Example 19 (unknown):
```unknown
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored
  cax = scatter(
```

Example 20 (python):
```python
sc.pl.draw_graph(adata,color=["seg","milestones"],legend_loc='on data',legend_fontsize=10)
```

Example 21 (unknown):
```unknown
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
```

Example 22 (python):
```python
adata.obs.milestones.value_counts()
```

Example 23 (unknown):
```unknown
47    3001
5     1111
28     875
2      635
11     624
32     522
44     341
20     340
29     207
16      18
Name: milestones, dtype: int64
```

Example 24 (python):
```python
data = adata.uns['graph']['pp_seg']
data
```

Example 25 (unknown):
```unknown
# save dataframe
data.to_csv('./graph_pp_seg.csv',index=None)
```

Example 26 (python):
```python
import networkx as nx
import matplotlib.pyplot as plt
from networkx import *

# node color: age
# node size: population size

# Create a directed graph
G = nx.from_pandas_edgelist(data, 'from', 'to', create_using=nx.DiGraph())

# Creating a dictionary for edge labels
edge_labels = dict(zip(zip(data['from'], data['to']), data['n']))

# Position of the node
coord = pd.DataFrame(adata.obsm['X_draw_graph_fa'])
coord['milestones'] = list(adata.obs['milestones'])
coord = coord.groupby('milestones').mean()
coord_dict = {}
for milestone in coord.index:
    coord_dict[int(milestone)] = np.array(coord.loc[milestone])

# size of the node: cell number
size_dict_pre = adata.obs.milestones.value_counts().to_dict()
sizes = []
for key in list(G):
    sizes.append(size_dict_pre[str(key)])
sizes = np.array(sizes)

# color of the node: psendotime
adata.obs["week_continuous"] = [int(x.split('W')[0]) for x in adata.obs["week"]]
df = adata.obs[['week_continuous','milestones']]
df = df.groupby('milestones').mean()
df['week_scale'] = (df['week_continuous'] - df['week_continuous'].min())/(df['week_continuous'].max() - df['week_continuous'].min())
df.index = df.index.astype('int')
colors = np.array(df.loc[list(G)]['week_scale'])

# Draw the network
plt.figure(figsize=(5, 5))
nx.draw(G, coord_dict, with_labels=True, node_color=colors,
        node_size=sizes*1.2, arrows=True, edgecolors='black',
        cmap='Greens',
       )

# Add edge labels
nx.draw_networkx_edge_labels(G, coord_dict, edge_labels=edge_labels,
                             font_size=10, font_color='red'
                            )
plt.savefig(f"{figdir}/MoMP_milestone_segment_plots.pdf", bbox_inches="tight",dpi=300)
plt.show()
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Supplementary_materials/Macrophages/Mo-MP_no7_MP-LYVE1_fitting-genes-pathway-analysis.html

**Contents:**
- Functions¶
- Variables¶
- Read in anndata object¶
- Read in pseudotime-fitting output¶
- Read in differential test output (vs root)¶
- Select genes¶
- Gene enrichment analysis¶
- Plotting gene transition¶

15922 rows × 8 columns

15922 rows × 9 columns

111454 rows × 7 columns

**Examples:**

Example 1 (python):
```python
import scanpy as sc
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import anndata
import scipy
import os

import matplotlib as mpl
# enables correct plotting of text
mpl.rcParams['pdf.fonttype'] = 42
mpl.rcParams['ps.fonttype'] = 42
```

Example 2 (python):
```python
import session_info
session_info.show()
```

Example 3 (python):
```python
-----
anndata             0.9.2
matplotlib          3.7.2
numpy               1.24.4
pandas              1.5.3
scanpy              1.9.3
scipy               1.10.1
seaborn             0.12.2
session_info        1.0.0
-----
```

Example 4 (python):
```python
PIL                         10.0.0
anyio                       NA
arrow                       1.2.3
asttokens                   NA
attr                        23.1.0
attrs                       23.1.0
babel                       2.12.1
backcall                    0.2.0
brotli                      1.0.9
certifi                     2023.11.17
cffi                        1.15.1
charset_normalizer          3.2.0
colorama                    0.4.6
comm                        0.1.4
cycler                      0.10.0
cython_runtime              NA
dateutil                    2.8.2
debugpy                     1.6.8
decorator                   5.1.1
defusedxml                  0.7.1
executing                   1.2.0
fastjsonschema              NA
fqdn                        NA
h5py                        3.9.0
idna                        3.4
igraph                      0.10.6
importlib_metadata          NA
importlib_resources         NA
ipykernel                   6.25.1
isoduration                 NA
jedi                        0.19.0
jinja2                      3.1.2
joblib                      1.3.2
json5                       NA
jsonpointer                 2.0
jsonschema                  4.19.0
jsonschema_specifications   NA
jupyter_events              0.7.0
jupyter_server              2.7.1
jupyterlab_server           2.24.0
kiwisolver                  1.4.5
leidenalg                   0.10.1
llvmlite                    0.40.1
markupsafe                  2.1.3
mpl_toolkits                NA
natsort                     8.4.0
nbformat                    5.9.2
numba                       0.57.1
overrides                   NA
packaging                   23.1
parso                       0.8.3
patsy                       0.5.3
pexpect                     4.8.0
pickleshare                 0.7.5
pkg_resources               NA
plotly                      5.16.1
prometheus_client           NA
prompt_toolkit              3.0.39
psutil                      5.9.5
ptyprocess                  0.7.0
pure_eval                   0.2.2
pydev_ipython               NA
pydevconsole                NA
pydevd                      2.9.5
pydevd_file_utils           NA
pydevd_plugins              NA
pydevd_tracing              NA
pygments                    2.16.1
pyparsing                   3.0.9
pythonjsonlogger            NA
pytz                        2023.3
referencing                 NA
requests                    2.31.0
rfc3339_validator           0.1.4
rfc3986_validator           0.1.1
rpds                        NA
ruamel                      NA
send2trash                  NA
six                         1.16.0
sklearn                     1.3.0
sniffio                     1.3.0
socks                       1.7.1
sphinxcontrib               NA
stack_data                  0.6.2
statsmodels                 0.14.0
texttable                   1.6.7
threadpoolctl               3.2.0
tornado                     6.3.3
traitlets                   5.9.0
typing_extensions           NA
uri_template                NA
urllib3                     2.0.4
wcwidth                     0.2.6
webcolors                   1.13
websocket                   1.6.2
yaml                        6.0.1
zipp                        NA
zmq                         25.1.1
zstandard                   0.19.0
```

Example 5 (unknown):
```unknown
-----
IPython             8.12.0
jupyter_client      8.3.0
jupyter_core        4.12.0
jupyterlab          4.0.5
notebook            7.0.2
-----
Python 3.8.17 (default, Jul  5 2023, 21:04:15) [GCC 11.2.0]
Linux-5.4.0-169-generic-x86_64-with-glibc2.17
-----
Session information updated at 2024-02-07 14:28
```

Example 6 (python):
```python
import importlib.util
import sys
spec = importlib.util.spec_from_file_location("module.name", "/nfs/team205/kk18/function/python/utils.py")
utils = importlib.util.module_from_spec(spec)
sys.modules["module.name"] = utils
spec.loader.exec_module(utils)
```

Example 7 (python):
```python
sc.settings.set_figure_params(dpi=300,dpi_save=300)
```

Example 8 (unknown):
```unknown
figdir = '/nfs/team205/heart/fetal_ms_figures/raw_plots_kk'
outdir = '/nfs/team205/heart/anndata_objects/Foetal/Immune/MoMP_Aug2023/seg_3-4-5-6-7-8-9'
os.listdir(outdir)
```

Example 9 (python):
```python
['adata.h5ad', 'differential-test_result.pkl', 'fitting_result.pkl']
```

Example 10 (python):
```python
adata = sc.read_h5ad(f'{outdir}/adata.h5ad')
adata
```

Example 11 (python):
```python
AnnData object with n_obs × n_vars = 3886 × 142283
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score', '_indices', '_scvi_batch', '_scvi_labels', 'fine_grain_withcycling', 't', 'seg', 'edge', 't_sd', 'milestones', 'fine_grain_mod'
    var: 'modality', 'gene_name_scRNA-0', 'gene_id', 'n_cells'
    uns: '_scvi_manager_uuid', '_scvi_uuid', 'adata_auc_gene', 'adata_auc_region', 'adata_gex', 'adata_gex_lognorm', 'adata_peak_tfidf', 'draw_graph', 'fine_grain_colors', 'fine_grain_mod_colors', 'fine_grain_withcycling_colors', 'graph', 'milestones_colors', 'neighbors', 'ppt', 'pseudotime_list', 'seg_colors', 'umap', 'week_colors'
    obsm: 'MultiVI_latent', 'X_R', 'X_draw_graph_fa', 'X_palantir', 'X_umap', 'X_umap13', '_scvi_extra_categorical_covs', '_scvi_extra_continuous_covs'
    obsp: 'connectivities', 'distances'
```

Example 12 (python):
```python
# scale t
# min-max scaling
adata.obs['t_scaled'] = (adata.obs['t']-adata.obs['t'].min())/(adata.obs['t'].max()-adata.obs['t'].min())
```

Example 13 (python):
```python
# summarise gene expression
adata_gex = adata.uns['adata_gex_lognorm'].copy()
print(adata_gex.X.data[:10])

# scale t
# min-max scaling
adata_gex.obs['t_scaled'] = (adata_gex.obs['t']-adata_gex.obs['t'].min())/(adata_gex.obs['t'].max()-adata_gex.obs['t'].min())
```

Example 14 (unknown):
```unknown
[1.0841573 2.6877544 1.0841573 2.1779046 1.0841573 2.1779046 1.0841573
 1.0841573 2.9238348 2.3781319]
```

Example 15 (python):
```python
group_mean_exp = utils.grouped_obs_mean(adata_gex, group_key='milestones')
group_mean_exp
```

Example 16 (unknown):
```unknown
# read in
fit_res_dict = utils.read_pkl(f'{outdir}/fitting_result.pkl')
fit_res = fit_res_dict['gene'].copy()
fit_res
```

Example 17 (unknown):
```unknown
# read in
diff_res_dict = utils.read_pkl(f'{outdir}/differential-test_result.pkl')
diff_res = diff_res_dict['gene'].copy()

# add mean gene expression
df = group_mean_exp.reset_index().melt(id_vars=['index'])
df.rename(columns={'index':'names','variable':'group','value':'mean_gene_exp'},
          inplace=True)
# add
diff_res = diff_res.merge(df, how='left', on=['group','names'])

diff_res
```

Example 18 (unknown):
```unknown
diff_res[diff_res['names']=='DOCK4']
```

Example 19 (unknown):
```unknown
set(diff_res['group'])
```

Example 20 (unknown):
```unknown
{'16', '20', '28', '32', '44', '47', '5'}
```

Example 21 (python):
```python
sns.scatterplot(data=diff_res, x='logfoldchanges',y='scores')
```

Example 22 (unknown):
```unknown
<Axes: xlabel='logfoldchanges', ylabel='scores'>
```

Example 23 (python):
```python
sns.scatterplot(data=diff_res,x='logfoldchanges',y='pvals_adj')
```

Example 24 (unknown):
```unknown
<Axes: xlabel='logfoldchanges', ylabel='pvals_adj'>
```

Example 25 (python):
```python
sns.scatterplot(data=diff_res,x='logfoldchanges',y='mean_gene_exp')
```

Example 26 (unknown):
```unknown
<Axes: xlabel='logfoldchanges', ylabel='mean_gene_exp'>
```

Example 27 (python):
```python
diff_logFC_thresh = 0.5
diff_pval_thresh = 0.05
diff_meanexp_thresh = 0.5
fit_pval_thresh = 0.05
fit_amplitude_thresh = 0.5

up_genes = list(set(diff_res[(diff_res['logfoldchanges']>diff_logFC_thresh)&\
                             (diff_res['pvals_adj']<diff_pval_thresh)&\
                             (diff_res['mean_gene_exp']>diff_meanexp_thresh)
                            ]['names']))
print(f'Upregulated gene number: {len(up_genes)}')
fit_genes = list(set(fit_res[(fit_res['fdr']<fit_pval_thresh)&(fit_res['A']>fit_amplitude_thresh)]['gene_name_scRNA-0']))
print(f'Fitted gene number: {len(fit_genes)}')

# take intersection
genes_of_interest = list(set(up_genes).intersection(fit_genes))
print(f'Final gene number: {len(genes_of_interest)}')
```

Example 28 (unknown):
```unknown
Upregulated gene number: 996
Fitted gene number: 2948
Final gene number: 956
```

Example 29 (python):
```python
# ligands
ligands = list(pd.read_table('/nfs/team205/kk18/data/geneset/HGNC/HGNC_542_Receptor_ligands.txt')['Approved symbol'])
print(len(ligands))
genes_of_interest_ligands = list(set(genes_of_interest).intersection(ligands))
print(len(genes_of_interest_ligands))
```

Example 30 (unknown):
```unknown
genes_of_interest_ligands
```

Example 31 (unknown):
```unknown
['LGALS1',
 'CFLAR',
 'TNFSF13B',
 'PSAP',
 'SPP1',
 'ICAM1',
 'RGS2',
 'ANXA1',
 'LGALS9',
 'CCL3',
 'ADM',
 'CCL4',
 'APP',
 'CCL4L2',
 'CXCL8',
 'VEGFA',
 'CCL3L1',
 'IL1B',
 'EDA']
```

Example 32 (python):
```python
import gseapy as gp
```

Example 33 (python):
```python
geneset_list = ['MSigDB_Hallmark_2020','KEGG_2021_Human','GO_Biological_Process_2023','Reactome_2022']
enr = gp.enrichr(gene_list=genes_of_interest, # or "./tests/data/gene_list.txt",
                 gene_sets=geneset_list,
                 organism='human', # don't forget to set organism to the one you desired! e.g. Yeast
                 outdir=None, # don't write to disk
                )

geneset = 'MSigDB_Hallmark_2020'
df = enr.results[enr.results['Gene_set'] == geneset]
df = df[df['Adjusted P-value'] < 0.1]
df = df.sort_values('Combined Score',ascending=False)

# plot
n_rank = 10
plt.rcParams['axes.grid'] = False
plt.rcParams['figure.figsize'] = 4,3
plt.barh(width=df[:n_rank]['Combined Score'],
         y=[x.split(' (')[0] for x in df[:n_rank]['Term']],
         color='darkred')
plt.gca().invert_yaxis()
plt.xlabel('Score')
plt.title(f'All genes\n({geneset})')
plt.margins(y=0.02)
plt.savefig(f"{figdir}/MoMP_LYVE1-up-genes_all-genes_pathways.pdf", bbox_inches="tight",dpi=300)
plt.show()

for geneset in geneset_list:
    print(geneset)
    df = enr.results[enr.results['Gene_set'] == geneset]
    display(df.sort_values('Combined Score',ascending=False)[:30])
```

Example 34 (unknown):
```unknown
MSigDB_Hallmark_2020
```

Example 35 (unknown):
```unknown
KEGG_2021_Human
```

Example 36 (unknown):
```unknown
GO_Biological_Process_2023
```

Example 37 (unknown):
```unknown
Reactome_2022
```

Example 38 (unknown):
```unknown
df = enr.results[enr.results['Gene_set'] == 'MSigDB_Hallmark_2020']
genes = df.set_index('Term').loc['TNF-alpha Signaling via NF-kB','Genes']
genes = genes.split(';')
genes
```

Example 39 (python):
```python
# ligands

geneset_list = ['MSigDB_Hallmark_2020','KEGG_2021_Human','GO_Biological_Process_2023','Reactome_2022']
enr = gp.enrichr(gene_list=genes_of_interest_ligands, # or "./tests/data/gene_list.txt",
                 gene_sets=geneset_list,
                 organism='human', # don't forget to set organism to the one you desired! e.g. Yeast
                 outdir=None, # don't write to disk
                )

geneset = 'MSigDB_Hallmark_2020'
df = enr.results[enr.results['Gene_set'] == geneset]
df = df[df['Adjusted P-value'] < 0.1]
df = df.sort_values('Combined Score',ascending=False)

# plot
n_rank = 10
plt.rcParams['axes.grid'] = False
plt.rcParams['figure.figsize'] = 4,3
plt.barh(width=df[:n_rank]['Combined Score'],
         y=[x.split(' (')[0] for x in df[:n_rank]['Term']],
         color='darkred')
plt.gca().invert_yaxis()
plt.xlabel('Score')
plt.title(f'Ligand genes\n({geneset})')
plt.margins(y=0.02)
plt.savefig(f"{figdir}/MoMP_LYVE1-up-genes_ligands_pathways.pdf", bbox_inches="tight",dpi=300)
plt.show()

for geneset in geneset_list:
    print(geneset)
    df = enr.results[enr.results['Gene_set'] == geneset]
    display(df.sort_values('Combined Score',ascending=False)[:30])
```

Example 40 (unknown):
```unknown
MSigDB_Hallmark_2020
```

Example 41 (unknown):
```unknown
KEGG_2021_Human
```

Example 42 (unknown):
```unknown
GO_Biological_Process_2023
```

Example 43 (unknown):
```unknown
Reactome_2022
```

Example 44 (python):
```python
adata_gex.X.data[:5]
```

Example 45 (unknown):
```unknown
array([1.0841573, 2.6877544, 1.0841573, 2.1779046, 1.0841573],
      dtype=float32)
```

Example 46 (python):
```python
n_pseudotime_bin=8
adata_gex.obs['t_scaled_binned'] = pd.cut(adata_gex.obs['t_scaled'],bins=n_pseudotime_bin)
exp_df = adata_gex[:,genes_of_interest_ligands].to_df()
exp_df['t_scaled_binned'] = adata_gex.obs['t_scaled_binned']
# groupby
exp_df = exp_df.groupby('t_scaled_binned').mean()
# min-max normalise
exp_df = (exp_df-exp_df.min())/(exp_df.max()-exp_df.min())
exp_df = exp_df.T.copy()

plt.rcParams['axes.grid'] = False
g = sns.clustermap(exp_df, col_cluster=False, row_cluster=True,
                   figsize=(3,7),
                  cmap='coolwarm',center=0,
               cbar_pos=(-0.05, 0.7, 0.03, 0.1),
                   xticklabels=True,yticklabels=True)
plt.savefig(f'{figdir}/MP_ligand_genes_heatmap.pdf',bbox_inches="tight",dpi=300)
plt.show()
# get gene order
row_order = g.dendrogram_row.reordered_ind
gene_order = list(exp_df.index[row_order])
```

Example 47 (python):
```python
# smoothing expression
from scipy.interpolate import UnivariateSpline

plt.rcParams['axes.grid'] = False
fig, ax = plt.subplots(figsize=(4,3))
cmap = 'RdPu'

adata_df = adata_gex[:,genes_of_interest_ligands].to_df()
adata_df['t_scaled'] = adata_gex.obs['t_scaled'].copy()
adata_df.sort_values('t_scaled',ascending=True,inplace=True)

pseudotime = adata_df['t_scaled'].values  # Convert to numpy array
pseudotime_reshaped = pseudotime.reshape(-1, 1)
pseudotime_plot = np.linspace(min(pseudotime), max(pseudotime), len(pseudotime))
pseudotime_plot_reshaped = pseudotime_plot.reshape(-1, 1)

gene_order_rev = gene_order.copy()
gene_order_rev.reverse()
for i, gene in enumerate(gene_order_rev):
    spl_res = UnivariateSpline(adata_df['t_scaled'], adata_df[gene], k=2)
    gene_smoothed_values = spl_res(adata_df['t_scaled'])
    # min-max normalise
    gene_smoothed_values_norm = (gene_smoothed_values-np.min(gene_smoothed_values))/(np.max(gene_smoothed_values)-np.min(gene_smoothed_values))
    norm = plt.Normalize(gene_smoothed_values_norm.min(), gene_smoothed_values_norm.max())
    sc = plt.scatter(pseudotime_plot, np.ones_like(pseudotime_plot) * (i + 1), 
                     c=gene_smoothed_values_norm, cmap=cmap, norm=norm, marker='s')

plt.colorbar(sc, label='Smoothed gene expression')

# Remove the box around the plot
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)

# Adjust x-axis position
# ax.set_ylim(0.5, len(sorted_genes) + 0.5)

plt.xlabel('Pseudotime')
# plt.ylabel('Gene')
plt.yticks(np.arange(1, len(gene_order_rev) + 1), gene_order_rev, fontsize='xx-small', fontstyle='italic')
# plt.title(title)
# plt.savefig(fname=fname, bbox_inches='tight')
plt.show()
```

Example 48 (python):
```python
np.min(gene_smoothed_values)
```

Example 49 (unknown):
```unknown
-0.020486605442069778
```

Example 50 (unknown):
```unknown
len(gene_order)
```

Example 51 (python):
```python
pseudotime = adata_df['t_scaled'].values  # Convert to numpy array
pseudotime_reshaped = pseudotime.reshape(-1, 1)
pseudotime_plot = np.linspace(min(pseudotime), max(pseudotime), len(pseudotime))
pseudotime_plot_reshaped = pseudotime_plot.reshape(-1, 1)

nrows = 4
ncols = 6
fig,axs=plt.subplots(nrows=nrows,ncols=ncols,figsize=(20,10)) 
for i,gene in enumerate(gene_order):
    spl_res = UnivariateSpline(adata_df['t_scaled'], adata_df[gene], k=2, s=None)
    gene_smoothed_values = spl_res(adata_df['t_scaled'])
    norm = plt.Normalize(gene_smoothed_values.min(), gene_smoothed_values.max())
    
    # plot
    axs[i//ncols,i%ncols].scatter(x=adata_df['t_scaled'],y=adata_df[gene])
    axs[i//ncols,i%ncols].plot(adata_df['t_scaled'],gene_smoothed_values,color='red')
    axs[i//ncols,i%ncols].set_title(gene)
    axs[i//ncols,i%ncols].set_xlabel('pseudotime')
    axs[i//ncols,i%ncols].set_ylabel('smoothed gene expression')
# fig.suptitle('')
fig.tight_layout()
```

Example 52 (unknown):
```unknown
* associated genes
* associated regions
* espressed genes
* expressed regions
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Supplementary_materials/Macrophages/Mo-MP_no1_MultiVI_prep-CPU.html

**Contents:**
- Variables¶
- Read in¶
- Prepare paired anndata¶
- Prepare RNA data¶

36601 rows × 3 columns

**Examples:**

Example 1 (python):
```python
import scanpy as sc
import pandas as pd
import numpy as np
import anndata
import scipy

from scipy.io import mmwrite, mmread
from scipy.sparse import csr_matrix

import sctk
```

Example 2 (python):
```python
import session_info
session_info.show()
```

Example 3 (unknown):
```unknown
-----
anndata             0.8.0
numpy               1.21.5
pandas              1.3.5
scanpy              1.8.2
scipy               1.7.3
sctk                NA
session_info        1.0.0
-----
```

Example 4 (python):
```python
PIL                 9.1.1
asciitree           NA
backcall            0.2.0
beta_ufunc          NA
binom_ufunc         NA
cached_property     1.5.1
cairo               1.20.0
cffi                1.14.4
cloudpickle         1.6.0
colorama            0.4.4
cycler              0.10.0
cython_runtime      NA
cytoolz             0.11.2
dask                2022.01.0
dateutil            2.8.1
decorator           4.4.2
defusedxml          0.6.0
fasteners           0.17.3
fsspec              2022.01.0
google              NA
h5py                3.1.0
igraph              0.10.2
importlib_metadata  1.7.0
ipykernel           5.4.2
ipython_genutils    0.2.0
ipywidgets          7.6.3
jedi                0.17.2
jinja2              3.1.2
joblib              1.0.0
kiwisolver          1.3.1
leidenalg           0.9.0
llvmlite            0.35.0
louvain             0.7.0
markupsafe          2.1.1
matplotlib          3.5.2
mpl_toolkits        NA
msgpack             1.0.3
natsort             7.1.0
nbinom_ufunc        NA
numba               0.52.0
numcodecs           0.9.1
numexpr             2.7.2
numpy_groupies      0.9.13
packaging           21.3
parso               0.7.1
pexpect             4.8.0
pickleshare         0.7.5
pkg_resources       NA
plotly              5.10.0
prompt_toolkit      3.0.10
psutil              5.9.0
ptyprocess          0.7.0
pyarrow             9.0.0
pygments            2.13.0
pyparsing           2.4.7
pytz                2020.5
scvelo              0.2.4
seaborn             0.11.1
setuptools_scm      NA
sinfo               0.3.1
six                 1.14.0
sklearn             1.0.2
sphinxcontrib       NA
statsmodels         0.12.1
storemagic          NA
tables              3.6.1
tblib               1.7.0
texttable           1.6.3
threadpoolctl       2.1.0
tlz                 0.11.2
toolz               0.11.2
tornado             6.1
traitlets           5.0.5
typing_extensions   NA
wcwidth             0.2.5
yaml                5.4.1
zarr                2.11.3
zipp                NA
zmq                 21.0.1
zope                NA
```

Example 5 (unknown):
```unknown
-----
IPython             7.19.0
jupyter_client      6.1.11
jupyter_core        4.7.0
notebook            6.2.0
-----
Python 3.7.9 | packaged by conda-forge | (default, Dec  9 2020, 21:08:20) [GCC 9.3.0]
Linux-4.15.0-158-generic-x86_64-with-debian-bullseye-sid
-----
Session information updated at 2023-03-15 10:59
```

Example 6 (python):
```python
sc.settings.set_figure_params(dpi=120)
```

Example 7 (unknown):
```unknown
cells_to_analyse = ['MonocytesMPOpos', 
                        'Monocytes', 
                        'MonocyteDerivedCells',
                        'MacrophagesCX3CR1pos', 
                        'MacrophagesTIMD4pos', 
                        'MacrophagesLYVE1pos',
                        'MacrophagesATF3pos'
                   ]

celltype_key = 'fine_grain'

region_to_analyse = ['whole sample']
```

Example 8 (unknown):
```unknown
!ls /nfs/team205/heart/anndata_objects/Foetal
```

Example 9 (python):
```python
14Feb_fromSemih		      metadata
archive			      multiome_ATAC
Feb28ObjectRaw.h5ad	      QC
Feb28ObjectRaw_Multiome.h5ad  RNA_foetal_post-scVI_raw.h5ad
Feb28ObjectRaw_Multiome.rds   RNA_foetal_raw.h5ad
Feb28ObjectRaw.obs.csv	      scVI
Feb28ObjectRaw_original.h5ad  visium_foetal_raw.h5ad
Immune
```

Example 10 (python):
```python
# rna
adata_rna = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/Feb28ObjectRaw.h5ad')
print(adata_rna.obs[celltype_key].value_counts())

# subset cells
adata_rna = adata_rna[adata_rna.obs[celltype_key].isin(cells_to_analyse)]

# subset regions
adata_rna = adata_rna[adata_rna.obs['region'].isin(region_to_analyse)]

# add modality to .var
adata_rna.var['modality']='Gene Expression'

adata_rna
```

Example 11 (unknown):
```unknown
GreatVesselAdventitialFibroblasts        61593
GreatVesselSmoothMuscleCells             22992
VentricularCardiomyocytesRightCompact    16215
VentricularCardiomyocytesLeftCompact     15591
MyocardialInterstitialFibroblasts1       13880
                                         ...  
MacrophagesTIMD4pos                        470
MonocytesMPOpos                            414
Megakaryocytes                             360
ChromaffinCells                            228
DendriticCellsMature                       184
Name: fine_grain, Length: 63, dtype: int64
```

Example 12 (unknown):
```unknown
AnnData object with n_obs × n_vars = 7905 × 36601
    obs: 'latent_RT_efficiency', 'latent_cell_probability', 'latent_scale', 'sangerID', 'combinedID', 'donor', 'region', 'age', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'donor_by_library-prep', 'multiplexed', 'SOC | status', 'SOC | log_prob_singleton', 'SOC | log_prob_doublet', 'batch_key', '_scvi_batch', 'FACSgate', 'fine_grain', 'mid_grain', 'coarse_grain', 'sex', 'week', 'trimester', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', '_scvi_labels', 'stress_score', 'hb1_score'
    var: 'gene_name_scRNA-0', 'gene_id', 'modality'
    uns: 'FACSgate_colors', '_scvi_manager_uuid', '_scvi_uuid', 'age_colors', 'cell_or_nuclei_colors', 'coarse_grain_colors', 'cycling_colors', 'donor_colors', 'fine_grain_colors', 'heart_or_greatvessels_colors', 'hvg', 'kit_10x_colors', 'modality_colors', 'neighbors', 'phase_colors', 'region_colors', 'sex_colors', 'trimester_colors', 'umap', 'week_colors'
    obsm: 'X_scVI', 'X_umap', '_scvi_extra_categorical_covs', '_scvi_extra_continuous', '_scvi_extra_continuous_covs', 'latent_gene_encoding'
    obsp: 'connectivities', 'distances'
```

Example 13 (python):
```python
pd.crosstab(adata_rna.obs[celltype_key],adata_rna.obs['region'])
```

Example 14 (python):
```python
pd.crosstab(adata_rna.obs[celltype_key],adata_rna.obs['kit_10x'])
```

Example 15 (python):
```python
pd.crosstab(adata_rna.obs[celltype_key],adata_rna.obs['cycling'])
```

Example 16 (python):
```python
adata_rna.obs['kit_10x'].value_counts()
```

Example 17 (unknown):
```unknown
5prime-v2      5775
Multiome-v1    2130
Name: kit_10x, dtype: int64
```

Example 18 (python):
```python
# atac
# ie only Multiome data
adata_atac = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/PeakMatrix/Foetal_Peaks.h5ad')

# subset cells
adata_atac = adata_atac[adata_atac.obs[celltype_key].isin(cells_to_analyse)]
print(adata_atac.obs[celltype_key].value_counts())

# subset regions
adata_atac = adata_atac[adata_atac.obs['region'].isin(region_to_analyse)]

# add modality to .var
adata_atac.var['modality']='Peaks'

adata_atac
```

Example 19 (unknown):
```unknown
MacrophagesLYVE1pos     4157
MacrophagesATF3pos       611
MacrophagesCX3CR1pos     551
MonocyteDerivedCells     386
MacrophagesTIMD4pos      303
Monocytes                117
MonocytesMPOpos           34
Name: fine_grain, dtype: int64
```

Example 20 (unknown):
```unknown
AnnData object with n_obs × n_vars = 1899 × 508040
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score'
    var: 'modality'
```

Example 21 (python):
```python
pd.crosstab(adata_atac.obs[celltype_key],adata_atac.obs['region'])
```

Example 22 (python):
```python
# shared barcodes, ie Multiome data (not all, since MultiomeRNA has a bit more data than MultiomeATAC)
barcodes_inter=list(set(adata_rna.obs_names).intersection(adata_atac.obs_names))
len(barcodes_inter)

# take shared cells
multiome_rna=adata_rna[barcodes_inter,:]
multiome_atac=adata_atac[barcodes_inter,:]

# concatenate
adata_paired=anndata.AnnData(
    X=scipy.sparse.hstack([multiome_rna.X,multiome_atac.X]).tocsr(), # concatenate sparse matrix and convert to Compressed Sparse Row format
    obs=multiome_rna.obs[['sangerID', 'combinedID', 'donor', 
                       'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 
                       'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts',
                       'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts',
                       'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt',
                       'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo',
                       'pct_counts_ribo', 'HB_score', 'multiplexed',
                        'batch_key','region','FACSgate','sex','week','trimester',
                     'fine_grain','mid_grain','coarse_grain',
                     'heart_or_greatvessels','cycling','S_score','G2M_score','phase','stress_score','hb1_score']],
    var=pd.concat([multiome_rna.var[['modality','gene_name_scRNA-0','gene_id']], \
                   multiome_atac.var[['modality']]
                  ],axis=0)
)
adata_paired
```

Example 23 (unknown):
```unknown
AnnData object with n_obs × n_vars = 1899 × 544641
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score'
    var: 'modality', 'gene_name_scRNA-0', 'gene_id'
```

Example 24 (python):
```python
# modify modality columns in .obs
adata_paired.obs['modality'] = 'Multiome'
adata_paired.obs['modality'].value_counts()
```

Example 25 (unknown):
```unknown
Multiome    1899
Name: modality, dtype: int64
```

Example 26 (python):
```python
# save
adata_paired.write('/nfs/team205/heart/anndata_objects/Foetal/Immune/MultiVI-pre_paired_foetal_MoMP_raw.h5ad')
```

Example 27 (python):
```python
adata_paired.obs.head()
```

Example 28 (python):
```python
# select single-cell data (3p and 5p)
adata_scrna = adata_rna[adata_rna.obs['kit_10x'].isin(['3prime-v3','5prime-v2'])]
adata_scrna.obs['kit_10x'].value_counts()
```

Example 29 (unknown):
```unknown
5prime-v2    5775
Name: kit_10x, dtype: int64
```

Example 30 (python):
```python
adata_scrna.X.data[:10]
```

Example 31 (unknown):
```unknown
array([1., 5., 1., 2., 1., 3., 1., 1., 2., 1.], dtype=float32)
```

Example 32 (python):
```python
# recreate anndata
adata_scrna=anndata.AnnData(
    X=adata_scrna.X, # concatenate sparse matrix and convert to Compressed Sparse Row format
    obs=adata_scrna.obs[['sangerID', 'combinedID', 'donor', 
                       'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 
                       'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts',
                       'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts',
                       'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt',
                       'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo',
                       'pct_counts_ribo', 'HB_score', 'multiplexed',
                        'batch_key','region','FACSgate','sex','week','trimester',
                     'fine_grain','mid_grain','coarse_grain',
                     'heart_or_greatvessels','cycling','S_score','G2M_score','phase','stress_score','hb1_score']],
    var=adata_scrna.var[['modality','gene_name_scRNA-0','gene_id']]
)
adata_scrna
```

Example 33 (unknown):
```unknown
AnnData object with n_obs × n_vars = 5775 × 36601
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score'
    var: 'modality', 'gene_name_scRNA-0', 'gene_id'
```

Example 34 (python):
```python
# modify modality columns in .obs
adata_scrna.obs['modality'] = 'scRNA'
adata_scrna.obs['modality'].value_counts()
```

Example 35 (unknown):
```unknown
scRNA    5775
Name: modality, dtype: int64
```

Example 36 (python):
```python
adata_scrna.var
```

Example 37 (python):
```python
# save
adata_scrna.write('/nfs/team205/heart/anndata_objects/Foetal/Immune/MultiVI-pre_scRNA_foetal_MoMP_raw.h5ad')
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Supplementary_materials/Macrophages/Mo-MP_no4_Milo_age-enrichment.html

**Contents:**
- Read in¶
- Define neighbourhoods¶
- Differential abundance testing with GLM¶

322 rows × 12 columns

**Examples:**

Example 1 (python):
```python
import scanpy as sc
import pandas as pd
import numpy as np
import anndata
import scipy
import matplotlib.pyplot as plt
import seaborn as sns

import milopy.core as milo
import milopy.plot as milopl
```

Example 2 (python):
```python
import session_info
session_info.show()
```

Example 3 (python):
```python
-----
anndata             0.8.0
matplotlib          3.5.2
milopy              0.1.1
numpy               1.21.5
pandas              1.3.5
scanpy              1.8.2
scipy               1.7.3
seaborn             0.11.1
session_info        1.0.0
-----
```

Example 4 (python):
```python
PIL                 9.1.1
asciitree           NA
backcall            0.2.0
beta_ufunc          NA
binom_ufunc         NA
cached_property     1.5.1
cairo               1.20.0
cffi                1.14.4
cloudpickle         1.6.0
colorama            0.4.4
cycler              0.10.0
cython_runtime      NA
cytoolz             0.11.2
dask                2022.01.0
dateutil            2.8.1
decorator           4.4.2
defusedxml          0.6.0
fasteners           0.17.3
fsspec              2022.01.0
google              NA
h5py                3.1.0
igraph              0.10.2
importlib_metadata  1.7.0
ipykernel           5.4.2
ipython_genutils    0.2.0
ipywidgets          7.6.3
jedi                0.17.2
jinja2              3.1.2
joblib              1.0.0
kiwisolver          1.3.1
leidenalg           0.9.0
llvmlite            0.35.0
louvain             0.7.0
markupsafe          2.1.1
mpl_toolkits        NA
msgpack             1.0.3
natsort             7.1.0
nbinom_ufunc        NA
numba               0.52.0
numcodecs           0.9.1
numexpr             2.7.2
packaging           21.3
parso               0.7.1
pexpect             4.8.0
pickleshare         0.7.5
pkg_resources       NA
plotly              5.10.0
prompt_toolkit      3.0.10
psutil              5.9.0
ptyprocess          0.7.0
pyarrow             9.0.0
pygments            2.13.0
pyparsing           2.4.7
pytz                2020.5
rpy2                3.4.2
setuptools_scm      NA
simplejson          3.19.1
sinfo               0.3.1
six                 1.14.0
sklearn             1.0.2
sphinxcontrib       NA
statsmodels         0.12.1
storemagic          NA
tables              3.6.1
tblib               1.7.0
texttable           1.6.3
threadpoolctl       2.1.0
tlz                 0.11.2
toolz               0.11.2
tornado             6.1
traitlets           5.0.5
typing_extensions   NA
tzlocal             NA
wcwidth             0.2.5
yaml                5.4.1
zarr                2.11.3
zipp                NA
zmq                 21.0.1
zope                NA
```

Example 5 (unknown):
```unknown
-----
IPython             7.19.0
jupyter_client      6.1.11
jupyter_core        4.7.0
notebook            6.2.0
-----
Python 3.7.9 | packaged by conda-forge | (default, Dec  9 2020, 21:08:20) [GCC 9.3.0]
Linux-4.15.0-213-generic-x86_64-with-debian-bullseye-sid
-----
Session information updated at 2023-11-18 20:24
```

Example 6 (python):
```python
## r2py setup
import anndata2ri
import rpy2.rinterface_lib.callbacks
import logging
rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)

anndata2ri.activate()
```

Example 7 (python):
```python
sc.settings.set_figure_params(dpi=200)
```

Example 8 (python):
```python
adata = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/Immune/MoMP_Aug2023/MultiVI-post_foetal_Mono-MP_withContCov-CCCov_13Mar2023_feature-select-0003.h5ad')

# modify obs names
adata.obs_names = [x.replace('_paired','').replace('_expression','') for x in adata.obs_names]

# subset only macrophages
adata = adata[adata.obs['fine_grain_mod'].isin(['MacrophagesATF3pos',
                                                         'MacrophagesCX3CR1pos',
                                                         'MacrophagesLYVE1pos',
                                                         'MacrophagesTIMD4pos',
                                                         'MacrophagesYSderived',])]
adata
```

Example 9 (python):
```python
View of AnnData object with n_obs × n_vars = 6777 × 142283
    obs: 'sangerID', 'combinedID', 'donor', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'multiplexed', 'batch_key', 'region', 'FACSgate', 'sex', 'week', 'trimester', 'fine_grain', 'mid_grain', 'coarse_grain', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', 'stress_score', 'hb1_score', '_indices', '_scvi_batch', '_scvi_labels', 'fine_grain_withcycling', 't', 'seg', 'edge', 't_sd', 'milestones', 'fine_grain_mod'
    var: 'modality', 'gene_name_scRNA-0', 'gene_id', 'n_cells'
    uns: '_scvi_manager_uuid', '_scvi_uuid', 'adata_auc_gene', 'adata_auc_region', 'adata_gex', 'adata_gex_lognorm', 'adata_peak_tfidf', 'draw_graph', 'fine_grain_colors', 'fine_grain_mod_colors', 'fine_grain_withcycling_colors', 'graph', 'milestones_colors', 'neighbors', 'ppt', 'pseudotime_list', 'seg_colors', 'umap', 'week_colors'
    obsm: 'MultiVI_latent', 'X_R', 'X_draw_graph_fa', 'X_palantir', 'X_umap', 'X_umap13', '_scvi_extra_categorical_covs', '_scvi_extra_continuous_covs'
    obsp: 'connectivities', 'distances'
```

Example 10 (python):
```python
adata.X.data[:10]
```

Example 11 (unknown):
```unknown
array([1., 1., 1., 2., 1., 1., 1., 2., 3., 1.], dtype=float32)
```

Example 12 (python):
```python
sc.pl.draw_graph(adata,color='fine_grain_mod', frameon=False)
```

Example 13 (python):
```python
sc.pl.draw_graph(adata,color='kit_10x', frameon=False)
```

Example 14 (unknown):
```unknown
/home/jovyan/my-conda-envs/cellpymc/lib/python3.7/site-packages/anndata/compat/_overloaded_dict.py:106: ImplicitModificationWarning: Trying to modify attribute `._uns` of view, initializing view as actual.
  self.data[key] = value
```

Example 15 (python):
```python
%%time
## Recompute KNN graph
sc.pp.neighbors(adata, use_rep = "MultiVI_latent", n_neighbors = 50, key_added="multivi")

## Construct neighbourhoods
milo.make_nhoods(adata, neighbors_key="multivi", prop=0.05)

## Count cells in neighbourhoods
milo.count_nhoods(adata, sample_col="sangerID")
```

Example 16 (unknown):
```unknown
CPU times: user 3.43 s, sys: 7.29 s, total: 10.7 s
Wall time: 2.3 s
```

Example 17 (unknown):
```unknown
/home/jovyan/my-conda-envs/cellpymc/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.
  warnings.warn("Transforming to str index.", ImplicitModificationWarning)
```

Example 18 (python):
```python
sns.set_context("talk")
plt.rcParams["figure.figsize"] = [8,6]
plt.hist(np.array(adata.obsm['nhoods'].sum(0)).ravel(), bins=50);
plt.xlabel('# cells in neighbourhood'); plt.ylabel('# neighbourhoods');
# plt.savefig(figdir+'nhood_size_hist.pdf',bbox_inches='tight');
```

Example 19 (python):
```python
adata.uns["nhood_adata"]
```

Example 20 (unknown):
```unknown
AnnData object with n_obs × n_vars = 322 × 21
    obs: 'index_cell', 'kth_distance'
    uns: 'sample_col'
```

Example 21 (python):
```python
adata.obs["week_continuous"] = [int(x.split('W')[0]) for x in adata.obs["week"]] #rather than adata.obs["week"].cat.codes
pd.crosstab(adata.obs['week_continuous'],adata.obs['week'])
```

Example 22 (python):
```python
milo.DA_nhoods(adata, design="~ week_continuous")
```

Example 23 (python):
```python
/home/jovyan/my-conda-envs/cellpymc/lib/python3.7/site-packages/rpy2/robjects/vectors.py:980: UserWarning: R object inheriting from "POSIXct" but without attribute "tzone".
  warnings.warn('R object inheriting from "POSIXct" but without '
/home/jovyan/my-conda-envs/cellpymc/lib/python3.7/site-packages/milopy/core.py:285: FutureWarning: In a future version of pandas all arguments of concat except for the argument 'objs' will be keyword-only
  nhood_adata.obs = pd.concat([nhood_adata.obs, res], 1)
```

Example 24 (python):
```python
import milopy.utils
milopy.utils.build_nhood_graph(adata,basis='X_draw_graph_fa')
```

Example 25 (python):
```python
plt.rcParams["figure.figsize"] = [5,5]
milopl.plot_nhood_graph(adata, 
                        alpha=0.1, ## SpatialFDR level (1%) 
                        min_size=3, ## Size of smallest dot
                        title = 'DA log-Fold Change\n(FDR < 0.1)'
                       )
```

Example 26 (python):
```python
milopy.utils.annotate_nhoods(adata, anno_col='fine_grain_mod')
```

Example 27 (python):
```python
da_results = adata.uns['nhood_adata'].obs.copy()
# drop nan
da_results.dropna(axis=0,inplace=True)
# reindex
da_results.index = list(range(len(da_results)))
# add 'Nhood' column
da_results['Nhood'] = da_results.index.astype('int')+1 
da_results['Nhood'] = da_results['Nhood'].astype('str')

da_results
```

Example 28 (python):
```python
plt.hist(adata.uns['nhood_adata'].obs["nhood_annotation_frac"]);
plt.xlabel("celltype fraction")
```

Example 29 (unknown):
```unknown
Text(0.5, 0, 'celltype fraction')
```

Example 30 (unknown):
```unknown
%load_ext rpy2.ipython
```

Example 31 (unknown):
```unknown
The rpy2.ipython extension is already loaded. To reload it, use:
  %reload_ext rpy2.ipython
```

Example 32 (javascript):
```javascript
%%R
library(miloR)
library(tidyverse)
library(ggbeeswarm)

# modifying color scale
plotDAbeeswarm_mod <- function(da.res, group.by=NULL, alpha=0.1, subset.nhoods=NULL){
  if (!is.null(group.by)) {
    if (!group.by %in% colnames(da.res)) {
      stop(group.by, " is not a column in da.res. Have you forgot to run annotateNhoods(x, da.res, ", group.by,")?")
    }
    if (is.numeric(da.res[,group.by])) {
      # stop(group.by, " is a numeric variable. Please bin to use for grouping.")
    }
    da.res <- mutate(da.res, group_by = da.res[,group.by])
  } else {
    da.res <- mutate(da.res, group_by = "g1")
  }

  if (!is.factor(da.res[,"group_by"])) {
    message("Converting group_by to factor...")
    da.res <- mutate(da.res, group_by = factor(group_by, levels=unique(group_by)))
    # anno_vec <- factor(anno_vec, levels=unique(anno_vec))
  }

  if (!is.null(subset.nhoods)) {
    da.res <- da.res[subset.nhoods,]
  }
  
  # Get position with ggbeeswarm
  beeswarm_pos <- ggplot_build(
    da.res %>%
      mutate(is_signif = ifelse(SpatialFDR < alpha, 1, 0)) %>%
      arrange(group_by) %>%
      ggplot(aes(group_by, logFC)) +
      geom_quasirandom()
  )
  
  pos_x <- beeswarm_pos$data[[1]]$x
  pos_y <- beeswarm_pos$data[[1]]$y
  
  n_groups <- unique(da.res$group_by) %>% length()
  
  da.res %>%
    mutate(is_signif = ifelse(SpatialFDR < alpha, 1, 0)) %>%
    mutate(logFC_color = ifelse(is_signif==1, logFC, NA)) %>%
    arrange(group_by) %>%
    mutate(Nhood=factor(Nhood, levels=unique(Nhood))) %>%
    mutate(pos_x = pos_x, pos_y=pos_y) %>%
    ggplot(aes(pos_x, pos_y, color=logFC_color)) +
    scale_color_gradient2(low = "blue",high = "red") +
    guides(color="none") +
    xlab(group.by) + ylab("Log Fold Change") +
    scale_x_continuous(
      breaks = seq(1,n_groups),
      labels = setNames(levels(da.res$group_by), seq(1,n_groups))
      ) +
    geom_point() +
    coord_flip() +
    theme_bw(base_size=22) +
    theme(strip.text.y =  element_text(angle=0))

}
```

Example 33 (unknown):
```unknown
%%R
head(da_results)
```

Example 34 (unknown):
```unknown
index_cell kth_distance
1 BHF_F_Hea11933670_BHF_F_Hea11596623_GATTGATGTTGAGCCG-1     1.456107
2 BHF_F_Hea11933668_BHF_F_Hea11596621_GCCAGGTTCACCGGTA-1     2.753402
3  7089STDY13216927_BHF_F_Hea13242534_CACCTGTTCCCGAAGC-1     1.462260
4 BHF_F_Hea11933667_BHF_F_Hea11596620_GTTGGAGCAACAGGAT-1     1.732683
5 BHF_F_Hea11933668_BHF_F_Hea11596621_TCGCGAGGTGTTTCAC-1     1.725373
6 BHF_F_Hea11933667_BHF_F_Hea11596620_GTGAGGAGTTGTTCAC-1     1.519573
       logFC   logCPM         F       PValue          FDR   SpatialFDR
1 -0.4854753 13.40691  6.165901 1.305019e-02 7.372215e-02 7.966068e-02
2  1.1289741 13.87388 35.805067 2.304259e-09 2.473238e-07 4.691488e-07
3  0.3866579 13.73443  4.126973 4.224791e-02 1.581840e-01 1.639379e-01
4  0.7250017 13.79357 17.435932 3.012828e-05 6.929505e-04 9.498297e-04
5  0.2479639 13.22380  1.731931 1.882138e-01 3.987162e-01 4.042951e-01
6  0.3353341 13.20775  2.624008 1.053098e-01 2.969872e-01 3.030582e-01
  Nhood_size    nhood_annotation nhood_annotation_frac Nhood
1         65 MacrophagesLYVE1pos             0.9076923     1
2         77 MacrophagesLYVE1pos             0.6103896     2
3         73 MacrophagesLYVE1pos             0.9726027     3
4         89 MacrophagesLYVE1pos             0.6516854     4
5         76 MacrophagesLYVE1pos             0.9473684     5
6         75 MacrophagesLYVE1pos             0.8800000     6
```

Example 35 (unknown):
```unknown
%%R -i da_results -w 600 -h 170
# da_results = da_results %>% filter(nhood_annotation %in% c(Monocytes))
plotDAbeeswarm_mod(da_results, group.by = "nhood_annotation")
```

---
