# Developing-Heart-Local-Improved - Cell Communication

**Pages:** 4

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Fig3_SAN_niche/snMultime_no3_create-network_TF-Region-TG_IonChannels.html

**Contents:**
- 1. Read in inputs¶
- Read in and prepare anndata contains gene expression¶
- Read in peak signatures¶
- Read in scenicplus output metadata¶
- 2. DEG: SANCM vs other CMs¶
- 3. DEG: Foetal vs Adult¶
- Read in adult anndata and concatenate with foetal data¶
- Prepare anndata for pseudobulk wilcoxon test¶
- Pseudobulk wilcoxon test¶
- 4. Obtain expressing genes and accessible regions¶
- Genes¶
- Peaks¶
- Prepare scenic metadata¶
- Select interactions¶
- TF-TF network¶
- Save networks¶
- Plot network¶

What this notebook does

**Examples:**

Example 1 (python):
```python
# Parameters
working_dir = "/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM"
path_to_adata_gex = "/nfs/team205/heart/anndata_objects/Foetal/Feb28ObjectRaw_Multiome.h5ad"
path_to_adata_peak_tfidf = "/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/PeakMatrix/Foetal_Peaks_tfidf.h5ad"
path_to_DAR_1 = "/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/MarkerPeaks/SANP-vs-aCM.csv"
path_to_DAR_2 = "/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/MarkerPeaks/aCM-vs-SANP.csv"
path_to_p2g_ArchR = ""
group_key = "fine_grain_aCMs_merged"
group_a = "SinoatrialNodeCardiomyocytes"
group_b = "AtrialCardiomyocytes"
subset_object = True
geneset_list = ["MSigDB_Hallmark_2020", "KEGG_2021_Human", "GO_Biological_Process_2023"]

path_to_eRegulon_signatures_gene = "/nfs/team205/heart/pyscenic/plus/eRegulon_signatures/Cardiomyocytes_Gene_eRegulon_signatures.pkl"
path_to_eRegulon_signatures_region = "/nfs/team205/heart/pyscenic/plus/eRegulon_signatures/Cardiomyocytes_Region_eRegulon_signatures.pkl"

path_to_gene_signatures='/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM/outs/gene_signatures.pkl'
path_to_peak_signatures='/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM/outs/peak_signatures.pkl'

path_to_eRegulon_metadata = "/nfs/team205/heart/pyscenic/plus/eRegulon_metadata/Cardiomyocytes_eRegulon_metadata_filtered.csv" # filtered metadata

###### parameters for TF selection for the population ######
cell_of_intest = 'SinoatrialNodeCardiomyocytes'
groups_for_signatures = ['SinoatrialNodeCardiomyocytes']

# proportion threshold for expressing genes
exp_prop_thresh = 0.1

# proportion threshold for accesible regions
access_prop_thresh = 0.01
```

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

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

# import hotspot
```

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

Example 4 (python):
```python
-----
anndata             0.9.2
matplotlib          3.7.2
numpy               1.24.4
pandas              1.5.3
scanpy              1.9.3
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
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
scipy                       1.10.1
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 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-03-25 16:32
```

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)

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

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

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

Example 10 (python):
```python
# read in adata
# whole data
adata = sc.read_h5ad(path_to_adata_gex)

# filter and log-normalise
sc.pp.filter_genes(adata, min_cells=3)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

print(adata.X.data[:10])
adata
```

Example 11 (unknown):
```unknown
[0.41374123 0.22811745 0.22811745 0.41374123 0.22811745 0.41374123
 0.22811745 0.41374123 0.22811745 0.22811745]
```

Example 12 (unknown):
```unknown
AnnData object with n_obs × n_vars = 211145 × 34377
    obs: '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', 'fine_grain_aCMs_merged', 'fine_grain_aCMs-node_merged'
    var: 'gene_id', 'n_cells'
    uns: 'log1p'
```

Example 13 (unknown):
```unknown
peak_signatures = utils.read_pkl(path_to_peak_signatures)
peak_signatures.keys()
```

Example 14 (unknown):
```unknown
dict_keys(['SinoatrialNodeCardiomyocytes_up', 'SinoatrialNodeCardiomyocytes_down'])
```

Example 15 (python):
```python
sig_peaks = peak_signatures['SinoatrialNodeCardiomyocytes_up']
print(len(sig_peaks))
sig_peaks[:5]
```

Example 16 (unknown):
```unknown
['chr1:55555893_55556393',
 'chr12:115794713_115795213',
 'chr3:113167147_113167647',
 'chr17:65228782_65229282',
 'chr11:116485972_116486472']
```

Example 17 (python):
```python
# scenicplus output
scenic_meta = pd.read_csv(path_to_eRegulon_metadata,index_col=0)
scenic_meta['Region'] = [x.replace('-','_') for x in scenic_meta['Region']]

# select active regulons
scenic_meta = scenic_meta[['_+_' in x for x in scenic_meta['Gene_signature_name']]]
scenic_meta.head()
```

Example 18 (unknown):
```unknown
# select based on R2G_rho
rho_thresh = 0.2
scenic_meta = scenic_meta[scenic_meta['R2G_rho']>rho_thresh]
```

Example 19 (python):
```python
logFC_threshold = 1
group_mean_threshold = 0.1

adata_sub = adata[adata.obs['coarse_grain']=='Cardiomyocytes']
adata_sub = adata_sub[adata_sub.obs['kit_10x']=='Multiome-v1']

res = utils.rank_genes_groups_with_mean(adata_sub,
                                 groupby = 'fine_grain',
                                  groups = ['SinoatrialNodeCardiomyocytes'],
                                  reference = 'rest',
                                  method = 't-test',
                                  pval_cutoff = 0.05,
                                        log2fc_min=logFC_threshold
                                 )
sns.scatterplot(x=res['logfoldchanges'],y=res['group_mean'])
plt.show()
```

Example 20 (python):
```python
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/scanpy/tools/_rank_genes_groups.py:580: ImplicitModificationWarning: Trying to modify attribute `._uns` of view, initializing view as actual.
  adata.uns[key_added] = {}
```

Example 21 (unknown):
```unknown
# filter genes based on the group mean expression
res = res[res['group_mean']>group_mean_threshold]
display(res.head())
degs_sancm_vs_others = list(res['names'])
len(degs_sancm_vs_others)
```

Example 22 (unknown):
```unknown
sig_genes = degs_sancm_vs_others.copy()
```

Example 23 (unknown):
```unknown
'ISL1' in sig_genes
```

Example 24 (python):
```python
# prepare adult anndata
adata_adult = sc.read_h5ad('/nfs/team205/heart/anndata_objects/8regions/RNA_adult-8reg_full_raw_cellstate-annotated.h5ad')

# subset SAN_P_cell
adata_adult = adata_adult[adata_adult.obs['cell_state']=='SAN_P_cell']
adata_adult = adata_adult[adata_adult.obs['kit_10x']=='Multiome-v1']

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

# rename
adata_adult.var['gene_id'] =  adata_adult.var_names.copy()
adata_adult.var.set_index('gene_name-new',inplace=True)
adata_adult.var_names = adata_adult.var_names.astype('str')
adata_adult.var_names_make_unique()

adata_adult
```

Example 25 (unknown):
```unknown
[0.911473   0.31626672 0.31626672 0.31626672 0.74952203 0.911473
 0.31626672 0.31626672 0.31626672 0.31626672]
```

Example 26 (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 27 (unknown):
```unknown
AnnData object with n_obs × n_vars = 245 × 32732
    obs: 'sangerID', 'combinedID', 'donor', 'donor_type', 'region', 'region_finest', 'age', 'gender', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'flushed', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'pct_counts_ribo', 'scrublet_score', 'scrublet_leiden', 'cluster_scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'batch_key', 'leiden_scVI', 'cell_type', 'cell_state_HCAv1', 'cell_state_scNym', 'cell_state_scNym_confidence', 'cell_state', 'latent_RT_efficiency', 'latent_cell_probability', 'latent_scale', 'n_counts', '_scvi_batch', '_scvi_labels', 'clus20', 'doublet_cls', 'original_or_new', 'batch', 'scANVI_predictions', 'leiden_scArches'
    var: 'gene_name_scRNA-0-original', 'gene_name_snRNA-1-original', 'gene_name_multiome-2-original', 'gene_id'
    uns: 'age_colors', 'cell_or_nuclei_colors', 'cell_state_colors', 'cell_type_colors', 'donor_colors', 'donor_type_colors', 'facility_colors', 'flushed_colors', 'gender_colors', 'kit_10x_colors', 'leiden', 'modality_colors', 'neighbors', 'original_or_new_colors', 'region_colors', 'region_finest_colors', 'scANVI_predictions_colors', 'umap', 'log1p'
    obsm: 'X_scArches', 'X_umap'
    obsp: 'connectivities', 'distances'
```

Example 28 (python):
```python
# concatenate
## shared features
shared = list(set(adata_sub.var_names).intersection(adata_adult.var_names))

adata_sub.obs['stage'] = 'foetal'
adata_adult.obs['stage'] = 'adult'
adata_adult.obs['week'] = 'adult'
adata_adult.obs['fine_grain'] = 'adult - SAN_P_cell'

adata_conc = adata_sub[:,shared].concatenate(adata_adult[:,shared],
                                             join='inner',
                                             batch_key=None,
                                             index_unique=None
                                            )
```

Example 29 (unknown):
```unknown
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/anndata/_core/anndata.py:1763: FutureWarning: The AnnData.concatenate method is deprecated in favour of the anndata.concat function. Please use anndata.concat instead.

See the tutorial for concat at: https://anndata.readthedocs.io/en/latest/concatenation.html
  warnings.warn(
```

Example 30 (python):
```python
# combine donor and stage key
adata_conc.obs['donor_new'] = adata_conc.obs['stage'].astype('str')+'_'+adata_conc.obs['donor'].astype('str')
set(adata_conc.obs['donor_new'])
```

Example 31 (unknown):
```unknown
{'adult_A61',
 'adult_AH1',
 'adult_AH2',
 'adult_AV10',
 'adult_AV14',
 'adult_AV3',
 'foetal_C104',
 'foetal_C82',
 'foetal_C85',
 'foetal_C87',
 'foetal_C92',
 'foetal_C98',
 'foetal_Hst33',
 'foetal_Hst36',
 'foetal_Hst39',
 'foetal_Hst40',
 'foetal_Hst41'}
```

Example 32 (python):
```python
# add drimester 
donor_trimester_df = adata_conc.obs[['donor','trimester']].copy()
donor_trimester_df = donor_trimester_df.drop_duplicates()
donor_trimester_df['trimester'] = [x if x==x else 'adult' for x in donor_trimester_df['trimester']]
donor_trimester_df
```

Example 33 (python):
```python
groupby = 'donor_new'
celltypes = ['SinoatrialNodeCardiomyocytes','adult - SAN_P_cell']
n_cell_thresh = 10

# subset anndata
ad = adata_conc[adata_conc.obs['fine_grain'].isin(celltypes)]
# remove group which only has less than a threshold
vc = ad.obs[groupby].value_counts()
grp_to_drop = list(vc.index[vc < 10])
ad = ad[~ad.obs[groupby].isin(grp_to_drop)]
print(f'dropped groups: {grp_to_drop}')

# mean expression per group
exp_mean = utils.grouped_obs_mean(ad, group_key=groupby)

# make anndata
import anndata
adata_pb = anndata.AnnData(
    X=exp_mean.T,
    obs=pd.DataFrame(exp_mean.columns).set_index(0),
    var=pd.DataFrame(exp_mean.index).set_index(0)
)

# add stage key
adata_pb.obs['stage'] = [x.split('_')[0] for x in adata_pb.obs_names]
adata_pb.obs['stage'] = adata_pb.obs['stage'].astype('category')
adata_pb.obs['stage'].cat.reorder_categories(['foetal','adult'],inplace=True)

# add donor and trimester key
adata_pb.obs['donor'] = [x.split('_')[1] for x in adata_pb.obs_names]
obs = adata_pb.obs.reset_index().merge(donor_trimester_df,how='left',on='donor').set_index(0)
obs.index.name = 'index'
adata_pb.obs = obs
adata_pb.obs['trimester'] = adata_pb.obs['trimester'].astype('category')
adata_pb.obs['trimester'].cat.reorder_categories(['firsttrimester','secondtrimester','adult'],inplace=True)


print(adata_pb.obs['stage'].value_counts())
print(adata_pb)
```

Example 34 (unknown):
```unknown
dropped groups: ['foetal_C82', 'foetal_Hst36']
foetal    8
adult     6
Name: stage, dtype: int64
AnnData object with n_obs × n_vars = 14 × 30749
    obs: 'stage', 'donor', 'trimester'
```

Example 35 (python):
```python
/tmp/ipykernel_59443/3068665136.py:27: 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.
  adata_pb.obs['stage'].cat.reorder_categories(['foetal','adult'],inplace=True)
/tmp/ipykernel_59443/3068665136.py:35: 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.
  adata_pb.obs['trimester'].cat.reorder_categories(['firsttrimester','secondtrimester','adult'],inplace=True)
```

Example 36 (python):
```python
# rank test
logFC_threshold_FvsA = 1
group_mean_threshold = 0.1

res_foetal_vs_adult = utils.rank_genes_groups_with_mean(adata_pb,
                                 groupby = 'stage',
                                  groups = ['foetal'],
                                  reference = 'adult',
                                  method = 'wilcoxon',
                                  pval_cutoff = 0.1,
                                        log2fc_min=logFC_threshold_FvsA
                                 )

sns.scatterplot(x=res_foetal_vs_adult['logfoldchanges'],y=res_foetal_vs_adult['group_mean'])
plt.show()
```

Example 37 (python):
```python
res_foetal_vs_adult = res_foetal_vs_adult[res_foetal_vs_adult['group_mean']>group_mean_threshold]
degs_foetal_vs_adult = list(res_foetal_vs_adult['names'])
print(len(degs_foetal_vs_adult))
display(res_foetal_vs_adult.head())
```

Example 38 (unknown):
```unknown
res_foetal_vs_adult[res_foetal_vs_adult['names']=='CACNA1D']
```

Example 39 (unknown):
```unknown
len(degs_foetal_vs_adult)
```

Example 40 (python):
```python
# get expression proportions
exp_mat = utils.get_exp_proportion(adata,groupby='fine_grain', use_rep='X')
exp_genes = list(exp_mat.index[exp_mat[cell_of_intest]>exp_prop_thresh])
len(exp_genes)
```

Example 41 (unknown):
```unknown
'ISL1' in exp_genes
```

Example 42 (python):
```python
adata_peak = sc.read_h5ad(path_to_adata_peak_tfidf)
adata_peak
```

Example 43 (unknown):
```unknown
AnnData object with n_obs × n_vars = 167022 × 508040
    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 44 (python):
```python
# get accessible proportions
access_mat = utils.get_exp_proportion(adata_peak,groupby='fine_grain', use_rep='X')
```

Example 45 (python):
```python
y = np.array([])
y = np.append(y,1)
y
```

Example 46 (unknown):
```unknown
array([1.])
```

Example 47 (python):
```python
x = np.linspace(0, 1, num=200)
y = np.array([])
for thresh in np.linspace(0, 1, num=200):
    y = np.append(y,sum(access_mat[cell_of_intest]>thresh))
y
```

Example 48 (unknown):
```unknown
array([4.54129e+05, 2.83575e+05, 1.83387e+05, 1.33237e+05, 1.08954e+05,
       8.68990e+04, 7.14970e+04, 6.31020e+04, 5.42520e+04, 4.73650e+04,
       4.19410e+04, 3.84920e+04, 3.46050e+04, 3.13530e+04, 2.92080e+04,
       2.66460e+04, 2.44760e+04, 2.24840e+04, 2.11130e+04, 1.95230e+04,
       1.81140e+04, 1.71810e+04, 1.59230e+04, 1.48690e+04, 1.39060e+04,
       1.32050e+04, 1.24110e+04, 1.16150e+04, 1.11050e+04, 1.04250e+04,
       9.75700e+03, 9.29600e+03, 8.66600e+03, 8.09900e+03, 7.58000e+03,
       7.23600e+03, 6.75900e+03, 6.33300e+03, 6.01400e+03, 5.58700e+03,
       5.19000e+03, 4.79500e+03, 4.52400e+03, 4.19800e+03, 3.91200e+03,
       3.67600e+03, 3.37500e+03, 3.10800e+03, 2.86700e+03, 2.67500e+03,
       2.45100e+03, 2.23600e+03, 2.09100e+03, 1.92000e+03, 1.74200e+03,
       1.63200e+03, 1.45200e+03, 1.32700e+03, 1.19600e+03, 1.08900e+03,
       9.63000e+02, 8.50000e+02, 7.88000e+02, 6.94000e+02, 6.04000e+02,
       5.33000e+02, 4.77000e+02, 4.09000e+02, 3.53000e+02, 3.20000e+02,
       2.65000e+02, 2.26000e+02, 2.02000e+02, 1.74000e+02, 1.28000e+02,
       1.05000e+02, 9.60000e+01, 7.80000e+01, 5.70000e+01, 4.80000e+01,
       3.60000e+01, 2.70000e+01, 2.50000e+01, 1.90000e+01, 1.60000e+01,
       1.50000e+01, 1.30000e+01, 9.00000e+00, 7.00000e+00, 7.00000e+00,
       6.00000e+00, 3.00000e+00, 3.00000e+00, 1.00000e+00, 1.00000e+00,
       1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00, 1.00000e+00,
       1.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
       0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00])
```

Example 49 (python):
```python
plt.rcParams['figure.figsize'] = 10,10
plt.scatter(x,y)
plt.show()
```

Example 50 (unknown):
```unknown
thresh_list = [0,0.01,0.02]
access_mat[cell_of_intest]>access_prop_thresh
```

Example 51 (unknown):
```unknown
chr1:817093_817593          False
chr1:819818_820318          False
chr1:820542_821042          False
chr1:825246_825746          False
chr1:826566_827066           True
                            ...  
chrX:155820072_155820572    False
chrX:155874490_155874990    False
chrX:155880548_155881048     True
chrX:155881064_155881564     True
chrX:155881618_155882118     True
Name: SinoatrialNodeCardiomyocytes, Length: 508040, dtype: bool
```

Example 52 (unknown):
```unknown
access_peaks = list(access_mat.index[access_mat[cell_of_intest]>access_prop_thresh])
len(access_peaks)
```

Example 53 (python):
```python
# add columns to scenicplus metadata
# scenic_meta['Gene_isTF'] = scenic_meta['Gene'].isin(list(set(scenic_meta['TF'])))
scenic_meta['TF_isExp'] = scenic_meta['TF'].isin(exp_genes)
scenic_meta['TF_isSignature'] = (scenic_meta['TF'].isin(sig_genes))&(scenic_meta['TF'].isin(exp_genes))
scenic_meta['Region_isAccess'] = scenic_meta['Region'].isin(access_peaks)
scenic_meta['Region_isSignature'] = (scenic_meta['Region'].isin(sig_peaks))&(scenic_meta['Region'].isin(access_peaks))
scenic_meta.head()
```

Example 54 (python):
```python
ionchannel = list(pd.read_csv('/nfs/team205/kk18/data/geneset/HGNC/HGNC_GID177_Ion-channels.txt',sep='\t')['Approved symbol'].unique())
# gpcr = list(pd.read_csv('/nfs/team205/kk18/data/geneset/HGNC/HGNC_GID139_G-protein-coupled-receptors.txt',sep='\t')['Approved symbol'].unique())
```

Example 55 (unknown):
```unknown
target_genes = list(set(sig_genes).intersection(ionchannel))
target_genes.sort()
target_genes
```

Example 56 (unknown):
```unknown
['CACNA1D',
 'CACNA1G',
 'CACNA2D2',
 'CACNA2D3',
 'CLIC5',
 'GABRB2',
 'GRIA3',
 'GRID2',
 'GRIK1',
 'HCN1',
 'HCN4',
 'KCNJ3',
 'KCNK1',
 'KCNN2',
 'KCNT2',
 'LRRC8C',
 'RYR3',
 'TRPC5',
 'TRPM3']
```

Example 57 (python):
```python
# select interactions

# select interaction which target "target_genes"
scenic_meta_sel = scenic_meta[(scenic_meta['Gene'].isin(target_genes))&\
                               scenic_meta['TF_isExp'] &\
                               scenic_meta['Region_isAccess'] &\
                              (scenic_meta['TF_isSignature']|scenic_meta['Region_isSignature'])
                             ]

# get TFs and target genes in the selected interaction
direct_TFs = list(set(scenic_meta_sel['TF']))
targets = list(set(scenic_meta_sel['Gene']))
# remove TFs already in target genes
TFs_sel = list(set(direct_TFs) - set(targets))

print(len(TFs_sel))
print(f'direct TFs: {direct_TFs}')
print('')

# select higher level interactions 
print('Start looking for higher level TFs')
higher_level_TFs = []
while len(TFs_sel)>0:
    scenic_meta_sel2 = scenic_meta[(scenic_meta['Gene'].isin(TFs_sel))&\
                               scenic_meta['TF_isExp'] &\
                               scenic_meta['Region_isAccess'] &\
                              (scenic_meta['TF_isSignature']|scenic_meta['Region_isSignature'])
                                  ]
    higher_level_TFs = higher_level_TFs + list(set(scenic_meta_sel2['TF']))
    if len(scenic_meta_sel2)>0:
        # concatenate
        scenic_meta_sel = pd.concat([scenic_meta_sel,scenic_meta_sel2])
        scenic_meta_sel = scenic_meta_sel.drop_duplicates()
        # get TFs and target genes in the selected interaction
        TFs_sel = list(set(scenic_meta_sel['TF']))
        targets = list(set(scenic_meta_sel['Gene']))
        # remove TFs already in target genes
        TFs_sel = list(set(TFs_sel) - set(targets))
        print(f'{len(TFs_sel)}: {TFs_sel}')
    else:
        print('no more higher level TFs')
        TFs_sel = []

scenic_meta_sel['direct_TF'] = scenic_meta_sel['TF'].isin(direct_TFs)
higher_level_TFs = list(set(higher_level_TFs)-set(direct_TFs))

# reindex
scenic_meta_sel.index = range(len(scenic_meta_sel))

print("")
print(f'Number of selected interactions: {len(scenic_meta_sel)}')
print(f"Number of TFs: {len(set(scenic_meta_sel['TF']))}")
```

Example 58 (unknown):
```unknown
29
direct TFs: ['RARB', 'FOS', 'MEF2C', 'NFIB', 'PKNOX2', 'MECP2', 'HMGA2', 'MEIS1', 'ZNF138', 'ISL1', 'NFIC', 'TEAD1', 'NR2F1', 'TBX3', 'SMAD9', 'TRPS1', 'TBX5', 'CREM', 'NR3C1', 'NFIA', 'NR2F2', 'TBX18', 'TBX20', 'SMAD3', 'ZNF536', 'POU6F2', 'THRB', 'SP1', 'SHOX2']

Start looking for higher level TFs
19: ['SP4', 'NFIB', 'BHLHE40', 'MECP2', 'MAX', 'ZNF138', 'TEAD1', 'FOXO1', 'CREM', 'GLIS3', 'NR3C1', 'ZNF529', 'JUN', 'ETV6', 'KLF13', 'ZNF41', 'SP1', 'TEAD3', 'MITF']
15: ['MECP2', 'MAX', 'CREM', 'ZNF138', 'ZNF41', 'SP4', 'TEAD1', 'NFIB', 'NR3C1', 'ZNF529', 'BHLHE40', 'JUN', 'KLF13', 'SP1', 'TEAD3']
no more higher level TFs

Number of selected interactions: 357
Number of TFs: 41
```

Example 59 (python):
```python
# get TF-TG dictionary
network = {}
for tf in set(scenic_meta_sel['TF']):
    df_sub = scenic_meta_sel[scenic_meta_sel['TF']==tf]
    network[tf] = list(set(df_sub['Gene']))
print(len(network.keys()))
```

Example 60 (python):
```python
# add metrics
scenic_meta_sel['TF_n_linked_region'] = np.nan
scenic_meta_sel['TF_n_target_genes'] = np.nan
scenic_meta_sel['TF_n_target_sig_regions'] = np.nan

for tf in set(scenic_meta_sel['TF']):
    tf_ind = scenic_meta_sel[scenic_meta_sel['TF']==tf].index
    targets = list(set(scenic_meta_sel.loc[tf_ind,'Gene']))
    scenic_meta_sel.loc[tf_ind,'TF_n_target_genes'] = len(targets)
    regions = list(set(scenic_meta_sel.loc[tf_ind,'Region']))
    regions = list(set(regions).intersection(sig_peaks))
    scenic_meta_sel.loc[tf_ind,'TF_n_target_sig_regions'] = len(regions)
    if tf in list(set(scenic_meta_sel['Gene'])):
        gene_ind = scenic_meta_sel[scenic_meta_sel['Gene']==tf].index
        n_region = len(set(scenic_meta_sel.loc[gene_ind,'Region']))
        scenic_meta_sel.loc[tf_ind,'TF_n_linked_region'] = n_region
    else:
        scenic_meta_sel.loc[tf_ind,'TF_n_linked_region'] = 0
        
scenic_meta_sel.head(3)
```

Example 61 (unknown):
```unknown
scenic_meta_sel_TF_metrices = scenic_meta_sel[['TF','direct_TF','TF_n_target_genes','TF_n_target_sig_regions','TF_n_linked_region']]
# df = df[df['direct_TF']]
scenic_meta_sel_TF_metrices.drop_duplicates(inplace=True)
scenic_meta_sel_TF_metrices.set_index('TF',inplace=True)
scenic_meta_sel_TF_metrices.sort_values('TF_n_target_genes',ascending=False)[:30]
```

Example 62 (python):
```python
/tmp/ipykernel_59443/2002270316.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  scenic_meta_sel_TF_metrices.drop_duplicates(inplace=True)
```

Example 63 (python):
```python
# select interactions those targets are TFs
scenic_meta_sel_TF = scenic_meta_sel[scenic_meta_sel['Gene'].isin(list(set(scenic_meta_sel['TF'])))]

network_TF = {}
for tf in set(scenic_meta_sel['TF']): # for all the TFs in the network
    if tf in set(scenic_meta_sel_TF['TF']): # whether the TF has TF target
        df_sub = scenic_meta_sel_TF[scenic_meta_sel_TF['TF']==tf]
        network_TF[tf] = list(set(df_sub['Gene']))
    else:
        network_TF[tf] = []
print(len(network_TF.keys()))
```

Example 64 (unknown):
```unknown
# calculate centrality
centrality_type='eigenvector'
centrality_dict = regulon_plot.get_centrality(network_TF,
                   directed=True,
                   centrality_type=centrality_type, # or 'degree','in_dgree','out_degree'
                   out_df=False # if False output dictionary
                  )
centrality_dict
```

Example 65 (unknown):
```unknown
{'RARB': 0.07132633275333206,
 'THRB': 0.35390860130150104,
 'FOS': 1.3629970460658244e-09,
 'TBX18': 0.3101423014476137,
 'GLIS3': 0.04934281480150286,
 'SHOX2': 0.6002172132880007,
 'SP4': 2.5997105523411263e-15,
 'TBX3': 0.15969133415071213,
 'MEF2C': 0.18247366637712445,
 'TBX5': 0.17125320995089946,
 'NFIA': 0.1040669178038443,
 'MEIS1': 0.1132634510375253,
 'NR2F2': 0.17125320995089946,
 'HMGA2': 0.04213053122609575,
 'NR2F1': 0.17125320995089946,
 'PKNOX2': 0.2615518061778349,
 'ZNF536': 0.1257520531317977,
 'POU6F2': 0.10722242062882631,
 'MITF': 0.05352293965092095,
 'NFIB': 2.5997105523411263e-15,
 'SMAD9': 0.031239281230728118,
 'ETV6': 0.18163891467430646,
 'TRPS1': 0.12645377861738516,
 'BHLHE40': 2.5997105523411263e-15,
 'MAX': 2.5997105523411263e-15,
 'ISL1': 0.28192058047477336,
 'TEAD1': 2.5997105523411263e-15,
 'FOXO1': 0.05093147832518298,
 'NFIC': 0.07907813980071299,
 'SMAD3': 0.03953906990035779,
 'TBX20': 0.011392408424827785,
 'ZNF529': 2.5997105523411263e-15,
 'JUN': 2.5997105523411263e-15,
 'KLF13': 2.5997105523411263e-15,
 'ZNF41': 2.5997105523411263e-15,
 'TEAD3': 2.5997105523411263e-15}
```

Example 66 (python):
```python
# add to dataframe and save
scenic_meta_sel_TF_metrices.loc[centrality_dict.keys(),f'centaliry_{centrality_type}'] = list(centrality_dict.values())
scenic_meta_sel_TF_metrices[f'centaliry_{centrality_type}'].fillna(0,inplace=True)
print(len(scenic_meta_sel_TF_metrices))
display(scenic_meta_sel_TF_metrices)

# save
scenic_meta_sel_TF_metrices.to_csv('/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM/TF_selection_df_target-IonChannels.csv')
```

Example 67 (python):
```python
/tmp/ipykernel_59443/1587124067.py:2: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  scenic_meta_sel_TF_metrices.loc[centrality_dict.keys(),f'centaliry_{centrality_type}'] = list(centrality_dict.values())
/tmp/ipykernel_59443/1587124067.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  scenic_meta_sel_TF_metrices[f'centaliry_{centrality_type}'].fillna(0,inplace=True)
```

Example 68 (unknown):
```unknown
# TF-TG network
utils.save_pkl(network,
               '/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM/TF-TG-network_target-IonChannels.pkl')

# TF-TF network
utils.save_pkl(network_TF,
               '/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM/TF-TF-network_target-IonChannels.pkl')
```

Example 69 (python):
```python
import networkx as nx
```

Example 70 (unknown):
```unknown
# prepare parameters
size_scale = 7000
base_size = 200

node_size_map = {}
for n in network.keys():
    if n in centrality_dict.keys():
        node_size_map[n]=int(centrality_dict[n]*size_scale + base_size)
    else: # means no connection with other TFs
        node_size_map[n]=base_size
        

node_size_map
```

Example 71 (unknown):
```unknown
{'RARB': 699,
 'FOS': 200,
 'SP4': 200,
 'MEF2C': 1477,
 'NFIB': 200,
 'PKNOX2': 2030,
 'BHLHE40': 200,
 'MECP2': 200,
 'MAX': 200,
 'HMGA2': 494,
 'MEIS1': 992,
 'ZNF138': 200,
 'ISL1': 2173,
 'NFIC': 753,
 'TEAD1': 200,
 'FOXO1': 556,
 'NR2F1': 1398,
 'TBX3': 1317,
 'SMAD9': 418,
 'TRPS1': 1085,
 'TBX5': 1398,
 'CREM': 200,
 'GLIS3': 545,
 'NR3C1': 200,
 'ZNF529': 200,
 'JUN': 200,
 'ETV6': 1471,
 'KLF13': 200,
 'NFIA': 928,
 'NR2F2': 1398,
 'TBX18': 2370,
 'TBX20': 279,
 'SMAD3': 476,
 'ZNF536': 1080,
 'POU6F2': 950,
 'ZNF41': 200,
 'THRB': 2677,
 'SP1': 200,
 'TEAD3': 200,
 'SHOX2': 4401,
 'MITF': 574}
```

Example 72 (python):
```python
# Create a directed graph
G = nx.DiGraph()

# Add nodes and edges
for tf, targets in network.items():
    G.add_node(tf, type='TF', size=node_size_map.get(tf, 1))  # Add TF node
    for tg in targets:
        if tg not in G:
            G.add_node(tg, type='TG', size=0)  # Add TG node if not already added
        G.add_edge(tf, tg)  # Add edge

# Define tf_nodes
tf_nodes = [node for node, attr in G.nodes(data=True) if attr['type'] == 'TF']
tg_nodes = [node for node, attr in G.nodes(data=True) if attr['type'] == 'TG']
tg_nodes.sort()

# Identifying direct and indirect TFs
direct_tfs = set()
for tf, targets in network.items():
    for tg in targets:
        if tg in tg_nodes:
            direct_tfs.add(tf)

indirect_tfs = list(set(tf_nodes) - direct_tfs)  # TFs that are not in direct TFs
direct_tfs = list(direct_tfs)
indirect_tfs.sort()
direct_tfs.sort()

# Manually creating circular positions for indirect and direct TFs, and linear positions for TGs
def create_circular_positions(node_list, center_x, center_y, radius):
    angle_step = 2 * np.pi / len(node_list)
    return {node: (center_x + np.cos(i * angle_step) * radius, center_y + np.sin(i * angle_step) * radius) 
            for i, node in enumerate(node_list)}

# Create circular positions for indirect plus direct TFs
all_tfs = direct_tfs+indirect_tfs
all_tfs.sort()
tf_pos = create_circular_positions(all_tfs, 0.5, 0.5, 0.35)

# Linear positions for TGs
tg_x = [1 for _ in tg_nodes]  # All TG nodes have the same x-coordinate
tg_y = np.linspace(0.9, 0.1, len(tg_nodes))  # Evenly spaced y-coordinates
tg_pos = {tg: (x, y) for tg, x, y in zip(tg_nodes, tg_x, tg_y)}

# Combine positions
pos = {**tf_pos, **tg_pos}

# Update node colors: indirect TFs to "lightskyblue"
for node in G.nodes:
    if node in indirect_tfs:
        G.nodes[node]['color'] = 'white'
    else:
        G.nodes[node]['color'] = 'peachpuff' if G.nodes[node]['type'] == 'TF' else 'white'

# node colors
node_colors = [data['color'] for node, data in G.nodes(data=True)]

# node label colors
label_colors = dict((node,'red') if node in degs_foetal_vs_adult else (node,'black') for node in G.nodes())

# Plotting with updated positions
plt.figure(figsize=(8, 8))
nx.draw(G, pos, with_labels=False, node_color=node_colors, node_size=[G.nodes[node]['size'] for node in G.nodes()], 
        edge_color='darkgrey', alpha=0.9, 
        edgecolors="black", linewidths=1, node_shape="o")
# Draw labels with specified colors
for node, color in label_colors.items():
    plt.text(pos[node][0], pos[node][1], s=node, 
             bbox=dict(facecolor='white', alpha=0),
             horizontalalignment='center', verticalalignment='center',
             color=color)

plt.title("SANPC GRN targeting ion channel genes")
plt.savefig(f"{figdir}/SANPC_GRN_IC-genes.pdf", bbox_inches="tight",dpi=300)
plt.show()
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Fig3_SAN_niche/snMultime_no3_create-network_TF-Region-TG_SANP-signatures.html

**Contents:**
- 1. Read in inputs¶
- Read in and prepare anndata contains gene expression¶
- Read in peak signatures¶
- Read in scenicplus output metadata¶
- 2. DEG: SANCM vs other CMs¶
- Prepare anndata for pseudobulk wilcoxon test¶
- Pseudobulk wilcoxon test¶
- 4. Obtain expressing genes and accessible regions¶
- Genes¶
- Peaks¶
- Prepare scenic metadata¶
- Select interactions¶
- TF-TF network¶
- Calculate centrality¶
- Save networks¶
- Plot network¶

What this notebook does

**Examples:**

Example 1 (python):
```python
# Parameters
working_dir = "/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM"
path_to_adata_gex = "/nfs/team205/heart/anndata_objects/Foetal/Feb28ObjectRaw_Multiome.h5ad"
path_to_adata_peak_tfidf = "/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/PeakMatrix/Foetal_Peaks_tfidf.h5ad"
path_to_DAR_1 = "/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/MarkerPeaks/SANP-vs-aCM.csv"
path_to_DAR_2 = "/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/MarkerPeaks/aCM-vs-SANP.csv"
path_to_p2g_ArchR = ""
group_key = "fine_grain_aCMs_merged"
group_a = "SinoatrialNodeCardiomyocytes"
group_b = "AtrialCardiomyocytes"
subset_object = True
geneset_list = ["MSigDB_Hallmark_2020", "KEGG_2021_Human", "GO_Biological_Process_2023"]

path_to_eRegulon_signatures_gene = "/nfs/team205/heart/pyscenic/plus/eRegulon_signatures/Cardiomyocytes_Gene_eRegulon_signatures.pkl"
path_to_eRegulon_signatures_region = "/nfs/team205/heart/pyscenic/plus/eRegulon_signatures/Cardiomyocytes_Region_eRegulon_signatures.pkl"

path_to_gene_signatures='/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM/outs/gene_signatures.pkl'
path_to_peak_signatures='/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM/outs/peak_signatures.pkl'

path_to_eRegulon_metadata = "/nfs/team205/heart/pyscenic/plus/eRegulon_metadata/Cardiomyocytes_eRegulon_metadata_filtered.csv" # filtered metadata

###### parameters for TF selection for the population ######
cell_of_intest = 'SinoatrialNodeCardiomyocytes'
groups_for_signatures = ['SinoatrialNodeCardiomyocytes']

# proportion threshold for expressing genes
exp_prop_thresh = 0.1

# proportion threshold for accesible regions
access_prop_thresh = 0.01
```

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

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

import os
```

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

Example 4 (python):
```python
-----
anndata             0.9.2
matplotlib          3.7.2
numpy               1.24.4
pandas              1.5.3
scanpy              1.9.3
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
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
scipy                       1.10.1
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 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-03-25 16:10
```

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)

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

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

Example 9 (unknown):
```unknown
figdir = '/nfs/team205/heart/fetal_ms_figures/raw_plots_kk'
table_dir = '/nfs/team205/heart/fetal_ms_figures/tables'
```

Example 10 (python):
```python
# read in adata
# whole data
adata = sc.read_h5ad(path_to_adata_gex)

# filter and log-normalise
sc.pp.filter_genes(adata, min_cells=3)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

print(adata.X.data[:10])
adata
```

Example 11 (unknown):
```unknown
[0.41374123 0.22811745 0.22811745 0.41374123 0.22811745 0.41374123
 0.22811745 0.41374123 0.22811745 0.22811745]
```

Example 12 (unknown):
```unknown
AnnData object with n_obs × n_vars = 211145 × 34377
    obs: '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', 'fine_grain_aCMs_merged', 'fine_grain_aCMs-node_merged'
    var: 'gene_id', 'n_cells'
    uns: 'log1p'
```

Example 13 (unknown):
```unknown
peak_signatures = utils.read_pkl(path_to_peak_signatures)
peak_signatures.keys()
```

Example 14 (unknown):
```unknown
dict_keys(['SinoatrialNodeCardiomyocytes_up', 'SinoatrialNodeCardiomyocytes_down'])
```

Example 15 (python):
```python
sig_peaks = peak_signatures['SinoatrialNodeCardiomyocytes_up']
print(len(sig_peaks))
sig_peaks[:5]
```

Example 16 (unknown):
```unknown
['chr1:55555893_55556393',
 'chr12:115794713_115795213',
 'chr3:113167147_113167647',
 'chr17:65228782_65229282',
 'chr11:116485972_116486472']
```

Example 17 (python):
```python
# scenicplus output
scenic_meta = pd.read_csv(path_to_eRegulon_metadata,index_col=0)
scenic_meta['Region'] = [x.replace('-','_') for x in scenic_meta['Region']]

# select active regulons
scenic_meta = scenic_meta[['_+_' in x for x in scenic_meta['Gene_signature_name']]]
scenic_meta.head()
```

Example 18 (unknown):
```unknown
# select based on R2G_rho
rho_thresh = 0.2
scenic_meta = scenic_meta[scenic_meta['R2G_rho']>rho_thresh]
```

Example 19 (python):
```python
logFC_threshold = 1
group_mean_threshold = 0.1

adata_sub = adata[adata.obs['coarse_grain']=='Cardiomyocytes']
adata_sub = adata_sub[adata_sub.obs['kit_10x']=='Multiome-v1']

res = utils.rank_genes_groups_with_mean(adata_sub,
                                 groupby = 'fine_grain',
                                  groups = ['SinoatrialNodeCardiomyocytes'],
                                  reference = 'rest',
                                  method = 't-test',
                                  pval_cutoff = 0.05,
                                        log2fc_min=logFC_threshold
                                 )
sns.scatterplot(x=res['logfoldchanges'],y=res['group_mean'])
plt.show()
```

Example 20 (python):
```python
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/scanpy/tools/_rank_genes_groups.py:580: ImplicitModificationWarning: Trying to modify attribute `._uns` of view, initializing view as actual.
  adata.uns[key_added] = {}
```

Example 21 (unknown):
```unknown
# filter genes based on the group mean expression
res = res[res['group_mean']>group_mean_threshold]
display(res.head())
degs_sancm_vs_others = list(res['names'])
len(degs_sancm_vs_others)
```

Example 22 (unknown):
```unknown
sig_genes = degs_sancm_vs_others.copy()
```

Example 23 (unknown):
```unknown
'ISL1' in sig_genes
```

Example 24 (python):
```python
# prepare adult anndata
adata_adult = sc.read_h5ad('/nfs/team205/heart/anndata_objects/8regions/RNA_adult-8reg_full_raw_cellstate-annotated.h5ad')

# subset SAN_P_cell
adata_adult = adata_adult[adata_adult.obs['cell_state']=='SAN_P_cell']
adata_adult = adata_adult[adata_adult.obs['kit_10x']=='Multiome-v1']

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

# rename
adata_adult.var['gene_id'] =  adata_adult.var_names.copy()
adata_adult.var.set_index('gene_name-new',inplace=True)
adata_adult.var_names = adata_adult.var_names.astype('str')
adata_adult.var_names_make_unique()

adata_adult
```

Example 25 (unknown):
```unknown
[0.911473   0.31626672 0.31626672 0.31626672 0.74952203 0.911473
 0.31626672 0.31626672 0.31626672 0.31626672]
```

Example 26 (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 27 (unknown):
```unknown
AnnData object with n_obs × n_vars = 245 × 32732
    obs: 'sangerID', 'combinedID', 'donor', 'donor_type', 'region', 'region_finest', 'age', 'gender', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'flushed', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'pct_counts_ribo', 'scrublet_score', 'scrublet_leiden', 'cluster_scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'batch_key', 'leiden_scVI', 'cell_type', 'cell_state_HCAv1', 'cell_state_scNym', 'cell_state_scNym_confidence', 'cell_state', 'latent_RT_efficiency', 'latent_cell_probability', 'latent_scale', 'n_counts', '_scvi_batch', '_scvi_labels', 'clus20', 'doublet_cls', 'original_or_new', 'batch', 'scANVI_predictions', 'leiden_scArches'
    var: 'gene_name_scRNA-0-original', 'gene_name_snRNA-1-original', 'gene_name_multiome-2-original', 'gene_id'
    uns: 'age_colors', 'cell_or_nuclei_colors', 'cell_state_colors', 'cell_type_colors', 'donor_colors', 'donor_type_colors', 'facility_colors', 'flushed_colors', 'gender_colors', 'kit_10x_colors', 'leiden', 'modality_colors', 'neighbors', 'original_or_new_colors', 'region_colors', 'region_finest_colors', 'scANVI_predictions_colors', 'umap', 'log1p'
    obsm: 'X_scArches', 'X_umap'
    obsp: 'connectivities', 'distances'
```

Example 28 (python):
```python
# concatenate
## shared features
shared = list(set(adata_sub.var_names).intersection(adata_adult.var_names))

adata_sub.obs['stage'] = 'foetal'
adata_adult.obs['stage'] = 'adult'
adata_adult.obs['week'] = 'adult'
adata_adult.obs['fine_grain'] = 'adult - SAN_P_cell'

adata_conc = adata_sub[:,shared].concatenate(adata_adult[:,shared],
                                             join='inner',
                                             batch_key=None,
                                             index_unique=None
                                            )
```

Example 29 (unknown):
```unknown
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/anndata/_core/anndata.py:1763: FutureWarning: The AnnData.concatenate method is deprecated in favour of the anndata.concat function. Please use anndata.concat instead.

See the tutorial for concat at: https://anndata.readthedocs.io/en/latest/concatenation.html
  warnings.warn(
```

Example 30 (python):
```python
# combine donor and stage key
adata_conc.obs['donor_new'] = adata_conc.obs['stage'].astype('str')+'_'+adata_conc.obs['donor'].astype('str')
set(adata_conc.obs['donor_new'])
```

Example 31 (unknown):
```unknown
{'adult_A61',
 'adult_AH1',
 'adult_AH2',
 'adult_AV10',
 'adult_AV14',
 'adult_AV3',
 'foetal_C104',
 'foetal_C82',
 'foetal_C85',
 'foetal_C87',
 'foetal_C92',
 'foetal_C98',
 'foetal_Hst33',
 'foetal_Hst36',
 'foetal_Hst39',
 'foetal_Hst40',
 'foetal_Hst41'}
```

Example 32 (python):
```python
# add drimester 
donor_trimester_df = adata_conc.obs[['donor','trimester']].copy()
donor_trimester_df = donor_trimester_df.drop_duplicates()
donor_trimester_df['trimester'] = [x if x==x else 'adult' for x in donor_trimester_df['trimester']]
donor_trimester_df
```

Example 33 (python):
```python
groupby = 'donor_new'
celltypes = ['SinoatrialNodeCardiomyocytes','adult - SAN_P_cell']
n_cell_thresh = 10

# subset anndata
ad = adata_conc[adata_conc.obs['fine_grain'].isin(celltypes)]
# remove group which only has less than a threshold
vc = ad.obs[groupby].value_counts()
grp_to_drop = list(vc.index[vc < 10])
ad = ad[~ad.obs[groupby].isin(grp_to_drop)]
print(f'dropped groups: {grp_to_drop}')

# mean expression per group
exp_mean = utils.grouped_obs_mean(ad, group_key=groupby)

# make anndata
import anndata
adata_pb = anndata.AnnData(
    X=exp_mean.T,
    obs=pd.DataFrame(exp_mean.columns).set_index(0),
    var=pd.DataFrame(exp_mean.index).set_index(0)
)

# add stage key
adata_pb.obs['stage'] = [x.split('_')[0] for x in adata_pb.obs_names]
adata_pb.obs['stage'] = adata_pb.obs['stage'].astype('category')
adata_pb.obs['stage'].cat.reorder_categories(['foetal','adult'],inplace=True)

# add donor and trimester key
adata_pb.obs['donor'] = [x.split('_')[1] for x in adata_pb.obs_names]
obs = adata_pb.obs.reset_index().merge(donor_trimester_df,how='left',on='donor').set_index(0)
obs.index.name = 'index'
adata_pb.obs = obs
adata_pb.obs['trimester'] = adata_pb.obs['trimester'].astype('category')
adata_pb.obs['trimester'].cat.reorder_categories(['firsttrimester','secondtrimester','adult'],inplace=True)


print(adata_pb.obs['stage'].value_counts())
print(adata_pb)
```

Example 34 (unknown):
```unknown
dropped groups: ['foetal_C82', 'foetal_Hst36']
foetal    8
adult     6
Name: stage, dtype: int64
AnnData object with n_obs × n_vars = 14 × 30749
    obs: 'stage', 'donor', 'trimester'
```

Example 35 (python):
```python
/tmp/ipykernel_257905/3068665136.py:27: 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.
  adata_pb.obs['stage'].cat.reorder_categories(['foetal','adult'],inplace=True)
/tmp/ipykernel_257905/3068665136.py:35: 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.
  adata_pb.obs['trimester'].cat.reorder_categories(['firsttrimester','secondtrimester','adult'],inplace=True)
```

Example 36 (python):
```python
# rank test
logFC_threshold_FvsA = 1
group_mean_threshold = 0.1

res_foetal_vs_adult = utils.rank_genes_groups_with_mean(adata_pb,
                                 groupby = 'stage',
                                  groups = ['foetal'],
                                  reference = 'adult',
                                  method = 'wilcoxon',
                                  pval_cutoff = 0.1,
                                        log2fc_min=logFC_threshold_FvsA
                                 )

sns.scatterplot(x=res_foetal_vs_adult['logfoldchanges'],y=res_foetal_vs_adult['group_mean'])
plt.show()
```

Example 37 (python):
```python
res_foetal_vs_adult = res_foetal_vs_adult[res_foetal_vs_adult['group_mean']>group_mean_threshold]
degs_foetal_vs_adult = list(res_foetal_vs_adult['names'])
print(len(degs_foetal_vs_adult))
display(res_foetal_vs_adult.head())
```

Example 38 (unknown):
```unknown
res_foetal_vs_adult[res_foetal_vs_adult['names']=='CACNA1D']
```

Example 39 (unknown):
```unknown
len(degs_foetal_vs_adult)
```

Example 40 (python):
```python
# get expression proportions
exp_mat = utils.get_exp_proportion(adata,groupby='fine_grain', use_rep='X')
exp_genes = list(exp_mat.index[exp_mat[cell_of_intest]>exp_prop_thresh])
len(exp_genes)
```

Example 41 (unknown):
```unknown
'ISL1' in exp_genes
```

Example 42 (python):
```python
adata_peak = sc.read_h5ad(path_to_adata_peak_tfidf)
adata_peak
```

Example 43 (unknown):
```unknown
AnnData object with n_obs × n_vars = 167022 × 508040
    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 44 (python):
```python
# get accessible proportions
access_mat = utils.get_exp_proportion(adata_peak,groupby='fine_grain', use_rep='X')
```

Example 45 (unknown):
```unknown
access_peaks = list(access_mat.index[access_mat[cell_of_intest]>access_prop_thresh])
len(access_peaks)
```

Example 46 (python):
```python
# add columns to scenicplus metadata
# scenic_meta['Gene_isTF'] = scenic_meta['Gene'].isin(list(set(scenic_meta['TF'])))
scenic_meta['TF_isExp'] = scenic_meta['TF'].isin(exp_genes)
scenic_meta['TF_isSignature'] = (scenic_meta['TF'].isin(sig_genes))&(scenic_meta['TF'].isin(exp_genes))
scenic_meta['Region_isAccess'] = scenic_meta['Region'].isin(access_peaks)
scenic_meta['Region_isSignature'] = (scenic_meta['Region'].isin(sig_peaks))&(scenic_meta['Region'].isin(access_peaks))
scenic_meta.head()
```

Example 47 (python):
```python
# save the modified metadata
scenic_meta.to_csv(os.path.join(table_dir,'csvs/scenicplus_eRegulon_metadata_CMs_SANPC-highlighted.csv'))
```

Example 48 (python):
```python
target_genes = sig_genes.copy() # all the signature genes
# removing TFs from the list
```

Example 49 (python):
```python
# select interactions

# select interaction which target "target_genes"
scenic_meta_sel = scenic_meta[(scenic_meta['Gene'].isin(target_genes))&\
                               scenic_meta['TF_isExp'] &\
                               scenic_meta['Region_isAccess'] &\
                              (scenic_meta['TF_isSignature']|scenic_meta['Region_isSignature'])
                             ]

# get TFs and target genes in the selected interaction
direct_TFs = list(set(scenic_meta_sel['TF']))
targets = list(set(scenic_meta_sel['Gene']))
# remove TFs already in target genes
TFs_sel = list(set(direct_TFs) - set(targets))

print(len(TFs_sel))
print(f'direct TFs: {direct_TFs}')
print('')

# select higher level interactions 
print('Start looking for higher level TFs')
higher_level_TFs = []
while len(TFs_sel)>0:
    scenic_meta_sel2 = scenic_meta[(scenic_meta['Gene'].isin(TFs_sel))&\
                               scenic_meta['TF_isExp'] &\
                               scenic_meta['Region_isAccess'] &\
                              (scenic_meta['TF_isSignature']|scenic_meta['Region_isSignature'])
                                  ]
    higher_level_TFs = higher_level_TFs + list(set(scenic_meta_sel2['TF']))
    if len(scenic_meta_sel2)>0:
        # concatenate
        scenic_meta_sel = pd.concat([scenic_meta_sel,scenic_meta_sel2])
        scenic_meta_sel = scenic_meta_sel.drop_duplicates()
        # get TFs and target genes in the selected interaction
        TFs_sel = list(set(scenic_meta_sel['TF']))
        targets = list(set(scenic_meta_sel['Gene']))
        # remove TFs already in target genes
        TFs_sel = list(set(TFs_sel) - set(targets))
        print(f'{len(TFs_sel)}: {TFs_sel}')
    else:
        print('no more higher level TFs')
        TFs_sel = []

scenic_meta_sel['direct_TF'] = scenic_meta_sel['TF'].isin(direct_TFs)
higher_level_TFs = list(set(higher_level_TFs)-set(direct_TFs))

# reindex
scenic_meta_sel.index = range(len(scenic_meta_sel))

print("")
print(f'Number of selected interactions: {len(scenic_meta_sel)}')
print(f"Number of TFs: {len(set(scenic_meta_sel['TF']))}")
```

Example 50 (unknown):
```unknown
103
direct TFs: ['YY1', 'MITF', 'NFIC', 'TBX5', 'TCF4', 'ZNF519', 'NCOA3', 'SP4', 'TBX18', 'MBNL2', 'RFX3', 'BACH2', 'MECP2', 'NR3C1', 'USF2', 'SP3', 'NR2F1', 'NR4A3', 'ELF2', 'GLIS3', 'BHLHE40', 'MEF2A', 'ESRRG', 'POU6F2', 'KAT7', 'GTF2I', 'ZNF536', 'RARA', 'PBX1', 'ATF1', 'TCF12', 'TBX3', 'MEF2D', 'E2F3', 'CTCF', 'PKNOX2', 'ZNF33B', 'ZNF529', 'ZNF260', 'JUND', 'ZNF148', 'ZNF83', 'ZBTB20', 'MAX', 'TFCP2', 'GABPA', 'POU2F1', 'TFDP1', 'TGIF2', 'ZNF644', 'CREB3L2', 'NFIA', 'THRB', 'SHOX2', 'ZNF100', 'GTF2IRD1', 'JUN', 'FOS', 'TBX20', 'IRF2', 'PHF21A', 'CREB1', 'ELF1', 'SATB1', 'TCF7L1', 'ZNF382', 'MEF2C', 'ZNF248', 'NFIB', 'SMAD9', 'HMGA2', 'CLOCK', 'ZNF138', 'GATA4', 'CREM', 'ZNF41', 'NR3C2', 'SMAD4', 'KLF3', 'PPARA', 'ZNF704', 'CREB5', 'SMAD3', 'GATA6', 'MEIS2', 'PKNOX1', 'NFE2L2', 'BNC2', 'NR4A1', 'FOXO1', 'LCORL', 'ISL1', 'PURA', 'KLF12', 'TRPS1', 'NCOA1', 'TCF7L2', 'MEIS1', 'NR6A1', 'SMAD5', 'VEZF1', 'E2F8', 'ESRRB', 'PBX3', 'CUX1', 'ZNF680', 'RFX2', 'MAZ', 'ZNF283', 'SP1', 'HAND2', 'SREBF2', 'ZNF337', 'KLF6', 'NR2F2', 'RXRA', 'TEAD1', 'ZNF611', 'ETV6', 'ZNF141', 'KLF13', 'TEAD3', 'RARB']

Start looking for higher level TFs
79: ['YY1', 'TCF4', 'ZNF519', 'NCOA3', 'SP4', 'RFX3', 'BACH2', 'MECP2', 'NR3C1', 'USF2', 'SP3', 'NR4A3', 'ELF2', 'BHLHE40', 'MEF2A', 'KAT7', 'RARA', 'ATF1', 'ZNF654', 'TCF12', 'MEF2D', 'E2F3', 'CTCF', 'ZNF33B', 'ZNF529', 'ZNF260', 'JUND', 'ZNF148', 'MAX', 'TFCP2', 'TFDP1', 'POU2F1', 'TGIF2', 'ZNF644', 'ZNF100', 'GTF2IRD1', 'JUN', 'IRF2', 'PHF21A', 'CREB1', 'ELF1', 'SATB1', 'ZNF382', 'ZNF248', 'NFIB', 'CLOCK', 'ZNF138', 'CREM', 'ZNF41', 'NR3C2', 'SMAD4', 'KLF3', 'PPARA', 'ZNF704', 'GATA6', 'PKNOX1', 'NFE2L2', 'LCORL', 'PURA', 'NCOA1', 'TCF7L2', 'NR6A1', 'SMAD5', 'VEZF1', 'E2F8', 'ESRRB', 'PBX3', 'ZNF680', 'RFX2', 'ZNF283', 'MAZ', 'SP1', 'SREBF2', 'ZNF337', 'RXRA', 'TEAD1', 'ZNF141', 'KLF13', 'TEAD3']
no more higher level TFs

Number of selected interactions: 3221
Number of TFs: 124
```

Example 51 (python):
```python
# get TF-TG dictionary
network = {}
for tf in set(scenic_meta_sel['TF']):
    df_sub = scenic_meta_sel[scenic_meta_sel['TF']==tf]
    network[tf] = list(set(df_sub['Gene']))
print(len(network.keys()))
```

Example 52 (unknown):
```unknown
# signature genes which are not included in the interaction
genes_not_incl = set(sig_genes) - set(scenic_meta_sel['Gene'])
genes_not_incl
```

Example 53 (unknown):
```unknown
{'AC004540.1',
 'AC004594.1',
 'AC004784.1',
 'AC005009.2',
 'AC005699.1',
 'AC006148.1',
 'AC006230.1',
 'AC007319.1',
 'AC007364.1',
 'AC008124.1',
 'AC008591.1',
 'AC009264.1',
 'AC010201.3',
 'AC010478.1',
 'AC012078.2',
 'AC012404.1',
 'AC012409.2',
 'AC012459.1',
 'AC015631.1',
 'AC015923.1',
 'AC020687.1',
 'AC020909.2',
 'AC025159.1',
 'AC040168.1',
 'AC053527.2',
 'AC068587.4',
 'AC069228.1',
 'AC073050.1',
 'AC073389.2',
 'AC079142.1',
 'AC090572.2',
 'AC090578.1',
 'AC092353.1',
 'AC092353.2',
 'AC092691.1',
 'AC092691.3',
 'AC100803.4',
 'AC104211.1',
 'AC104596.1',
 'AC112487.1',
 'AC117464.1',
 'AC119674.1',
 'AC139718.2',
 'AC245297.3',
 'ADAM22',
 'ADAMTS3',
 'ADARB2',
 'ADORA2A-AS1',
 'AF064860.1',
 'AGBL3',
 'AKT3',
 'AL022068.1',
 'AL050309.1',
 'AL078600.1',
 'AL136456.1',
 'AL139294.1',
 'AL513166.1',
 'AL590006.1',
 'AL596087.2',
 'ANKAR',
 'AP001347.1',
 'AP005262.2',
 'AR',
 'ARHGEF7',
 'ATF7IP2',
 'AZIN1-AS1',
 'BX284613.2',
 'C2CD6',
 'C3orf70',
 'C7',
 'CALN1',
 'CASC15',
 'CFAP221',
 'CNBD1',
 'CNTNAP3B',
 'COL6A3',
 'COL6A6',
 'CRYM',
 'DACH2',
 'DBN1',
 'DCN',
 'DGKB',
 'DOCK3',
 'DPH6-DT',
 'ENTPD1-AS1',
 'ERMARD',
 'ERVMER61-1',
 'FAM66C',
 'FAM66D',
 'FBLN1',
 'FN1',
 'FRG1-DT',
 'FSIP2',
 'FSIP2-AS1',
 'FUT8',
 'GDNF-AS1',
 'GEM',
 'GPR37',
 'GRB14',
 'GRIK1',
 'GRM5',
 'GRXCR2',
 'GSTA4',
 'GSTO1',
 'HBZ',
 'HOTAIRM1',
 'HS6ST3',
 'HULC',
 'INTS6L',
 'KHDRBS2',
 'KIAA0319',
 'KLHL4',
 'LINC00648',
 'LINC00886',
 'LINC00937',
 'LINC00970',
 'LINC01060',
 'LINC01117',
 'LINC01122',
 'LINC01170',
 'LINC01184',
 'LINC01237',
 'LINC01266',
 'LINC01355',
 'LINC01619',
 'LINC01880',
 'LINC02147',
 'LINC02197',
 'LINC02208',
 'LINC02268',
 'LINC02798',
 'LRCH2',
 'LTBP4',
 'LUM',
 'MAGI2-AS3',
 'MCPH1-AS1',
 'MDGA2',
 'MEF2C-AS1',
 'MIPOL1',
 'MIR100HG',
 'MIR4300HG',
 'MMP16',
 'MRPS30-DT',
 'MTCP1',
 'MYCBP2-AS1',
 'MYLK-AS1',
 'NCAM2',
 'NLGN4X',
 'NLGN4Y',
 'NME4',
 'NOMO1',
 'NOVA1',
 'NPTN',
 'NR2F1-AS1',
 'NR2F2-AS1',
 'OGN',
 'PCLO',
 'PLSCR1',
 'POSTN',
 'PRKD1',
 'PRKG1-AS1',
 'PTCHD1-AS',
 'PTK2B',
 'PTPRE',
 'PURG',
 'RAB31',
 'RBMS3',
 'RBMS3-AS3',
 'RBP1',
 'RCAN1',
 'RFX3-AS1',
 'RGS7BP',
 'RIPK2',
 'RMDN2',
 'RPS6KA5',
 'RPS6KA6',
 'SERPINI1',
 'SESTD1',
 'SH3BGRL',
 'SLC25A21',
 'SLC2A3',
 'SLC38A6',
 'SLC7A2',
 'SPARCL1',
 'SRGAP2B',
 'ST6GALNAC5',
 'ST8SIA4',
 'STS',
 'SYNPR',
 'TAPT1-AS1',
 'TBX5-AS1',
 'TENM3-AS1',
 'TFPI',
 'TMEM232',
 'TMEM72-AS1',
 'TRPC5',
 'TSTD3',
 'TTLL7',
 'UBA6-AS1',
 'UGT2B4',
 'UNC79',
 'UPRT',
 'USP2-AS1',
 'UTY',
 'VEGFC',
 'VLDLR-AS1',
 'VSIG10',
 'Z85996.3',
 'ZDHHC15',
 'ZNF25',
 'ZNF483',
 'ZRANB2-AS2'}
```

Example 54 (python):
```python
# add metrics
scenic_meta_sel['TF_n_linked_region'] = np.nan
scenic_meta_sel['TF_n_target_genes'] = np.nan
scenic_meta_sel['TF_n_target_sig_regions'] = np.nan

for tf in set(scenic_meta_sel['TF']):
    tf_ind = scenic_meta_sel[scenic_meta_sel['TF']==tf].index
    targets = list(set(scenic_meta_sel.loc[tf_ind,'Gene']))
    scenic_meta_sel.loc[tf_ind,'TF_n_target_genes'] = len(targets)
    regions = list(set(scenic_meta_sel.loc[tf_ind,'Region']))
    regions = list(set(regions).intersection(sig_peaks))
    scenic_meta_sel.loc[tf_ind,'TF_n_target_sig_regions'] = len(regions)
    if tf in list(set(scenic_meta_sel['Gene'])):
        gene_ind = scenic_meta_sel[scenic_meta_sel['Gene']==tf].index
        n_region = len(set(scenic_meta_sel.loc[gene_ind,'Region']))
        scenic_meta_sel.loc[tf_ind,'TF_n_linked_region'] = n_region
    else:
        scenic_meta_sel.loc[tf_ind,'TF_n_linked_region'] = 0
        
scenic_meta_sel.head(3)
```

Example 55 (unknown):
```unknown
scenic_meta_sel_TF_metrices = scenic_meta_sel[['TF','direct_TF','TF_n_target_genes','TF_n_target_sig_regions','TF_n_linked_region']]
# df = df[df['direct_TF']]
scenic_meta_sel_TF_metrices.drop_duplicates(inplace=True)
scenic_meta_sel_TF_metrices.set_index('TF',inplace=True)
scenic_meta_sel_TF_metrices.sort_values('TF_n_target_genes',ascending=False)[:30]
```

Example 56 (python):
```python
/tmp/ipykernel_257905/2002270316.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  scenic_meta_sel_TF_metrices.drop_duplicates(inplace=True)
```

Example 57 (python):
```python
# select interactions those targets are TFs
scenic_meta_sel_TF = scenic_meta_sel[scenic_meta_sel['Gene'].isin(list(set(scenic_meta_sel['TF'])))]

network_TF = {}
for tf in set(scenic_meta_sel['TF']): # for all the TFs in the network
    if tf in set(scenic_meta_sel_TF['TF']): # whether the TF has TF target
        df_sub = scenic_meta_sel_TF[scenic_meta_sel_TF['TF']==tf]
        network_TF[tf] = list(set(df_sub['Gene']))
    # else:
    #    network_TF[tf] = []
print(len(network_TF.keys()))
```

Example 58 (unknown):
```unknown
centrality_type='eigenvector'
centrality_dict = regulon_plot.get_centrality(network_TF,
                   directed=True,
                   centrality_type=centrality_type, # or 'degree','in_dgree','out_degree'
                   out_df=False # if False output dictionary
                  )
centrality_dict
```

Example 59 (unknown):
```unknown
{'MITF': 0.049283207312519395,
 'THRB': 0.3258743088735697,
 'TBX5': 0.1576876663350341,
 'MBNL2': 0.07520027692231576,
 'NR2F1': 0.1576876663350341,
 'SMAD3': 0.036407046988683106,
 'ZNF536': 0.11579080941535508,
 'ETV6': 0.16725068440370758,
 'TCF7L1': 0.11413507214995011,
 'PKNOX2': 0.24083340658969826,
 'SHOX2': 0.5526719860295682,
 'CUX1': 0.036407046988683106,
 'FOXO1': 0.04689702436756745,
 'NR2F2': 0.1576876663350341,
 'HAND2': 0.07281409397736383,
 'MEF2C': 0.16801931261233685,
 'TCF4': 2.393778722528883e-15,
 'MEIS2': 0.016832573884945252,
 'SP4': 2.393778722528883e-15,
 'TBX3': 0.14704164566243016,
 'TBX18': 0.28557488505513406,
 'CREB3L2': 0.1940947133237148,
 'NFIC': 0.07281409397736383,
 'GLIS3': 0.04543420423291107,
 'TBX20': 0.010489977378886729,
 'CREB5': 0.08572932912523773,
 'BHLHE40': 2.393778722528883e-15,
 'MEF2A': 2.393778722528883e-15,
 'PBX1': 0.02212373949323264,
 'ESRRG': 0.03394955944879274,
 'POU6F2': 0.09872897151884444,
 'ISL1': 0.2595887016636425,
 'ZNF654': 2.393778722528883e-15,
 'GATA4': 0.04481257611127607,
 'BNC2': 0.21245384340382825,
 'TRPS1': 0.11643694886151107,
 'ZNF529': 2.393778722528883e-15,
 'ZNF148': 2.393778722528883e-15,
 'MAX': 2.393778722528883e-15,
 'NR4A1': 0.12005192785758158,
 'NFIA': 0.0958234267017418,
 'SMAD9': 0.02876471253689032,
 'ZBTB20': 0.12144208013029432,
 'JUN': 2.393778722528883e-15,
 'FOS': 1.255029458877223e-09,
 'KLF6': 1.255029458877223e-09,
 'PHF21A': 2.393778722528883e-15,
 'GTF2I': 0.03879322993363506,
 'MEIS1': 0.10429147156003983,
 'HMGA2': 0.03879322993363506,
 'GABPA': 0.03879322993363506,
 'RARB': 0.06567633367767987,
 'NFIB': 2.393778722528883e-15,
 'ZNF41': 2.393778722528883e-15,
 'GATA6': 2.393778722528883e-15,
 'KLF12': 4.787557445057766e-14,
 'TCF7L2': 2.393778722528883e-15,
 'ZNF611': 0.024078988891304342,
 'PBX3': 2.393778722528883e-15,
 'ZNF83': 4.787557445057766e-14,
 'TEAD1': 2.393778722528883e-15,
 'KLF13': 2.393778722528883e-15,
 'TEAD3': 2.393778722528883e-15}
```

Example 60 (python):
```python
# add to dataframe and save
scenic_meta_sel_TF_metrices.loc[centrality_dict.keys(),f'centaliry_{centrality_type}'] = list(centrality_dict.values())
scenic_meta_sel_TF_metrices[f'centaliry_{centrality_type}'].fillna(0,inplace=True)
print(len(scenic_meta_sel_TF_metrices))
display(scenic_meta_sel_TF_metrices)

# save
scenic_meta_sel_TF_metrices.to_csv('/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM/TF_selection_df_target-SANPC-signatures.csv')
```

Example 61 (python):
```python
/tmp/ipykernel_257905/3865649271.py:2: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  scenic_meta_sel_TF_metrices.loc[centrality_dict.keys(),f'centaliry_{centrality_type}'] = list(centrality_dict.values())
/tmp/ipykernel_257905/3865649271.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  scenic_meta_sel_TF_metrices[f'centaliry_{centrality_type}'].fillna(0,inplace=True)
```

Example 62 (unknown):
```unknown
# TF-TG network
utils.save_pkl(network,
               '/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM/TF-TG-network_target-SANPC-signatures.pkl')

# TF-TF network
utils.save_pkl(network_TF,
               '/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM/TF-TF-network_target-SANPC-signatures.pkl')
```

Example 63 (python):
```python
import networkx as nx
```

Example 64 (unknown):
```unknown
# prepare parameters
size_scale = 7000
base_size = 200

node_size_map = {}
for n in network.keys():
    if n in centrality_dict.keys():
        node_size_map[n]=int(centrality_dict[n]*size_scale + base_size)
    else: # means no connection with other TFs
        node_size_map[n]=base_size
        

node_size_map
```

Example 65 (unknown):
```unknown
{'YY1': 200,
 'MITF': 544,
 'NFIC': 709,
 'TBX5': 1303,
 'TCF4': 200,
 'ZNF519': 200,
 'NCOA3': 200,
 'SP4': 200,
 'TBX18': 2199,
 'MBNL2': 726,
 'RFX3': 200,
 'BACH2': 200,
 'MECP2': 200,
 'NR3C1': 200,
 'USF2': 200,
 'SP3': 200,
 'NR2F1': 1303,
 'NR4A3': 200,
 'ELF2': 200,
 'GLIS3': 518,
 'BHLHE40': 200,
 'MEF2A': 200,
 'ESRRG': 437,
 'POU6F2': 891,
 'KAT7': 200,
 'GTF2I': 471,
 'ZNF536': 1010,
 'RARA': 200,
 'PBX1': 354,
 'ATF1': 200,
 'ZNF654': 200,
 'TCF12': 200,
 'TBX3': 1229,
 'MEF2D': 200,
 'E2F3': 200,
 'CTCF': 200,
 'PKNOX2': 1885,
 'ZNF33B': 200,
 'ZNF529': 200,
 'ZNF260': 200,
 'JUND': 200,
 'ZNF148': 200,
 'ZNF83': 200,
 'ZBTB20': 1050,
 'MAX': 200,
 'TFCP2': 200,
 'GABPA': 471,
 'POU2F1': 200,
 'TFDP1': 200,
 'TGIF2': 200,
 'ZNF644': 200,
 'CREB3L2': 1558,
 'NFIA': 870,
 'THRB': 2481,
 'SHOX2': 4068,
 'ZNF100': 200,
 'GTF2IRD1': 200,
 'JUN': 200,
 'FOS': 200,
 'TBX20': 273,
 'IRF2': 200,
 'PHF21A': 200,
 'CREB1': 200,
 'ELF1': 200,
 'SATB1': 200,
 'TCF7L1': 998,
 'ZNF382': 200,
 'MEF2C': 1376,
 'ZNF248': 200,
 'NFIB': 200,
 'SMAD9': 401,
 'HMGA2': 471,
 'CLOCK': 200,
 'ZNF138': 200,
 'GATA4': 513,
 'CREM': 200,
 'ZNF41': 200,
 'NR3C2': 200,
 'SMAD4': 200,
 'KLF3': 200,
 'PPARA': 200,
 'ZNF704': 200,
 'CREB5': 800,
 'SMAD3': 454,
 'GATA6': 200,
 'MEIS2': 317,
 'PKNOX1': 200,
 'NFE2L2': 200,
 'BNC2': 1687,
 'NR4A1': 1040,
 'FOXO1': 528,
 'LCORL': 200,
 'ISL1': 2017,
 'PURA': 200,
 'KLF12': 200,
 'TRPS1': 1015,
 'NCOA1': 200,
 'TCF7L2': 200,
 'MEIS1': 930,
 'NR6A1': 200,
 'SMAD5': 200,
 'VEZF1': 200,
 'E2F8': 200,
 'ESRRB': 200,
 'PBX3': 200,
 'CUX1': 454,
 'ZNF680': 200,
 'RFX2': 200,
 'MAZ': 200,
 'ZNF283': 200,
 'SP1': 200,
 'HAND2': 709,
 'SREBF2': 200,
 'ZNF337': 200,
 'KLF6': 200,
 'NR2F2': 1303,
 'RXRA': 200,
 'TEAD1': 200,
 'ZNF611': 368,
 'ETV6': 1370,
 'ZNF141': 200,
 'KLF13': 200,
 'TEAD3': 200,
 'RARB': 659}
```

Example 66 (unknown):
```unknown
regulon_plot.draw_networkx(
    network=network_TF,
    nodes_color_base='white',
    nodes_color_map=None, # list
    nodes_size_base=400,
    nodes_size_map=node_size_map,
     spring_layout_k=1.8,
    spring_layout_scale=1,
    figsize =(7,7),
    title='SANPC TF network',
    random_seed=1,
    numpy_random_seed=100,
    pos_seed=5,
    savefig=f"{figdir}/SANPC_TF-network.pdf",
    save_dpi=300
)
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Fig3_SAN_niche/snMultime_no3_create-network_TF-Region-TG_AxonGuidance.html

**Contents:**
- 1. Read in inputs¶
- Read in and prepare anndata contains gene expression¶
- Read in peak signatures¶
- Read in scenicplus output metadata¶
- 2. DEG: SANCM vs other CMs¶
- 3. DEG: Foetal vs Adult¶
- Read in adult anndata and concatenate with foetal data¶
- Prepare anndata for pseudobulk wilcoxon test¶
- Pseudobulk wilcoxon test¶
- 4. Obtain expressing genes and accessible regions¶
- Genes¶
- Peaks¶
- Prepare scenic metadata¶
- Select interactions¶
- TF-TF network¶
- Save networks¶
- Plot network¶

What this notebook does

**Examples:**

Example 1 (python):
```python
# Parameters
working_dir = "/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM"
path_to_adata_gex = "/nfs/team205/heart/anndata_objects/Foetal/Feb28ObjectRaw_Multiome.h5ad"
path_to_adata_peak_tfidf = "/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/PeakMatrix/Foetal_Peaks_tfidf.h5ad"
path_to_DAR_1 = "/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/MarkerPeaks/SANP-vs-aCM.csv"
path_to_DAR_2 = "/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output/MarkerPeaks/aCM-vs-SANP.csv"
path_to_p2g_ArchR = ""
group_key = "fine_grain_aCMs_merged"
group_a = "SinoatrialNodeCardiomyocytes"
group_b = "AtrialCardiomyocytes"
subset_object = True
geneset_list = ["MSigDB_Hallmark_2020", "KEGG_2021_Human", "GO_Biological_Process_2023"]

path_to_eRegulon_signatures_gene = "/nfs/team205/heart/pyscenic/plus/eRegulon_signatures/Cardiomyocytes_Gene_eRegulon_signatures.pkl"
path_to_eRegulon_signatures_region = "/nfs/team205/heart/pyscenic/plus/eRegulon_signatures/Cardiomyocytes_Region_eRegulon_signatures.pkl"

path_to_gene_signatures='/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM/outs/gene_signatures.pkl'
path_to_peak_signatures='/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM/outs/peak_signatures.pkl'

path_to_eRegulon_metadata = "/nfs/team205/heart/pyscenic/plus/eRegulon_metadata/Cardiomyocytes_eRegulon_metadata_filtered.csv" # filtered metadata

###### parameters for TF selection for the population ######
cell_of_intest = 'SinoatrialNodeCardiomyocytes'
groups_for_signatures = ['SinoatrialNodeCardiomyocytes']

# proportion threshold for expressing genes
exp_prop_thresh = 0.1

# proportion threshold for accesible regions
access_prop_thresh = 0.01
```

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

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

# import hotspot
```

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

Example 4 (python):
```python
-----
anndata             0.9.2
matplotlib          3.7.2
numpy               1.24.4
pandas              1.5.3
scanpy              1.9.3
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
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
scipy                       1.10.1
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 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-03-25 16:38
```

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)

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

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

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

Example 10 (python):
```python
# read in adata
# whole data
adata = sc.read_h5ad(path_to_adata_gex)

# filter and log-normalise
sc.pp.filter_genes(adata, min_cells=3)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

print(adata.X.data[:10])
adata
```

Example 11 (unknown):
```unknown
[0.41374123 0.22811745 0.22811745 0.41374123 0.22811745 0.41374123
 0.22811745 0.41374123 0.22811745 0.22811745]
```

Example 12 (unknown):
```unknown
AnnData object with n_obs × n_vars = 211145 × 34377
    obs: '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', 'fine_grain_aCMs_merged', 'fine_grain_aCMs-node_merged'
    var: 'gene_id', 'n_cells'
    uns: 'log1p'
```

Example 13 (unknown):
```unknown
peak_signatures = utils.read_pkl(path_to_peak_signatures)
peak_signatures.keys()
```

Example 14 (unknown):
```unknown
dict_keys(['SinoatrialNodeCardiomyocytes_up', 'SinoatrialNodeCardiomyocytes_down'])
```

Example 15 (python):
```python
sig_peaks = peak_signatures['SinoatrialNodeCardiomyocytes_up']
print(len(sig_peaks))
sig_peaks[:5]
```

Example 16 (unknown):
```unknown
['chr1:55555893_55556393',
 'chr12:115794713_115795213',
 'chr3:113167147_113167647',
 'chr17:65228782_65229282',
 'chr11:116485972_116486472']
```

Example 17 (python):
```python
# scenicplus output
scenic_meta = pd.read_csv(path_to_eRegulon_metadata,index_col=0)
scenic_meta['Region'] = [x.replace('-','_') for x in scenic_meta['Region']]

# select active regulons
scenic_meta = scenic_meta[['_+_' in x for x in scenic_meta['Gene_signature_name']]]
scenic_meta.head()
```

Example 18 (unknown):
```unknown
# select based on R2G_rho
rho_thresh = 0.2
scenic_meta = scenic_meta[scenic_meta['R2G_rho']>rho_thresh]
```

Example 19 (python):
```python
logFC_threshold = 1
group_mean_threshold = 0.1

adata_sub = adata[adata.obs['coarse_grain']=='Cardiomyocytes']
adata_sub = adata_sub[adata_sub.obs['kit_10x']=='Multiome-v1']

res = utils.rank_genes_groups_with_mean(adata_sub,
                                 groupby = 'fine_grain',
                                  groups = ['SinoatrialNodeCardiomyocytes'],
                                  reference = 'rest',
                                  method = 't-test',
                                  pval_cutoff = 0.05,
                                        log2fc_min=logFC_threshold
                                 )
sns.scatterplot(x=res['logfoldchanges'],y=res['group_mean'])
plt.show()
```

Example 20 (python):
```python
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/scanpy/tools/_rank_genes_groups.py:580: ImplicitModificationWarning: Trying to modify attribute `._uns` of view, initializing view as actual.
  adata.uns[key_added] = {}
```

Example 21 (unknown):
```unknown
# filter genes based on the group mean expression
res = res[res['group_mean']>group_mean_threshold]
display(res.head())
degs_sancm_vs_others = list(res['names'])
len(degs_sancm_vs_others)
```

Example 22 (unknown):
```unknown
sig_genes = degs_sancm_vs_others.copy()
```

Example 23 (unknown):
```unknown
'ISL1' in sig_genes
```

Example 24 (python):
```python
# prepare adult anndata
adata_adult = sc.read_h5ad('/nfs/team205/heart/anndata_objects/8regions/RNA_adult-8reg_full_raw_cellstate-annotated.h5ad')

# subset SAN_P_cell
adata_adult = adata_adult[adata_adult.obs['cell_state']=='SAN_P_cell']
adata_adult = adata_adult[adata_adult.obs['kit_10x']=='Multiome-v1']

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

# rename
adata_adult.var['gene_id'] =  adata_adult.var_names.copy()
adata_adult.var.set_index('gene_name-new',inplace=True)
adata_adult.var_names = adata_adult.var_names.astype('str')
adata_adult.var_names_make_unique()

adata_adult
```

Example 25 (unknown):
```unknown
[0.911473   0.31626672 0.31626672 0.31626672 0.74952203 0.911473
 0.31626672 0.31626672 0.31626672 0.31626672]
```

Example 26 (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 27 (unknown):
```unknown
AnnData object with n_obs × n_vars = 245 × 32732
    obs: 'sangerID', 'combinedID', 'donor', 'donor_type', 'region', 'region_finest', 'age', 'gender', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'flushed', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'pct_counts_ribo', 'scrublet_score', 'scrublet_leiden', 'cluster_scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'batch_key', 'leiden_scVI', 'cell_type', 'cell_state_HCAv1', 'cell_state_scNym', 'cell_state_scNym_confidence', 'cell_state', 'latent_RT_efficiency', 'latent_cell_probability', 'latent_scale', 'n_counts', '_scvi_batch', '_scvi_labels', 'clus20', 'doublet_cls', 'original_or_new', 'batch', 'scANVI_predictions', 'leiden_scArches'
    var: 'gene_name_scRNA-0-original', 'gene_name_snRNA-1-original', 'gene_name_multiome-2-original', 'gene_id'
    uns: 'age_colors', 'cell_or_nuclei_colors', 'cell_state_colors', 'cell_type_colors', 'donor_colors', 'donor_type_colors', 'facility_colors', 'flushed_colors', 'gender_colors', 'kit_10x_colors', 'leiden', 'modality_colors', 'neighbors', 'original_or_new_colors', 'region_colors', 'region_finest_colors', 'scANVI_predictions_colors', 'umap', 'log1p'
    obsm: 'X_scArches', 'X_umap'
    obsp: 'connectivities', 'distances'
```

Example 28 (python):
```python
# concatenate
## shared features
shared = list(set(adata_sub.var_names).intersection(adata_adult.var_names))

adata_sub.obs['stage'] = 'foetal'
adata_adult.obs['stage'] = 'adult'
adata_adult.obs['week'] = 'adult'
adata_adult.obs['fine_grain'] = 'adult - SAN_P_cell'

adata_conc = adata_sub[:,shared].concatenate(adata_adult[:,shared],
                                             join='inner',
                                             batch_key=None,
                                             index_unique=None
                                            )
```

Example 29 (unknown):
```unknown
/nfs/team205/kk18/miniconda3/envs/vitro/lib/python3.8/site-packages/anndata/_core/anndata.py:1763: FutureWarning: The AnnData.concatenate method is deprecated in favour of the anndata.concat function. Please use anndata.concat instead.

See the tutorial for concat at: https://anndata.readthedocs.io/en/latest/concatenation.html
  warnings.warn(
```

Example 30 (python):
```python
# combine donor and stage key
adata_conc.obs['donor_new'] = adata_conc.obs['stage'].astype('str')+'_'+adata_conc.obs['donor'].astype('str')
set(adata_conc.obs['donor_new'])
```

Example 31 (unknown):
```unknown
{'adult_A61',
 'adult_AH1',
 'adult_AH2',
 'adult_AV10',
 'adult_AV14',
 'adult_AV3',
 'foetal_C104',
 'foetal_C82',
 'foetal_C85',
 'foetal_C87',
 'foetal_C92',
 'foetal_C98',
 'foetal_Hst33',
 'foetal_Hst36',
 'foetal_Hst39',
 'foetal_Hst40',
 'foetal_Hst41'}
```

Example 32 (python):
```python
# add drimester 
donor_trimester_df = adata_conc.obs[['donor','trimester']].copy()
donor_trimester_df = donor_trimester_df.drop_duplicates()
donor_trimester_df['trimester'] = [x if x==x else 'adult' for x in donor_trimester_df['trimester']]
donor_trimester_df
```

Example 33 (python):
```python
groupby = 'donor_new'
celltypes = ['SinoatrialNodeCardiomyocytes','adult - SAN_P_cell']
n_cell_thresh = 10

# subset anndata
ad = adata_conc[adata_conc.obs['fine_grain'].isin(celltypes)]
# remove group which only has less than a threshold
vc = ad.obs[groupby].value_counts()
grp_to_drop = list(vc.index[vc < 10])
ad = ad[~ad.obs[groupby].isin(grp_to_drop)]
print(f'dropped groups: {grp_to_drop}')

# mean expression per group
exp_mean = utils.grouped_obs_mean(ad, group_key=groupby)

# make anndata
import anndata
adata_pb = anndata.AnnData(
    X=exp_mean.T,
    obs=pd.DataFrame(exp_mean.columns).set_index(0),
    var=pd.DataFrame(exp_mean.index).set_index(0)
)

# add stage key
adata_pb.obs['stage'] = [x.split('_')[0] for x in adata_pb.obs_names]
adata_pb.obs['stage'] = adata_pb.obs['stage'].astype('category')
adata_pb.obs['stage'].cat.reorder_categories(['foetal','adult'],inplace=True)

# add donor and trimester key
adata_pb.obs['donor'] = [x.split('_')[1] for x in adata_pb.obs_names]
obs = adata_pb.obs.reset_index().merge(donor_trimester_df,how='left',on='donor').set_index(0)
obs.index.name = 'index'
adata_pb.obs = obs
adata_pb.obs['trimester'] = adata_pb.obs['trimester'].astype('category')
adata_pb.obs['trimester'].cat.reorder_categories(['firsttrimester','secondtrimester','adult'],inplace=True)


print(adata_pb.obs['stage'].value_counts())
print(adata_pb)
```

Example 34 (unknown):
```unknown
dropped groups: ['foetal_C82', 'foetal_Hst36']
foetal    8
adult     6
Name: stage, dtype: int64
AnnData object with n_obs × n_vars = 14 × 30749
    obs: 'stage', 'donor', 'trimester'
```

Example 35 (python):
```python
/tmp/ipykernel_58626/3068665136.py:27: 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.
  adata_pb.obs['stage'].cat.reorder_categories(['foetal','adult'],inplace=True)
/tmp/ipykernel_58626/3068665136.py:35: 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.
  adata_pb.obs['trimester'].cat.reorder_categories(['firsttrimester','secondtrimester','adult'],inplace=True)
```

Example 36 (python):
```python
# rank test
logFC_threshold_FvsA = 1
group_mean_threshold = 0.1

res_foetal_vs_adult = utils.rank_genes_groups_with_mean(adata_pb,
                                 groupby = 'stage',
                                  groups = ['foetal'],
                                  reference = 'adult',
                                  method = 'wilcoxon',
                                  pval_cutoff = 0.1,
                                        log2fc_min=logFC_threshold_FvsA
                                 )

sns.scatterplot(x=res_foetal_vs_adult['logfoldchanges'],y=res_foetal_vs_adult['group_mean'])
plt.show()
```

Example 37 (python):
```python
res_foetal_vs_adult = res_foetal_vs_adult[res_foetal_vs_adult['group_mean']>group_mean_threshold]
degs_foetal_vs_adult = list(res_foetal_vs_adult['names'])
print(len(degs_foetal_vs_adult))
display(res_foetal_vs_adult.head())
```

Example 38 (unknown):
```unknown
res_foetal_vs_adult[res_foetal_vs_adult['names']=='CACNA1D']
```

Example 39 (unknown):
```unknown
len(degs_foetal_vs_adult)
```

Example 40 (python):
```python
# get expression proportions
exp_mat = utils.get_exp_proportion(adata,groupby='fine_grain', use_rep='X')
exp_genes = list(exp_mat.index[exp_mat[cell_of_intest]>exp_prop_thresh])
len(exp_genes)
```

Example 41 (unknown):
```unknown
'ISL1' in exp_genes
```

Example 42 (python):
```python
adata_peak = sc.read_h5ad(path_to_adata_peak_tfidf)
adata_peak
```

Example 43 (unknown):
```unknown
AnnData object with n_obs × n_vars = 167022 × 508040
    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 44 (python):
```python
# get accessible proportions
access_mat = utils.get_exp_proportion(adata_peak,groupby='fine_grain', use_rep='X')
```

Example 45 (unknown):
```unknown
access_peaks = list(access_mat.index[access_mat[cell_of_intest]>access_prop_thresh])
len(access_peaks)
```

Example 46 (python):
```python
# add columns to scenicplus metadata
# scenic_meta['Gene_isTF'] = scenic_meta['Gene'].isin(list(set(scenic_meta['TF'])))
scenic_meta['TF_isExp'] = scenic_meta['TF'].isin(exp_genes)
scenic_meta['TF_isSignature'] = (scenic_meta['TF'].isin(sig_genes))&(scenic_meta['TF'].isin(exp_genes))
scenic_meta['Region_isAccess'] = scenic_meta['Region'].isin(access_peaks)
scenic_meta['Region_isSignature'] = (scenic_meta['Region'].isin(sig_peaks))&(scenic_meta['Region'].isin(access_peaks))
scenic_meta.head()
```

Example 47 (unknown):
```unknown
target_genes = ['SEMA3A', 'SLIT2', 'SLIT3', 'EFNB2', 'CXCL12', 'FGF10', 'VEGFC', 'NELL1']
# based on the previous analysis
```

Example 48 (python):
```python
# select interactions

# select interaction which target "target_genes"
scenic_meta_sel = scenic_meta[(scenic_meta['Gene'].isin(target_genes))&\
                               scenic_meta['TF_isExp'] &\
                               scenic_meta['Region_isAccess'] &\
                              (scenic_meta['TF_isSignature']|scenic_meta['Region_isSignature'])
                             ]

# get TFs and target genes in the selected interaction
direct_TFs = list(set(scenic_meta_sel['TF']))
targets = list(set(scenic_meta_sel['Gene']))
# remove TFs already in target genes
TFs_sel = list(set(direct_TFs) - set(targets))

print(len(TFs_sel))
print(f'direct TFs: {direct_TFs}')
print('')

# select higher level interactions 
print('Start looking for higher level TFs')
higher_level_TFs = []
while len(TFs_sel)>0:
    scenic_meta_sel2 = scenic_meta[(scenic_meta['Gene'].isin(TFs_sel))&\
                               scenic_meta['TF_isExp'] &\
                               scenic_meta['Region_isAccess'] &\
                              (scenic_meta['TF_isSignature']|scenic_meta['Region_isSignature'])
                                  ]
    higher_level_TFs = higher_level_TFs + list(set(scenic_meta_sel2['TF']))
    if len(scenic_meta_sel2)>0:
        # concatenate
        scenic_meta_sel = pd.concat([scenic_meta_sel,scenic_meta_sel2])
        scenic_meta_sel = scenic_meta_sel.drop_duplicates()
        # get TFs and target genes in the selected interaction
        TFs_sel = list(set(scenic_meta_sel['TF']))
        targets = list(set(scenic_meta_sel['Gene']))
        # remove TFs already in target genes
        TFs_sel = list(set(TFs_sel) - set(targets))
        print(f'{len(TFs_sel)}: {TFs_sel}')
    else:
        print('no more higher level TFs')
        TFs_sel = []

scenic_meta_sel['direct_TF'] = scenic_meta_sel['TF'].isin(direct_TFs)
higher_level_TFs = list(set(higher_level_TFs)-set(direct_TFs))

# reindex
scenic_meta_sel.index = range(len(scenic_meta_sel))

print("")
print(f'Number of selected interactions: {len(scenic_meta_sel)}')
print(f"Number of TFs: {len(set(scenic_meta_sel['TF']))}")
```

Example 49 (unknown):
```unknown
14
direct TFs: ['SMAD9', 'MEF2C', 'BNC2', 'NR2F1', 'NFIB', 'PKNOX2', 'THRB', 'HMGA2', 'TBX5', 'POU6F2', 'NR2F2', 'FOS', 'TRPS1', 'TBX18']

Start looking for higher level TFs
11: ['RARB', 'NFIA', 'SHOX2', 'TBX3', 'JUN', 'NFIB', 'ZNF41', 'MAX', 'ETV6', 'MITF', 'ISL1']
13: ['KLF13', 'ZNF529', 'FOXO1', 'NFIB', 'JUN', 'GLIS3', 'ZNF41', 'MAX', 'TEAD1', 'MEIS1', 'SP4', 'TEAD3', 'BHLHE40']
10: ['KLF13', 'ZNF529', 'NFIB', 'JUN', 'ZNF41', 'MAX', 'TEAD1', 'TEAD3', 'SP4', 'BHLHE40']
no more higher level TFs

Number of selected interactions: 275
Number of TFs: 33
```

Example 50 (python):
```python
# get TF-TG dictionary
network = {}
for tf in set(scenic_meta_sel['TF']):
    df_sub = scenic_meta_sel[scenic_meta_sel['TF']==tf]
    network[tf] = list(set(df_sub['Gene']))
print(len(network.keys()))
```

Example 51 (python):
```python
# add metrics
scenic_meta_sel['TF_n_linked_region'] = np.nan
scenic_meta_sel['TF_n_target_genes'] = np.nan
scenic_meta_sel['TF_n_target_sig_regions'] = np.nan

for tf in set(scenic_meta_sel['TF']):
    tf_ind = scenic_meta_sel[scenic_meta_sel['TF']==tf].index
    targets = list(set(scenic_meta_sel.loc[tf_ind,'Gene']))
    scenic_meta_sel.loc[tf_ind,'TF_n_target_genes'] = len(targets)
    regions = list(set(scenic_meta_sel.loc[tf_ind,'Region']))
    regions = list(set(regions).intersection(sig_peaks))
    scenic_meta_sel.loc[tf_ind,'TF_n_target_sig_regions'] = len(regions)
    if tf in list(set(scenic_meta_sel['Gene'])):
        gene_ind = scenic_meta_sel[scenic_meta_sel['Gene']==tf].index
        n_region = len(set(scenic_meta_sel.loc[gene_ind,'Region']))
        scenic_meta_sel.loc[tf_ind,'TF_n_linked_region'] = n_region
    else:
        scenic_meta_sel.loc[tf_ind,'TF_n_linked_region'] = 0
        
scenic_meta_sel.head(3)
```

Example 52 (unknown):
```unknown
scenic_meta_sel_TF_metrices = scenic_meta_sel[['TF','direct_TF','TF_n_target_genes','TF_n_target_sig_regions','TF_n_linked_region']]
# df = df[df['direct_TF']]
scenic_meta_sel_TF_metrices.drop_duplicates(inplace=True)
scenic_meta_sel_TF_metrices.set_index('TF',inplace=True)
scenic_meta_sel_TF_metrices.sort_values('TF_n_target_genes',ascending=False)[:30]
```

Example 53 (python):
```python
/tmp/ipykernel_58626/2002270316.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  scenic_meta_sel_TF_metrices.drop_duplicates(inplace=True)
```

Example 54 (python):
```python
# select interactions those targets are TFs
scenic_meta_sel_TF = scenic_meta_sel[scenic_meta_sel['Gene'].isin(list(set(scenic_meta_sel['TF'])))]

network_TF = {}
for tf in set(scenic_meta_sel['TF']): # for all the TFs in the network
    if tf in set(scenic_meta_sel_TF['TF']): # whether the TF has TF target
        df_sub = scenic_meta_sel_TF[scenic_meta_sel_TF['TF']==tf]
        network_TF[tf] = list(set(df_sub['Gene']))
    else:
        network_TF[tf] = []
print(len(network_TF.keys()))
```

Example 55 (unknown):
```unknown
# calculate centrality
centrality_type='eigenvector'
centrality_dict = regulon_plot.get_centrality(network_TF,
                   directed=True,
                   centrality_type=centrality_type, # or 'degree','in_dgree','out_degree'
                   out_df=False # if False output dictionary
                  )
centrality_dict
```

Example 56 (unknown):
```unknown
{'NFIA': 0.10256624163167062,
 'THRB': 0.3488051331071077,
 'TBX3': 0.15738853720889912,
 'SMAD9': 0.03078880142438707,
 'SHOX2': 0.5915619010224106,
 'ISL1': 0.2778552011352862,
 'TBX18': 0.305669956425896,
 'KLF13': 2.56222194632898e-15,
 'ZNF529': 2.56222194632898e-15,
 'PKNOX2': 0.2577801506738265,
 'NR2F1': 0.16878368729177762,
 'ETV6': 0.17901962703762447,
 'TBX5': 0.16878368729177762,
 'NR2F2': 0.16878368729177762,
 'TRPS1': 0.12463027719678846,
 'ZNF41': 2.56222194632898e-15,
 'BNC2': 0.2274036004364482,
 'GLIS3': 0.04863127661046923,
 'MEF2C': 0.1798423413704834,
 'RARB': 0.07029778563893391,
 'HMGA2': 0.04152299632772208,
 'MITF': 0.05275112280562701,
 'POU6F2': 0.10567624115934543,
 'MEIS1': 0.11163015809739468,
 'JUN': 2.56222194632898e-15,
 'SP4': 2.56222194632898e-15,
 'BHLHE40': 2.56222194632898e-15,
 'FOS': 1.3433422197969282e-09,
 'NFIB': 2.56222194632898e-15,
 'MAX': 2.56222194632898e-15,
 'TEAD1': 2.56222194632898e-15,
 'FOXO1': 0.050197031129579024,
 'TEAD3': 2.56222194632898e-15}
```

Example 57 (python):
```python
# add to dataframe and save
scenic_meta_sel_TF_metrices.loc[centrality_dict.keys(),f'centaliry_{centrality_type}'] = list(centrality_dict.values())
scenic_meta_sel_TF_metrices[f'centaliry_{centrality_type}'].fillna(0,inplace=True)
print(len(scenic_meta_sel_TF_metrices))
display(scenic_meta_sel_TF_metrices)

# save
scenic_meta_sel_TF_metrices.to_csv('/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM/TF_selection_df_target-AxonGuidances.csv')
```

Example 58 (python):
```python
/tmp/ipykernel_58626/851482954.py:2: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  scenic_meta_sel_TF_metrices.loc[centrality_dict.keys(),f'centaliry_{centrality_type}'] = list(centrality_dict.values())
/tmp/ipykernel_58626/851482954.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  scenic_meta_sel_TF_metrices[f'centaliry_{centrality_type}'].fillna(0,inplace=True)
```

Example 59 (unknown):
```unknown
# TF-TG network
utils.save_pkl(network,
               '/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM/TF-TG-network_target-AxonGuidances.pkl')

# TF-TF network
utils.save_pkl(network_TF,
               '/nfs/team205/heart/Vitro/TF_selection/foetal_SANPcell_vs_workingCM/TF-TF-network_target-AxonGuidances.pkl')
```

Example 60 (python):
```python
import networkx as nx
```

Example 61 (unknown):
```unknown
# prepare parameters
size_scale = 7000
base_size = 200

node_size_map = {}
for n in network.keys():
    if n in centrality_dict.keys():
        node_size_map[n]=int(centrality_dict[n]*size_scale + base_size)
    else: # means no connection with other TFs
        node_size_map[n]=base_size
        

node_size_map
```

Example 62 (unknown):
```unknown
{'NFIA': 917,
 'SHOX2': 4340,
 'KLF13': 200,
 'ZNF529': 200,
 'PKNOX2': 2004,
 'ZNF41': 200,
 'THRB': 2641,
 'SMAD9': 415,
 'MEF2C': 1458,
 'TBX3': 1301,
 'JUN': 200,
 'MITF': 569,
 'SP4': 200,
 'BHLHE40': 200,
 'FOS': 200,
 'TRPS1': 1072,
 'NFIB': 200,
 'MAX': 200,
 'TEAD1': 200,
 'HMGA2': 490,
 'ETV6': 1453,
 'POU6F2': 939,
 'NR2F2': 1381,
 'RARB': 692,
 'BNC2': 1791,
 'FOXO1': 551,
 'NR2F1': 1381,
 'GLIS3': 540,
 'TBX5': 1381,
 'MEIS1': 981,
 'TEAD3': 200,
 'TBX18': 2339,
 'ISL1': 2144}
```

Example 63 (python):
```python
# Create a directed graph
G = nx.DiGraph()

# Add nodes and edges
for tf, targets in network.items():
    G.add_node(tf, type='TF', size=node_size_map.get(tf, 1))  # Add TF node
    for tg in targets:
        if tg not in G:
            G.add_node(tg, type='TG', size=0)  # Add TG node if not already added
        G.add_edge(tf, tg)  # Add edge

# Define tf_nodes
tf_nodes = [node for node, attr in G.nodes(data=True) if attr['type'] == 'TF']
tg_nodes = [node for node, attr in G.nodes(data=True) if attr['type'] == 'TG']
tg_nodes.sort()

# Identifying direct and indirect TFs
direct_tfs = set()
for tf, targets in network.items():
    for tg in targets:
        if tg in tg_nodes:
            direct_tfs.add(tf)

indirect_tfs = list(set(tf_nodes) - direct_tfs)  # TFs that are not in direct TFs
direct_tfs = list(direct_tfs)
indirect_tfs.sort()
direct_tfs.sort()

# Manually creating circular positions for indirect and direct TFs, and linear positions for TGs
def create_circular_positions(node_list, center_x, center_y, radius):
    angle_step = 2 * np.pi / len(node_list)
    return {node: (center_x + np.cos(i * angle_step) * radius, center_y + np.sin(i * angle_step) * radius) 
            for i, node in enumerate(node_list)}

# Create circular positions for indirect plus direct TFs
all_tfs = direct_tfs+indirect_tfs
all_tfs.sort()
tf_pos = create_circular_positions(all_tfs, 0.5, 0.5, 0.25)

# Linear positions for TGs
tg_x = [1 for _ in tg_nodes]  # All TG nodes have the same x-coordinate
tg_y = np.linspace(0.7, 0.3, len(tg_nodes))  # Evenly spaced y-coordinates
tg_pos = {tg: (x, y) for tg, x, y in zip(tg_nodes, tg_x, tg_y)}

# Combine positions
pos = {**tf_pos, **tg_pos}

# Update node colors: indirect TFs to "lightskyblue"
for node in G.nodes:
    if node in indirect_tfs:
        G.nodes[node]['color'] = 'white'
    else:
        G.nodes[node]['color'] = 'peachpuff' if G.nodes[node]['type'] == 'TF' else 'white'

# node colors
node_colors = [data['color'] for node, data in G.nodes(data=True)]

# node label colors
label_colors = dict((node,'red') if node in degs_foetal_vs_adult else (node,'black') for node in G.nodes())

# Plotting with updated positions
plt.figure(figsize=(8, 5.5))
nx.draw(G, pos, with_labels=False, node_color=node_colors, node_size=[G.nodes[node]['size'] for node in G.nodes()], 
        edge_color='darkgrey', alpha=0.9, 
        edgecolors="black", linewidths=1, node_shape="o")
# Draw labels with specified colors
for node, color in label_colors.items():
    plt.text(pos[node][0], pos[node][1], s=node, 
             bbox=dict(facecolor='white', alpha=0),
             horizontalalignment='center', verticalalignment='center',
             color=color)

plt.title("SANPC GRN targeting axon guidance genes")
plt.savefig(f"{figdir}/SANPC_GRN_Axon-guidance-genes.pdf", bbox_inches="tight",dpi=300)
plt.show()
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Fig3_SAN_niche/cellphonedb-scnRNAseq_plot_SAN.html

**Contents:**
- Import modules¶
- Variables¶
- Read in gene expression data¶
- Read in cpdb outputs¶
- Update the fine-grained cell type name for plotting¶
- Plot¶
  - Receptor expressed in PSNS¶

5 rows × 3981 columns

**Examples:**

Example 1 (python):
```python
import scanpy as sc
import pandas as pd
pd.set_option('display.max_rows', 100)
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns

import ktplotspy as kpy

import os 
from matplotlib.patches import Rectangle
import itertools
```

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

Example 3 (python):
```python
-----
ktplotspy           0.2.4
matplotlib          3.7.5
numpy               1.24.4
pandas              2.0.3
scanpy              1.9.8
seaborn             0.13.2
session_info        1.0.0
-----
```

Example 4 (python):
```python
Bio                 1.83
PIL                 10.4.0
anndata             0.9.2
asttokens           NA
backcall            0.2.0
backports           NA
certifi             2024.12.14
charset_normalizer  3.4.1
comm                0.2.2
cycler              0.12.1
cython_runtime      NA
dateutil            2.9.0.post0
debugpy             1.8.12
decorator           5.1.1
executing           2.2.0
h5py                3.11.0
idna                3.10
importlib_metadata  NA
importlib_resources NA
ipykernel           6.29.5
jaraco              NA
jedi                0.19.2
joblib              1.4.2
kiwisolver          1.4.7
llvmlite            0.41.1
mizani              0.9.3
more_itertools      10.3.0
mpl_toolkits        NA
natsort             8.4.0
numba               0.58.1
packaging           24.2
parso               0.8.4
patsy               1.0.1
pexpect             4.9.0
pickleshare         0.7.5
pkg_resources       NA
platformdirs        4.3.6
plotnine            0.12.4
prompt_toolkit      3.0.50
psutil              6.1.1
ptyprocess          0.7.0
pure_eval           0.2.3
pycircos            NA
pydev_ipython       NA
pydevconsole        NA
pydevd              3.2.3
pydevd_file_utils   NA
pydevd_plugins      NA
pydevd_tracing      NA
pygments            2.19.1
pyparsing           3.1.4
pytz                2024.2
requests            2.32.3
scipy               1.10.1
six                 1.17.0
sklearn             1.3.2
stack_data          0.6.3
statsmodels         0.14.1
threadpoolctl       3.5.0
tornado             6.4.2
traitlets           5.14.3
typing_extensions   NA
urllib3             2.2.3
wcwidth             0.2.13
zipp                NA
zmq                 26.2.0
```

Example 5 (unknown):
```unknown
-----
IPython             8.12.3
jupyter_client      8.6.3
jupyter_core        5.7.2
-----
Python 3.8.20 | packaged by conda-forge | (default, Sep 30 2024, 17:52:49) [GCC 13.3.0]
Linux-5.15.0-126-generic-x86_64-with-glibc2.17
-----
Session information updated at 2025-04-27 11:07
```

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

Example 7 (unknown):
```unknown
figdir = '/nfs/team205/heart/fetal_ms_figures/raw_plots_kk/Aug2024'
table_dir = '/nfs/team205/heart/fetal_ms_figures/tables/Aug2024'
```

Example 8 (python):
```python
path_to_cpdb = f'/nfs/team205/heart/cellphonedb/outputs/foetal/stat/all_cells_min0p05'
path_to_adata = f'/nfs/team205/heart/cellphonedb/inputs/foetal/log_norm_counts_all.h5ad'
celltype_col = 'fine_grain'
```

Example 9 (python):
```python
# read in adata
adata = sc.read_h5ad(path_to_adata)
adata
```

Example 10 (unknown):
```unknown
AnnData object with n_obs × n_vars = 297473 × 34805
    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', 'cell_id'
    var: 'gene_name_scRNA-0', 'gene_id', 'n_cells'
    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 11 (unknown):
```unknown
timestamp = [x for x in os.listdir(path_to_cpdb) if '_means_' in x][0].split('_means_')[1].replace('.txt','')
timestamp
```

Example 12 (unknown):
```unknown
'01_26_2024_03:48:49'
```

Example 13 (python):
```python
cpdbout_means =  pd.read_table(f'{path_to_cpdb}/statistical_analysis_means_{timestamp}.txt',index_col=0)
cpdbout_pvals =  pd.read_table(f'{path_to_cpdb}/statistical_analysis_pvalues_{timestamp}.txt',index_col=0)
cpdbout_sigmeans =  pd.read_table(f'{path_to_cpdb}/statistical_analysis_significant_means_{timestamp}.txt',index_col=0)
interaction_metadata = pd.read_csv('/nfs/team205/heart/cellphonedb/customDB/interaction_meta/interaction_metadata_v5.0.0_2024-1-26.csv')
```

Example 14 (python):
```python
/tmp/ipykernel_705256/3444635095.py:1: DtypeWarning: Columns (4) have mixed types. Specify dtype option on import or set low_memory=False.
/tmp/ipykernel_705256/3444635095.py:2: DtypeWarning: Columns (4) have mixed types. Specify dtype option on import or set low_memory=False.
```

Example 15 (python):
```python
# add gene names to cpdb outputs
cpdbout_means['gene_name_1'] = list(interaction_metadata.set_index('id_cp_interaction').loc[cpdbout_means.index,'gene_name_1'])
cpdbout_means['gene_name_2'] = list(interaction_metadata.set_index('id_cp_interaction').loc[cpdbout_means.index,'gene_name_2'])
cpdbout_pvals['gene_name_1'] = list(interaction_metadata.set_index('id_cp_interaction').loc[cpdbout_pvals.index,'gene_name_1'])
cpdbout_pvals['gene_name_2'] = list(interaction_metadata.set_index('id_cp_interaction').loc[cpdbout_pvals.index,'gene_name_2'])
cpdbout_sigmeans['gene_name_1'] = list(interaction_metadata.set_index('id_cp_interaction').loc[cpdbout_sigmeans.index,'gene_name_1'])
cpdbout_sigmeans['gene_name_2'] = list(interaction_metadata.set_index('id_cp_interaction').loc[cpdbout_sigmeans.index,'gene_name_2'])
```

Example 16 (python):
```python
# insert gene_name columns to the cpdb output
# to select interactions based on the list
cpdbout_means[['gene_name_1','gene_name_2']] = interaction_metadata.set_index('id_cp_interaction').loc[cpdbout_means.index,['gene_name_1','gene_name_2']]
cpdbout_pvals[['gene_name_1','gene_name_2']] = interaction_metadata.set_index('id_cp_interaction').loc[cpdbout_pvals.index,['gene_name_1','gene_name_2']]
cpdbout_sigmeans[['gene_name_1','gene_name_2']] = interaction_metadata.set_index('id_cp_interaction').loc[cpdbout_sigmeans.index,['gene_name_1','gene_name_2']]
```

Example 17 (python):
```python
# query a gene
gene = 'ADGRG1'
mask1 = [gene in x.split(',') for x in interaction_metadata['gene_name_1']]
mask2 = [gene in x.split(',') for x in interaction_metadata['gene_name_2']]
mask = [any(x) for x in zip(mask1,mask2)]
interaction_metadata[mask][['id_cp_interaction','gene_name_1','gene_name_2']]
```

Example 18 (python):
```python
df = pd.read_csv('/nfs/team205/heart/anndata_objects/Foetal/finegrain_name_mapping.csv')
celltype_mapping = df.set_index('Full_name')['Short_name_finalised'].to_dict()

# additional for Xenium
celltype_mapping['SinoatrialNodePacemakerCells']='SANPC'
```

Example 19 (python):
```python
# update cell type name in adata
adata.obs.replace({celltype_col:celltype_mapping},inplace=True)

# update cell type name in means dataframe
updated_columns = []
for col in cpdbout_means.columns:
    for key, value in celltype_mapping.items():
        col = col.replace(key, value)
    updated_columns.append(col)
cpdbout_means.columns = updated_columns.copy()

# update cell type name in pvals dataframe
updated_columns = []
for col in cpdbout_pvals.columns:
    for key, value in celltype_mapping.items():
        col = col.replace(key, value)
    updated_columns.append(col)
cpdbout_pvals.columns = updated_columns.copy()

# update cell type name in sig_means dataframe
updated_columns = []
for col in cpdbout_sigmeans.columns:
    for key, value in celltype_mapping.items():
        col = col.replace(key, value)
    updated_columns.append(col)
cpdbout_sigmeans.columns = updated_columns.copy()
```

Example 20 (unknown):
```unknown
cpdbout_means.head()
```

Example 21 (unknown):
```unknown
receptor_genes = [
    # Netrin receptors
    "DCC", "NEO1", "UNC5A", "UNC5B", "UNC5C", "UNC5D",
    # Slit (Robo) receptors
    "ROBO1", "ROBO2", "ROBO3", "ROBO4",
    # Semaphorin receptors
    "NRP1", "NRP2",
    "PLXNA1", "PLXNA2", "PLXNA3", "PLXNA4",
    "PLXNB1", "PLXNB2", "PLXNB3",
    "PLXNC1", "PLXND1",
    # Ephrin (Eph) receptors
    "EPHA1", "EPHA2", "EPHA3", "EPHA4", "EPHA5", "EPHA6", "EPHA7", "EPHA8", "EPHA10",
    "EPHB1", "EPHB2", "EPHB3", "EPHB4", "EPHB6",
    # Neurotrophin / GDNF family receptors
    "NTRK1", "NTRK2", "NTRK3", "NGFR", "RET",
    "GFRA1", "GFRA2", "GFRA3", "GFRA4",
    # Neuregulin (ErbB) receptors
    "ERBB2", "ERBB3", "ERBB4",
    # Additional innervation‑related receptors
    "EDNRB", "LRP8", "VLDLR", "L1CAM", "NCAM1",
    # Others
    "CXCR4","CNTNAP4","FGFR2"
]
receptor_celltype = 'ParaN'

# ligand_genes = ['SEMA3A', 'SLIT2', 'SLIT3', 'EFNB2', 'CXCL12', 'FGF10', 'VEGFC', 'NELL1']
ligand_celltypes = ['SANPCHn','SANPCHd','SANPCT','aCMR']
```

Example 22 (python):
```python
def receiving_signal(
    cpdbout_means,
    cpdbout_pvals,
    receptor_celltype,
    receptor_genes,
    ligand_celltypes=None
):
    # Subset results
    
    ### get columns
    cellpair_cols = [x for x in cpdbout_means.columns if '|' in x]
    other_cols = [x for x in cpdbout_means.columns if '|' not in x]
    
    ### get dataframe which has receptor_celltype and the receptor_genes at the "b/2" side (majority will be like this)
    mask_row = (
        cpdbout_means['gene_name_2']
        .str.split(',')
        .apply(lambda names: any(g in receptor_genes for g in names))
    )
    if ligand_celltypes!=None:
        cellpair_cols_sel = [x for x in cellpair_cols if (f'|{receptor_celltype}' in x)&(x.split('|')[0] in ligand_celltypes)]
    else:
        cellpair_cols_sel = [x for x in cellpair_cols if (f'|{receptor_celltype}' in x)]
    col_sel = other_cols + cellpair_cols_sel
    means_2 = cpdbout_means.loc[mask_row, col_sel]
    pvals_2 = cpdbout_pvals.loc[mask_row, col_sel]

    ### get dataframe which has receptor_celltype and the receptor_genes at the "a/1" side
    mask_row = (
        cpdbout_means['gene_name_1']
        .str.split(',')
        .apply(lambda names: any(g in receptor_genes for g in names))
    )
    if ligand_celltypes!=None:
        cellpair_cols_sel = [x for x in cellpair_cols if (f'{receptor_celltype}|' in x)&(x.split('|')[1] in ligand_celltypes)]
    else:
        cellpair_cols_sel = [x for x in cellpair_cols if (f'{receptor_celltype}|' in x)]
    col_sel = other_cols + cellpair_cols_sel
    means_1 = cpdbout_means.loc[mask_row, col_sel]
    pvals_1 = cpdbout_pvals.loc[mask_row, col_sel]
    # remove rows which are already exist in means_2
    index_keep = list(set(means_1.index)-set(means_2.index))
    means_1 = means_1.loc[index_keep]
    pvals_1 = pvals_1.loc[index_keep]
    # reverse gene_a and gene_b
    mapper = {'partner_a':'partner_b',
              'partner_b':'partner_a',
              'gene_a':'gene_b',
              'gene_b':'gene_a',
              'receptor_a':'receptor_b',
              'receptor_b':'receptor_a',
              'directionality':'directionality_reverse'}
    means_1.rename(columns=mapper,inplace=True)
    pvals_1.rename(columns=mapper,inplace=True)

    new_columns = []
    for col in means_1.columns:
        if '|' in col:
            cell_a = col.split('|')[0]
            cell_b = col.split('|')[1]
            new_columns.append(f'{cell_b}|{cell_a}')
        else:
            new_columns.append(col)
    means_1.columns = new_columns
    pvals_1.columns = new_columns

    new_interacting_pair = []
    for ind in means_1.index:
        partner_a = means_1.loc[ind,'partner_a']
        partner_b = means_1.loc[ind,'partner_b']
        if 'complex:' in partner_a:
            partner_a = partner_a.replace('complex:','')
        else:
            partner_a = means_1.loc[ind,'gene_a']
        if 'complex:' in partner_b:
            partner_b = partner_b.replace('complex:','')
        else:
            partner_b = means_1.loc[ind,'gene_b']
        new_interacting_pair.append(f'{partner_a}_{partner_b}')
    means_1['interacting_pair'] = new_interacting_pair
    pvals_1['interacting_pair'] = new_interacting_pair

    ### concatenate
    means = pd.concat([means_2,means_1])
    pvals = pd.concat([pvals_2,pvals_1])

    ### remove NaN
    means = means.dropna(axis=0, how='all')
    means = means.dropna(axis=1, how='all')
    pvals = pvals.dropna(axis=0, how='all')
    pvals = pvals.dropna(axis=1, how='all')

    ### formatting
    means = means.reset_index()
    pvals = pvals.reset_index()
    # means.drop(labels=['gene_name_1','gene_name_2','directionality_reverse'],axis=1,inplace=True)
    # pvals.drop(labels=['gene_name_1','gene_name_2','directionality_reverse'],axis=1,inplace=True)
    means.drop(labels=['gene_name_1','gene_name_2'],axis=1,inplace=True)
    pvals.drop(labels=['gene_name_1','gene_name_2'],axis=1,inplace=True)
    
    return means, pvals
```

Example 23 (python):
```python
means,pvals = receiving_signal(
    cpdbout_means,
    cpdbout_pvals,
    receptor_celltype,
    receptor_genes,
    ligand_celltypes=['SANPC','aCMR','vCMCR','vCMTR']
)

plot = kpy.plot_cpdb(
    adata=adata,
    cell_type1='.',
    cell_type2='.',  # this means all cell-types
    means=means,
    pvals=pvals,
    alpha=0.05,
    celltype_key=celltype_col,
    genes=None,
    figsize=(5.2, 9),
    title='SAN, CMs-to-ParaN,\nall AxonG',
    highlight_size =1.5,
    standard_scale=False
)
plot.save(f"{figdir}/snRNAseq_CellPhoneDB_SAN_CM-AxonG-to-ParaN_all-ligands.pdf", bbox_inches="tight",dpi=300)
plot
```

Example 24 (unknown):
```unknown
/software/cellgen/team205/kk18/envs/cpdb/lib/python3.8/site-packages/plotnine/ggplot.py:587: PlotnineWarning: Saving 5.2 x 9 in image.
/software/cellgen/team205/kk18/envs/cpdb/lib/python3.8/site-packages/plotnine/ggplot.py:588: PlotnineWarning: Filename: /nfs/team205/heart/fetal_ms_figures/raw_plots_kk/Aug2024/snRNAseq_CellPhoneDB_SAN_CM-AxonG-to-ParaN_all-ligands.pdf
```

Example 25 (unknown):
```unknown
<Figure Size: (520 x 900)>
```

---
