# Developing-Heart-Local-Improved - Trisomy21

**Pages:** 9

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Fig5_Trisomy21/no12_GSEA_finegrain.html

**Contents:**
- Read in DEG results¶
- GSEA analysis¶
- Euploid-enriched nhoods vs Others¶

20860 rows × 10 columns

**Examples:**

Example 1 (python):
```python
import scanpy as sc
import pandas as pd
import numpy as np
import anndata
import scipy
import matplotlib.pyplot as plt
import seaborn as sns
import os
import warnings
warnings.simplefilter("ignore", UserWarning)
warnings.simplefilter(action='ignore', category=pd.errors.PerformanceWarning)
```

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

Example 3 (unknown):
```unknown
os.getcwd()
```

Example 4 (unknown):
```unknown
'/nfs/team205/kk18/notebooks/Foetal/Trisomy21/RNA/milo_Mar2025_refine-nhoods'
```

Example 5 (python):
```python
deg_res = pd.read_csv('/nfs/team205/kk18/notebooks/Foetal/Trisomy21/RNA/milo_Mar2025_refine-nhoods/milo_nhood-reassigned_pseudobulk_edgeR_res_tri21-nhoods_vs_euploid-nhoods_finegrain.csv',
                          index_col=0)
print(deg_res.shape)
deg_res.head()
```

Example 6 (python):
```python
import blitzgsea as blitz
```

Example 7 (python):
```python
%%time
geneset_list = ['MSigDB_Hallmark_2020','KEGG_2021_Human','GO_Biological_Process_2023','Reactome_2022'] # 
deg_fdr_thresh = 0.1
deg_p_thresh = None

################################################
count1 = 0
for celltype in deg_res['celltype'].unique():
    print(celltype)
    # subset
    res_sub = deg_res[deg_res['celltype']==celltype]
    # filter or not based on stats of the DEtest
    if deg_fdr_thresh!=None:
        res_sub = res_sub[res_sub['FDR']<deg_fdr_thresh]
    elif deg_p_thresh!=None:
        res_sub = res_sub[res_sub['PValue']<deg_p_thresh]
    else:
        pass
    # GSEA
    if res_sub.shape[0]>=5:
        res_sub = res_sub[['gene_name','logFC']]
        res_sub.columns = [0,1]
        count2 = 0
        for i, geneset in enumerate(geneset_list):
            # use enrichr submodule to retrieve gene set library
            library = blitz.enrichr.get_library(geneset)
            # run enrichment analysis
            if __name__ == "__main__":  # make sure process is main, when run in a script it can cause errors otherwise
                result = blitz.gsea(res_sub, library)
                result = result.reset_index()
            if result.shape[0]>0:
                count2 += 1
                result['geneset'] = geneset
                # concatenate
                if count2==1:
                    result_concat = result.copy()
                else:
                    result_concat = pd.concat([result_concat,result])
            else:
                print(f'gsea, {geneset}: no results')
        if count2==0: # no results at all
            continue
    else:
        print(f'significant genes({list(res_sub["gene_name"])}) is less than 5')
        continue
    
    # concatenate
    count1 += 1
    result_concat['celltype'] = celltype
    if count1==1:
        result_concat_concat = result_concat.copy()
    else:
        result_concat_concat = pd.concat([result_concat_concat,result_concat])
    del result_concat

# save
result_concat_concat.to_csv('/nfs/team205/kk18/notebooks/Foetal/Trisomy21/RNA/milo_Mar2025_refine-nhoods/GSEA_trisomy21-vs-euploid-nhoods_finegrain.csv')
result_concat_concat
```

Example 8 (unknown):
```unknown
AtrialCardiomyocytesRight
AtrialCardiomyocytesLeft
CoronaryCapillaryEndothelialCells
EndocardialCells
SubEpicardialFibroblasts
GreatVesselAdventitialFibroblasts
GreatVesselSmoothMuscleCells
MacrophagesATF3pos
MacrophagesLYVE1pos
MyocardialInterstitialFibroblasts
VentricularCardiomyocytesCycling
VentricularCardiomyocytesLeftCompact
VentricularCardiomyocytesLeftTrabeculated
VentricularCardiomyocytesRightCompact
CPU times: user 21.6 s, sys: 1.09 s, total: 22.6 s
Wall time: 37.6 s
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Fig5_Trisomy21/no10_scoring_maturation-compaction-signatures.html

**Contents:**
- Function¶
- Read in adata¶
- Normalisation¶
- vCM Maturation signature genes¶

**Examples:**

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

import warnings
warnings.simplefilter("ignore", UserWarning)
warnings.filterwarnings("ignore", category=DeprecationWarning)
```

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

Example 3 (python):
```python
import importlib.util
import sys
spec = importlib.util.spec_from_file_location("module.name", "/rfs/project/rfs-iCNyzSAaucw/kk837/function/python/utils.py")
utils = importlib.util.module_from_spec(spec)
sys.modules["module.name"] = utils
spec.loader.exec_module(utils)
```

Example 4 (unknown):
```unknown
figdir = '/rfs/project/rfs-iCNyzSAaucw/kk837/notebooks/Foetal/fetal_ms_figures'
```

Example 5 (python):
```python
from scipy.stats import ranksums
from statsmodels.stats.multitest import multipletests

def plot_boxplot(adata, 
                 celltype,
                 score_name,
                 logFC_direction = 'euploid', # euploid or tri21 or both
                 fdr_thresh=0.1,
                 palette = ['gray','lightgray'],
                 savepath = None
                ):
    
    # subset adata
    ad = adata[adata.obs['nhood_annotation']==celltype]
    # target nhoods
    mask_sig = ad.obs['SpatialFDR']<fdr_thresh
    mask_eup = ad.obs['logFC']<0
    mask_t21 = ad.obs['logFC']>0
    nhoods_eup = list(ad.obs[mask_sig&mask_eup].index)
    nhoods_t21 = list(ad.obs[mask_sig&mask_t21].index)
    # add nhood labels
    ad.obs['group_to_compare'] = 'reference'
    ad.obs.loc[nhoods_eup,'group_to_compare'] = 'euploid_nhoods'
    ad.obs.loc[nhoods_t21,'group_to_compare'] = 'tri21_nhoods'
    # subset adata
    if logFC_direction=='both':
        ad.obs['group_to_compare'] = ad.obs['group_to_compare'].astype('category')
        ad.obs['group_to_compare'] = ad.obs['group_to_compare'].cat.reorder_categories(
            ['euploid_nhoods','reference','tri21_nhoods'])
    elif logFC_direction in ['tri21','euploid']:
        ad = ad[ad.obs['group_to_compare'].isin(['reference',f'{logFC_direction}_nhoods'])]
        ad.obs['group_to_compare'] = ad.obs['group_to_compare'].astype('category')
        ad.obs['group_to_compare'] = ad.obs['group_to_compare'].cat.reorder_categories(
            [f'{logFC_direction}_nhoods','reference'])
    else:
        return print('ERROR: logFC_direction is not appropriate')
    
    # prepare datafrane for plotting
    data = ad.obs.copy()
    
    ### plot ###
    split_col = 'group_to_compare'
    # Create a catplot with boxplots for each gene
    ncols=1
    g = sns.catplot(x=split_col, y=score_name, col='nhood_annotation', data=data, kind='box', palette=palette,
                    flierprops=dict(marker=''),
                    col_wrap=ncols, sharey=False, sharex=False, height=5, aspect=0.6)
    # Overlay swarmplot on each Axes
    # We'll use the 'hue' parameter to match the subsets and turn off dodge to align the points with the boxes.
    for celltype, ax in zip(data['nhood_annotation'].unique(), g.axes.flat):
        subset = data[data['nhood_annotation'] == celltype]
        sns.swarmplot(x=split_col, y=score_name, data=subset, color='k', size=2, ax=ax, dodge=False)

    # Update titles to remove "gene ="
    for ax, celltype in zip(g.axes.flat, data['nhood_annotation'].unique()):
        middle_ind = round(len(celltype)/2)
        celltype_title = f'{celltype[:middle_ind]}\n{celltype[middle_ind:]}'
        ax.set_title(celltype_title)

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

    # Adjust the spacing of the plots
    g.fig.subplots_adjust(wspace=0.4, hspace=0.4)
    
    # Show the plot
    if savepath!=None:
        plt.savefig(savepath,bbox_inches="tight",dpi=300)
    plt.ylabel(score_name)
    plt.show()

    # statistical test
    if logFC_direction=='both':
        # Wilcoxon Rank-Sum (Mann-Whitney U) Test
        # euploid vs reference
        ranksum_eup_vs_ref = ranksums(data[data['group_to_compare']=='reference'][score_name],
                                      data[data['group_to_compare']=='euploid_nhoods'][score_name])
        # trisomy21 vs reference
        ranksum_t21_vs_ref = ranksums(data[data['group_to_compare']=='reference'][score_name],
                                      data[data['group_to_compare']=='tri21_nhoods'][score_name])
        # trisomy21 vs euploid
        ranksum_t21_vs_eup = ranksums(data[data['group_to_compare']=='euploid_nhoods'][score_name],
                                      data[data['group_to_compare']=='tri21_nhoods'][score_name])
        # multipletest correction
        pvals_adjusted = multipletests([ranksum_eup_vs_ref[1],
                                        ranksum_t21_vs_ref[1],
                                        ranksum_t21_vs_eup[1]],alpha=0.05, method='fdr_bh')[1]
        print(f'adjusted p-values of eup-vs-ref, t21-vs-ref, t21-vs-eup: {pvals_adjusted}')
        
    elif logFC_direction in ['tri21','euploid']:
        print(ranksums(data[data['group_to_compare']=='reference'][score_name],
             data[data['group_to_compare']==f'{logFC_direction}_nhoods'][score_name]))
```

Example 6 (python):
```python
# pseudobulked per nhood
# nhood-by-gene
adata = sc.read_h5ad('/rfs/project/rfs-iCNyzSAaucw/kk837/data_objects/Foetal/trisomy21/Euploid_T21Hearts_Aug2024_sel_subsampled-per-cellstate-donor_pseudobulked-nhood-sum.h5ad')
adata
```

Example 7 (unknown):
```unknown
AnnData object with n_obs × n_vars = 6693 × 36601
    obs: 'index_cell', 'kth_distance', 'logFC', 'logCPM', 'F', 'PValue', 'FDR', 'SpatialFDR', 'Nhood_size', 'nhood_annotation', 'nhood_annotation_frac'
    var: 'gene_id', 'gene_name_scRNA-0-0', 'feature_type-1', 'mt-1', 'ribo-1', 'n_cells_by_counts-1', 'mean_counts-1', 'log1p_mean_counts-1', 'pct_dropout_by_counts-1', 'total_counts-1', 'log1p_total_counts-1', 'gene_name_GenomicFeatures', 'gene_biotype_GenomicFeatures', 'gene_seq_start_GenomicFeatures', 'gene_seq_end_GenomicFeatures', 'seq_name_GenomicFeatures', 'seq_strand_GenomicFeatures', 'seq_coord_system_GenomicFeatures', 'description_GenomicFeatures', 'gene_id_version_GenomicFeatures', 'canonical_transcript_GenomicFeatures', 'gene_name'
```

Example 8 (python):
```python
sc.pp.filter_genes(adata, min_cells=3)
# log1pCPM
sc.pp.normalize_total(adata, target_sum=1e6)
sc.pp.log1p(adata)
```

Example 9 (python):
```python
# read in gene signature
path_to_gene_signatures = '/rfs/project/rfs-iCNyzSAaucw/kk837/notebooks/Foetal/vCMLeft/gene_signatures.pkl'
target_signatures = ['compact_specific_up','common_up']

gene_signatures = utils.read_pkl(path_to_gene_signatures)
print(gene_signatures.keys())

sig_genes = []
for sig in target_signatures:
    print(sig)
    sig_genes = sig_genes + gene_signatures[sig]
print(len(sig_genes))
sig_genes[:5]
```

Example 10 (unknown):
```unknown
dict_keys(['compact_specific_up', 'trab_specific_up', 'common_up', 'compact_specific_down', 'trab_specific_down', 'common_down'])
compact_specific_up
common_up
466
```

Example 11 (unknown):
```unknown
['KIDINS220', 'ANKRD18A', 'MT-ATP6', 'ZNF385D', 'VEGFA']
```

Example 12 (python):
```python
# set some variables
score_name = 'vCMcompact_maturation_score'
celltypes_to_analyse = ['VentricularCardiomyocytesLeftCompact',
                         'VentricularCardiomyocytesLeftTrabeculated',
                         'VentricularCardiomyocytesRightCompact',
                         'VentricularCardiomyocytesRightTrabeculated',
                        'VentricularCardiomyocytesCycling']

if celltypes_to_analyse!=None:
    print('subsetting')
    adata_sub = adata[adata.obs['nhood_annotation'].isin(celltypes_to_analyse)]
else:
    print('just removing "Mixed"')
    adata_sub = adata[adata.obs['nhood_annotation']!='Mixed']
```

Example 13 (python):
```python
# obtain gene score
sc.tl.score_genes(adata_sub, gene_list=sig_genes, 
                  score_name=score_name, random_state=0)
```

Example 14 (python):
```python
# boxplot
plot_boxplot(adata_sub, 
                 celltype='VentricularCardiomyocytesLeftCompact',
                 score_name=score_name,
                 logFC_direction = 'both', # euploid or tri21
                 fdr_thresh=0.1,
                 palette = ['lightblue','lightgray','lightpink'],
                 savepath = f'{figdir}/Tri21_Milo-nhoods_Boxplot_vCMleftcompact_vCM-compaction-maturation.pdf'
                )
```

Example 15 (python):
```python
/tmp/ipykernel_2045939/1831124056.py:45: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `x` variable to `hue` and set `legend=False` for the same effect.

  g = sns.catplot(x=split_col, y=score_name, col='nhood_annotation', data=data, kind='box', palette=palette,
```

Example 16 (unknown):
```unknown
adjusted p-values of eup-vs-ref, t21-vs-ref, t21-vs-eup: [2.24125208e-15 5.07696046e-03 2.04802511e-08]
```

Example 17 (python):
```python
# boxplot
plot_boxplot(adata_sub, 
                 celltype='VentricularCardiomyocytesRightCompact',
                  score_name=score_name,
                 logFC_direction = 'both', # euploid or tri21
                 fdr_thresh=0.1,
                 palette = ['lightblue','lightgray','lightpink'],
                 savepath = f'{figdir}/Tri21_Milo-nhoods_Boxplot_vCMrightcompact_vCM-compaction-maturation.pdf'
                )
```

Example 18 (python):
```python
/tmp/ipykernel_2045939/1831124056.py:45: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `x` variable to `hue` and set `legend=False` for the same effect.

  g = sns.catplot(x=split_col, y=score_name, col='nhood_annotation', data=data, kind='box', palette=palette,
```

Example 19 (unknown):
```unknown
adjusted p-values of eup-vs-ref, t21-vs-ref, t21-vs-eup: [8.18327508e-09 4.29946937e-03 3.74536337e-03]
```

Example 20 (unknown):
```unknown
3.74536337e-03
```

Example 21 (unknown):
```unknown
0.00374536337
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Fig5_Trisomy21/no08_reassign-nhoods_significant-nhoods.html

**Contents:**
- Read in adata¶
- Identify significant nhoods¶
- Subset cells and nhoods¶
- Assess cells which belong to multipe nhoods¶
- Re-assign cells to a neighbourhood¶
- Save¶

--> most of the nhoods contain repetitive cells

**Examples:**

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

import warnings
warnings.simplefilter("ignore", UserWarning)
warnings.simplefilter(action='ignore', category=pd.errors.PerformanceWarning)
```

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

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

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

Example 5 (python):
```python
# read in post-Milo single-cell adata
adata = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/trisomy21/Euploid_T21Hearts_Aug2024_sel_subsampled-per-cellstate-donor.h5ad')
adata
```

Example 6 (python):
```python
AnnData object with n_obs × n_vars = 83477 × 36601
    obs: 'latent_RT_efficiency', 'latent_cell_probability', 'latent_scale', 'sangerID', 'combinedID', 'donor', 'region', 'age', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'donor_by_library-prep', 'multiplexed', 'SOC | status', 'SOC | log_prob_singleton', 'SOC | log_prob_doublet', 'batch_key', '_scvi_batch', 'FACSgate', 'fine_grain', 'mid_grain', 'coarse_grain', 'sex', 'week', 'trimester', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', '_scvi_labels', 'stress_score', 'hb1_score', 'euploid_or_t21', 'donor_type', 'region_finest', 'flushed', 'scrublet_leiden', 'cluster_scrublet_score', 'leiden_scVI_0.1', 'leiden_scVI_0.2', 'leiden_scVI_0.3', 'leiden_scVI_0.5', 'leiden_scVI_0.7', 'leiden_scVI_1', 'celltypist_coarse', 'conf_score_coarse', 'cells_to_remove', 'celltypist_fine', 'conf_score_fine', 'coarse_grain_merged', 'fine_grain_merged', 'nhood_ixs_random', 'nhood_ixs_refined', 'nhood_kth_distance'
    var: 'gene_id', 'gene_name_scRNA-0-0', 'feature_type-1', 'mt-1', 'ribo-1', 'n_cells_by_counts-1', 'mean_counts-1', 'log1p_mean_counts-1', 'pct_dropout_by_counts-1', 'total_counts-1', 'log1p_total_counts-1', 'gene_name_GenomicFeatures', 'gene_biotype_GenomicFeatures', 'gene_seq_start_GenomicFeatures', 'gene_seq_end_GenomicFeatures', 'seq_name_GenomicFeatures', 'seq_strand_GenomicFeatures', 'seq_coord_system_GenomicFeatures', 'description_GenomicFeatures', 'gene_id_version_GenomicFeatures', 'canonical_transcript_GenomicFeatures', 'gene_name'
    uns: 'coarse_grain_merged_colors', 'correcting-donor-region_n-layers-2', 'donor_colors', 'euploid_or_t21_colors', 'fine_grain_merged_colors', 'mid_grain_colors', 'nhood_adata', 'nhood_neighbors_key', 'region_colors', 'scVI_latent_correcting-donor-region_n-layers-2-forMilo', 'umap'
    obsm: 'X_umap_correcting-donor-region_n-layers-2', 'latent_gene_encoding', 'nhoods', 'nhoods_reassigned', 'scVI_latent_correcting-donor-region_n-layers-2'
    obsp: 'correcting-donor-region_n-layers-2_connectivities', 'correcting-donor-region_n-layers-2_distances', 'scVI_latent_correcting-donor-region_n-layers-2-forMilo_connectivities', 'scVI_latent_correcting-donor-region_n-layers-2-forMilo_distances'
```

Example 7 (python):
```python
# get significant nhoods
DA_fdr_thresh = 0.1
fcneg_label = 'euploid'
fcpos_label = 'trisomy21'

##############################
nhood_adata = adata.uns['nhood_adata'].copy()
mask_sig = nhood_adata.obs['SpatialFDR'] < DA_fdr_thresh
mask_fcneg = nhood_adata.obs['logFC'] < 0
mask_fcpos = nhood_adata.obs['logFC'] > 0
nhoods_fcneg = list(nhood_adata.obs_names[mask_sig&mask_fcneg])
print(f'FC-negative nhoods: {len(nhoods_fcneg)}')
print(nhoods_fcneg[:5])
nhoods_fcpos = list(nhood_adata.obs_names[mask_sig&mask_fcpos])
print(f'FC-positive nhoods: {len(nhoods_fcpos)}')
print(nhoods_fcpos[:5])

# add to nhood_adata
adata.uns['nhood_adata'].obs['deg_group'] = 'nan'
adata.uns['nhood_adata'].obs.loc[nhoods_fcneg,'deg_group'] = fcneg_label
adata.uns['nhood_adata'].obs.loc[nhoods_fcpos,'deg_group'] = fcpos_label
adata.uns['nhood_adata'].obs['deg_group'].value_counts()
```

Example 8 (unknown):
```unknown
FC-negative nhoods: 238
['nhood58', 'nhood157', 'nhood161', 'nhood178', 'nhood185']
FC-positive nhoods: 339
['nhood2', 'nhood20', 'nhood36', 'nhood53', 'nhood105']
```

Example 9 (unknown):
```unknown
deg_group
nan          6116
trisomy21     339
euploid       238
Name: count, dtype: int64
```

Example 10 (unknown):
```unknown
# gent nhood index
mhoods_sel = nhoods_fcneg + nhoods_fcpos
nhood_index = [int(x.replace('nhood','')) for x in mhoods_sel]
nhood_index.sort()
```

Example 11 (python):
```python
# get cell obsnames
# which belongs to the significant nhoods obtained above
mask = adata.obsm['nhoods'][:,nhood_index].toarray().sum(axis=1)>0 # get cells which has at least one neighbourhod to belong
cell_obsnames = adata.obs_names[mask]
print(len(cell_obsnames))
```

Example 12 (python):
```python
# subset
## cell
adata_sub = adata[cell_obsnames]
## nhoods, in the two different data storage
adata_sub.obsm['nhoods'] = adata_sub.obsm['nhoods'][:,nhood_index]
print(adata_sub.obsm['nhoods'].shape)
adata_sub.uns['nhood_adata'] = adata_sub.uns['nhood_adata'][nhood_index]
print(adata_sub.uns['nhood_adata'].shape)
```

Example 13 (unknown):
```unknown
(14846, 577)
(577, 24)
```

Example 14 (python):
```python
mask = adata_sub.obsm['nhoods'].toarray().sum(axis=1)>1 # get cells which has at least two neighbourhod to belong
cells_rep = adata_sub.obs_names[mask]
print(f'number of cells: {adata_sub.shape[0]}')
print(f'fraction of repetitive cells: {len(cells_rep)/adata_sub.shape[0]}')
```

Example 15 (unknown):
```unknown
number of cells: 14846
fraction of repetitive cells: 0.3483766671157214
```

Example 16 (unknown):
```unknown
len(cells_rep)
```

Example 17 (python):
```python
print(adata_sub.obsm['nhoods'].shape)
plt.hist(adata_sub.obsm['nhoods'].toarray().sum(axis=1),bins=10)
plt.xlabel('number of nhoods')
plt.show()
```

Example 18 (unknown):
```unknown
(14846, 577)
```

Example 19 (python):
```python
# obtain nhoods which contains repetitive cells
mask = adata_sub[cells_rep].obsm['nhoods'].toarray().sum(axis=0)>0
nhoods_rep = np.nonzero(mask)[0]
print(len(nhoods_rep))
nhoods_rep[:5]
```

Example 20 (unknown):
```unknown
array([0, 1, 2, 3, 4])
```

Example 21 (python):
```python
# get centroid for each nhood
cell_by_nhood = adata_sub.obsm['nhoods'][:,nhoods_rep].copy()
n_cell_per_nhood = cell_by_nhood.toarray().sum(axis=0)
scVI_latent_rep = 'scVI_latent_correcting-donor-region_n-layers-2'
cell_by_scvi = adata_sub.obsm[scVI_latent_rep].copy()
nhood_by_scvi_sum = np.array(np.matmul(cell_by_nhood.T.todense(),cell_by_scvi))
nhood_by_scvi_mean = (nhood_by_scvi_sum.T/n_cell_per_nhood).T # i.e. centroid
nhood_by_scvi_mean = pd.DataFrame(nhood_by_scvi_mean,
                                  index=[f'nhood{str(x)}' for x in nhoods_rep])
print(nhood_by_scvi_mean.shape)
nhood_by_scvi_mean.head()
```

Example 22 (python):
```python
# for each repetitive cell, calculate distance to the centroids of the nhoods, and select the most closest one
import numpy as np
def assign_cells_to_neighborhoods(cell_nhood, cell_latent, nhood_centroids):
    """
    Assigns each cell to the closest neighborhood (nhood) centroid among the neighborhoods 
    the cell belongs to.
    
    Parameters:
      cell_nhood: 2D binary NumPy array of shape (n_cells, n_nhoods).
                  cell_nhood[i, j] is 1 if cell i belongs to nhood j.
      cell_latent: 2D NumPy array of shape (n_cells, latent_dim) containing the cell coordinates.
      nhood_centroids: 2D NumPy array of shape (n_nhoods, latent_dim) containing the centroids.
    
    Returns:
      A 1D array of length n_cells where each element is the index of the nearest nhood centroid 
      (among those the cell belongs to) for that cell.
      
    Explanation
    1. Distance Computation: The squared Euclidean distance between each cell (vector in latent space) and each nhood centroid is computed using a vectorized approach. This avoids explicit loops and leverages fast linear algebra routines.
    2. Masking: Since each cell only belongs to certain nhoods (as indicated by the cell_nhood binary matrix), we set the distances for non-member nhoods to infinity (np.inf). This ensures that the minimization only considers the valid nhoods.
    3. Choosing the Minimum: Finally, we use np.argmin along the nhood axis (axis=1) to find the index of the closest centroid among those the cell belongs to.  
    """
    # Compute squared Euclidean distances between each cell and each nhood centroid.
    # Using the identity: ||a-b||^2 = ||a||^2 + ||b||^2 - 2 * a.dot(b)
    cell_sq = np.sum(cell_latent**2, axis=1, keepdims=True)   # shape: (n_cells, 1)
    nhood_sq = np.sum(nhood_centroids**2, axis=1)  # shape: (n_nhoods,)
    distances = cell_sq + nhood_sq - 2 * (cell_latent @ nhood_centroids.T)  # shape: (n_cells, n_nhoods)
    
    # Mask distances for nhoods that the cell does not belong to.
    # For these entries, set the distance to a very high value (np.inf) so they are not chosen.
    mask = cell_nhood.astype(bool)
    distances[~mask] = np.inf
    
    # For each cell, find the index of the nhood with the minimal distance.
    closest_nhood = np.argmin(distances, axis=1)
    
    return closest_nhood
```

Example 23 (python):
```python
%%time
cell_by_nhood = adata_sub[cells_rep].obsm['nhoods'][:,nhoods_rep].toarray().copy()
print(cell_by_nhood.shape)
cell_by_scVI = adata_sub[cells_rep].obsm[scVI_latent_rep].toarray().copy()
print(cell_by_scVI.shape)

assignments = assign_cells_to_neighborhoods(cell_by_nhood, cell_by_scVI, np.array(nhood_by_scvi_mean))
assignments = nhood_by_scvi_mean.index[assignments]
print(len(assignments))
assignments[:5]
```

Example 24 (unknown):
```unknown
(5172, 559)
(5172, 50)
5172
CPU times: user 98.8 ms, sys: 2.6 ms, total: 101 ms
Wall time: 104 ms
```

Example 25 (unknown):
```unknown
Index(['nhood566', 'nhood398', 'nhood395', 'nhood79', 'nhood6'], dtype='object')
```

Example 26 (python):
```python
# reassign nhoods
cell_index = np.nonzero(adata_sub.obs_names.isin(cells_rep))[0]
nhood_reassigned_index = [int(x.replace('nhood','')) for x in assignments]

cell_by_nhood = adata_sub.obsm['nhoods'].toarray().copy()
cell_by_nhood[cell_index,:] = 0
cell_by_nhood[cell_index,nhood_reassigned_index] = 1
from scipy.sparse import csr_matrix
adata_sub.obsm['nhoods_reassigned'] = csr_matrix(cell_by_nhood)
adata_sub.obsm['nhoods_reassigned']
```

Example 27 (python):
```python
<14846x577 sparse matrix of type '<class 'numpy.float32'>'
	with 14846 stored elements in Compressed Sparse Row format>
```

Example 28 (python):
```python
# check
print(np.max(adata_sub.obsm['nhoods_reassigned'].toarray().sum(axis=1)))
plt.hist(adata_sub.obsm['nhoods_reassigned'].toarray().sum(axis=1))
plt.xlabel('number of nhoods')
plt.show()
```

Example 29 (python):
```python
# n_cells per nhood
plt.hist(adata_sub.obsm['nhoods'].toarray().sum(axis=0),bins=100)
plt.title('original')
plt.show()

plt.hist(adata_sub.obsm['nhoods_reassigned'].toarray().sum(axis=0),bins=100)
plt.title('reassigned')
plt.show()
```

Example 30 (python):
```python
# save
adata_sub.write('/nfs/team205/heart/anndata_objects/Foetal/trisomy21/Euploid_T21Hearts_Aug2024_sel_subsampled-per-cellstate-donor_cells-in-sig-nhoods-sfdr0p1.h5ad')
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Fig5_Trisomy21/no03-2_fraction-of-compactCM.html

**Contents:**
- Import modules¶
- Read in anndata¶
- Calculate and compare fraction¶
- Left vCMs¶
- Right vCMs¶

• Calculating fraction using NOT subsampled adata

**Examples:**

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

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 numpy as np
import numpy.random as random
import pandas as pd
import scanpy as sc
```

Example 2 (python):
```python
sc.settings.verbosity = 3             # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.settings.set_figure_params(dpi=80, facecolor='white', color_map = 'RdPu',)
```

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

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

Example 5 (python):
```python
import os
os.getcwd()
```

Example 6 (unknown):
```unknown
'/nfs/team205/kk18/notebooks/Foetal/Trisomy21/RNA'
```

Example 7 (python):
```python
adata = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/trisomy21/Euploid_T21Hearts_Aug2024_sel.h5ad')
adata
```

Example 8 (unknown):
```unknown
AnnData object with n_obs × n_vars = 147494 × 36601
    obs: 'latent_RT_efficiency', 'latent_cell_probability', 'latent_scale', 'sangerID', 'combinedID', 'donor', 'region', 'age', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'donor_by_library-prep', 'multiplexed', 'SOC | status', 'SOC | log_prob_singleton', 'SOC | log_prob_doublet', 'batch_key', '_scvi_batch', 'FACSgate', 'fine_grain', 'mid_grain', 'coarse_grain', 'sex', 'week', 'trimester', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', '_scvi_labels', 'stress_score', 'hb1_score', 'euploid_or_t21', 'donor_type', 'region_finest', 'flushed', 'scrublet_leiden', 'cluster_scrublet_score', 'leiden_scVI_0.1', 'leiden_scVI_0.2', 'leiden_scVI_0.3', 'leiden_scVI_0.5', 'leiden_scVI_0.7', 'leiden_scVI_1', 'celltypist_coarse', 'conf_score_coarse', 'cells_to_remove', 'celltypist_fine', 'conf_score_fine', 'coarse_grain_merged', 'fine_grain_merged'
    var: 'gene_id', 'gene_name_scRNA-0-0', 'feature_type-1', 'ambient_expression-0-1', 'ambient_expression-1-1', 'ambient_expression-10-1', 'ambient_expression-11-1', 'ambient_expression-12-1', 'ambient_expression-13-1', 'ambient_expression-14-1', 'ambient_expression-2-1', 'ambient_expression-3-1', 'ambient_expression-4-1', 'ambient_expression-5-1', 'ambient_expression-6-1', 'ambient_expression-7-1', 'ambient_expression-8-1', 'ambient_expression-9-1', 'mt-1', 'ribo-1', 'n_cells_by_counts-1', 'mean_counts-1', 'log1p_mean_counts-1', 'pct_dropout_by_counts-1', 'total_counts-1', 'log1p_total_counts-1'
    obsm: 'X_scVI', 'X_umap', 'latent_gene_encoding'
```

Example 9 (python):
```python
from scipy.stats import ranksums
```

Example 10 (python):
```python
cells_to_include = ['VentricularCardiomyocytesLeftCompact','VentricularCardiomyocytesLeftTrabeculated']
group_col = 'fine_grain_merged'

ad = adata[adata.obs[group_col].isin(cells_to_include)]
ad.obs['condition_donor'] = ad.obs['euploid_or_t21'].astype('str')+'_'+ad.obs['donor'].astype('str')

# sangerID proportion
crosstab=pd.crosstab(ad.obs['condition_donor'],ad.obs[group_col])
# Normalize each row to sum to 1 (i.e., to get fractions)
crosstab_normalized = crosstab.div(crosstab.sum(axis=1), axis=0)

# Plotting the stacked bar chart
ax = crosstab_normalized.plot(kind='bar',stacked=True,figsize=(4, 4),legend=False)

# Adding labels and title
plt.rcParams['axes.grid'] = False
plt.ylabel('Fraction')
plt.title('Normalized Stacked Bar Chart')
plt.legend(bbox_to_anchor=(0.7, -1, 0.5, 0.5))
plt.xlabel(None)
# Show the plot
plt.show()
```

Example 11 (python):
```python
cells_to_include = ['VentricularCardiomyocytesLeftCompact','VentricularCardiomyocytesLeftTrabeculated']
group_col = 'fine_grain_merged'
split_col = 'euploid_or_t21'
celltype = 'VentricularCardiomyocytesLeftCompact'
ncols = 6


ad = adata[adata.obs[group_col].isin(cells_to_include)]
ad.obs['condition_donor'] = ad.obs['euploid_or_t21'].astype('str')+'_'+ad.obs['donor'].astype('str')

# sangerID proportion
crosstab=pd.crosstab(ad.obs['condition_donor'],ad.obs[group_col])
# Normalize each row to sum to 1 (i.e., to get fractions)
crosstab_normalized = crosstab.div(crosstab.sum(axis=1), axis=0)

# prepare data for plotting
data = crosstab_normalized.reset_index().melt(id_vars=['condition_donor'])
data['euploid_or_t21'] = [x.split('_')[0] for x in data['condition_donor']]
data = data[data[group_col]==celltype]
palette = ['gray','lightgray']

# Create a catplot with boxplots for each gene
g = sns.catplot(x=split_col, y='value', col='fine_grain_merged', data=data, kind='box', palette=palette,
                col_wrap=ncols, sharey=False, sharex=False, height=3, aspect=0.6)

# Overlay swarmplot on each Axes
# We'll use the 'hue' parameter to match the subsets and turn off dodge to align the points with the boxes.
for celltype, ax in zip(data['fine_grain_merged'].unique(), g.axes.flat):
    subset = data[data['fine_grain_merged'] == celltype]
    sns.swarmplot(x=split_col, y='value', data=subset, color='k', size=5, ax=ax, dodge=False)

# Update titles to remove "gene ="
for ax, celltype in zip(g.axes.flat, data['fine_grain_merged'].unique()):
    middle_ind = round(len(celltype)/2)
    celltype_title = f'{celltype[:middle_ind]}\n{celltype[middle_ind:]}'
    ax.set_title(celltype_title)

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

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

# Manually create a legend for the groups.
# Since there are no legend handles for the swarmplot, we create them manually.
from matplotlib.patches import Patch
legend_patches = [Patch(facecolor=palette[0], edgecolor='k', label='euploid'),
                  Patch(facecolor=palette[1], edgecolor='k', label='tri21'),
                  ]

# Add a single legend for the entire figure with our custom patches
g.fig.legend(handles=legend_patches, labels=['euploid', 'tri21'], title='', loc='center', bbox_to_anchor=(0.3, 0.8))

# Show the plot
plt.savefig(f'{figdir}/Tri21_BoxPlot_vCMLeft-fraction-compact-trab.pdf',bbox_inches="tight",dpi=300)
plt.show()

# statistical test
ranksums(data[data['euploid_or_t21']=='euploid']['value'],
         data[data['euploid_or_t21']=='tri21']['value']
        )
```

Example 12 (python):
```python
/tmp/ipykernel_3787554/1312438397.py:24: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `x` variable to `hue` and set `legend=False` for the same effect.

  g = sns.catplot(x=split_col, y='value', col='fine_grain_merged', data=data, kind='box', palette=palette,
```

Example 13 (unknown):
```unknown
RanksumsResult(statistic=2.1213203435596424, pvalue=0.03389485352468927)
```

Example 14 (python):
```python
cells_to_include = ['VentricularCardiomyocytesRightCompact','VentricularCardiomyocytesRightTrabeculated']
group_col = 'fine_grain_merged'

ad = adata[adata.obs[group_col].isin(cells_to_include)]
ad.obs['condition_donor'] = ad.obs['euploid_or_t21'].astype('str')+'_'+ad.obs['donor'].astype('str')

# sangerID proportion
crosstab=pd.crosstab(ad.obs['condition_donor'],ad.obs[group_col])
# Normalize each row to sum to 1 (i.e., to get fractions)
crosstab_normalized = crosstab.div(crosstab.sum(axis=1), axis=0)

# Plotting the stacked bar chart
ax = crosstab_normalized.plot(kind='bar',stacked=True,figsize=(4, 4),legend=False)

# Adding labels and title
plt.rcParams['axes.grid'] = False
plt.ylabel('Fraction')
plt.title('Normalized Stacked Bar Chart')
plt.legend(bbox_to_anchor=(0.7, -1, 0.5, 0.5))
plt.xlabel(None)
# Show the plot
plt.show()
```

Example 15 (python):
```python
cells_to_include = ['VentricularCardiomyocytesRightCompact','VentricularCardiomyocytesRightTrabeculated']
group_col = 'fine_grain_merged'
split_col = 'euploid_or_t21'
celltype = 'VentricularCardiomyocytesRightCompact'


ad = adata[adata.obs[group_col].isin(cells_to_include)]
ad.obs['condition_donor'] = ad.obs['euploid_or_t21'].astype('str')+'_'+ad.obs['donor'].astype('str')

# sangerID proportion
crosstab=pd.crosstab(ad.obs['condition_donor'],ad.obs[group_col])
# Normalize each row to sum to 1 (i.e., to get fractions)
crosstab_normalized = crosstab.div(crosstab.sum(axis=1), axis=0)

# prepare data for plotting
data = crosstab_normalized.reset_index().melt(id_vars=['condition_donor'])
data['euploid_or_t21'] = [x.split('_')[0] for x in data['condition_donor']]
data = data[data[group_col]==celltype]
palette = ['gray','lightgray']

# Create a catplot with boxplots for each gene
g = sns.catplot(x=split_col, y='value', col='fine_grain_merged', data=data, kind='box', palette=palette,
                col_wrap=ncols, sharey=False, sharex=False, height=3, aspect=0.6)

# Overlay swarmplot on each Axes
# We'll use the 'hue' parameter to match the subsets and turn off dodge to align the points with the boxes.
for celltype, ax in zip(data['fine_grain_merged'].unique(), g.axes.flat):
    subset = data[data['fine_grain_merged'] == celltype]
    sns.swarmplot(x=split_col, y='value', data=subset, color='k', size=5, ax=ax, dodge=False)

# Update titles to remove "gene ="
for ax, celltype in zip(g.axes.flat, data['fine_grain_merged'].unique()):
    middle_ind = round(len(celltype)/2)
    celltype_title = f'{celltype[:middle_ind]}\n{celltype[middle_ind:]}'
    ax.set_title(celltype_title)

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

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

# Manually create a legend for the groups.
# Since there are no legend handles for the swarmplot, we create them manually.
from matplotlib.patches import Patch
legend_patches = [Patch(facecolor=palette[0], edgecolor='k', label='euploid'),
                  Patch(facecolor=palette[1], edgecolor='k', label='tri21'),
                  ]

# Add a single legend for the entire figure with our custom patches
g.fig.legend(handles=legend_patches, labels=['euploid', 'tri21'], title='', loc='center', bbox_to_anchor=(0.3, 0.8))

# Show the plot
plt.savefig(f'{figdir}/Tri21_BoxPlot_vCMRight-fraction-compact-trab.pdf',bbox_inches="tight",dpi=300)
plt.show()

# statistical test
ranksums(data[data['euploid_or_t21']=='euploid']['value'],
         data[data['euploid_or_t21']=='tri21']['value']
        )
```

Example 16 (python):
```python
/tmp/ipykernel_3787554/2915419533.py:22: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `x` variable to `hue` and set `legend=False` for the same effect.

  g = sns.catplot(x=split_col, y='value', col='fine_grain_merged', data=data, kind='box', palette=palette,
```

Example 17 (unknown):
```unknown
RanksumsResult(statistic=2.1213203435596424, pvalue=0.03389485352468927)
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Fig5_Trisomy21/no07_milo.html

**Contents:**
- Read in adata¶
- Milo¶
- scVI - n_layers:2¶
- Plot¶

6693 rows × 8 columns

--> logFC<0 is euploid

577 rows × 11 columns

**Examples:**

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

import warnings
warnings.simplefilter("ignore", UserWarning)
```

Example 2 (python):
```python
'''
# Configure rpy2 to Use the Correct R Installation:
import os
os.environ['R_HOME'] = '/nfs/team205/kk18/miniconda3/envs/vitro/lib/R'
# check
from rpy2.rinterface_lib import openrlib
print(openrlib.R_HOME)
'''
```

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

Example 4 (python):
```python
# import session_info
# session_info.show()
```

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

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

Example 7 (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
# atlas
adata = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/trisomy21/Euploid_T21Hearts_Aug2024_sel_subsampled-per-cellstate-donor.h5ad')
adata
```

Example 9 (python):
```python
AnnData object with n_obs × n_vars = 83477 × 36601
    obs: 'latent_RT_efficiency', 'latent_cell_probability', 'latent_scale', 'sangerID', 'combinedID', 'donor', 'region', 'age', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'donor_by_library-prep', 'multiplexed', 'SOC | status', 'SOC | log_prob_singleton', 'SOC | log_prob_doublet', 'batch_key', '_scvi_batch', 'FACSgate', 'fine_grain', 'mid_grain', 'coarse_grain', 'sex', 'week', 'trimester', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', '_scvi_labels', 'stress_score', 'hb1_score', 'euploid_or_t21', 'donor_type', 'region_finest', 'flushed', 'scrublet_leiden', 'cluster_scrublet_score', 'leiden_scVI_0.1', 'leiden_scVI_0.2', 'leiden_scVI_0.3', 'leiden_scVI_0.5', 'leiden_scVI_0.7', 'leiden_scVI_1', 'celltypist_coarse', 'conf_score_coarse', 'cells_to_remove', 'celltypist_fine', 'conf_score_fine', 'coarse_grain_merged', 'fine_grain_merged', 'nhood_ixs_random', 'nhood_ixs_refined', 'nhood_kth_distance'
    var: 'gene_id', 'gene_name_scRNA-0-0', 'feature_type-1', 'mt-1', 'ribo-1', 'n_cells_by_counts-1', 'mean_counts-1', 'log1p_mean_counts-1', 'pct_dropout_by_counts-1', 'total_counts-1', 'log1p_total_counts-1', 'gene_name_GenomicFeatures', 'gene_biotype_GenomicFeatures', 'gene_seq_start_GenomicFeatures', 'gene_seq_end_GenomicFeatures', 'seq_name_GenomicFeatures', 'seq_strand_GenomicFeatures', 'seq_coord_system_GenomicFeatures', 'description_GenomicFeatures', 'gene_id_version_GenomicFeatures', 'canonical_transcript_GenomicFeatures', 'gene_name'
    uns: 'coarse_grain_merged_colors', 'correcting-donor-region_n-layers-2', 'donor_colors', 'euploid_or_t21_colors', 'fine_grain_merged_colors', 'mid_grain_colors', 'nhood_adata', 'nhood_neighbors_key', 'region_colors', 'scVI_latent_correcting-donor-region_n-layers-2-forMilo', 'umap'
    obsm: 'X_umap_correcting-donor-region_n-layers-2', 'latent_gene_encoding', 'nhoods', 'scVI_latent_correcting-donor-region_n-layers-2'
    obsp: 'correcting-donor-region_n-layers-2_connectivities', 'correcting-donor-region_n-layers-2_distances', 'scVI_latent_correcting-donor-region_n-layers-2-forMilo_connectivities', 'scVI_latent_correcting-donor-region_n-layers-2-forMilo_distances'
```

Example 10 (unknown):
```unknown
basis = 'umap_correcting-donor-region_n-layers-2'
```

Example 11 (python):
```python
pd.crosstab(adata.obs['donor'],adata.obs['sangerID']).T
```

Example 12 (python):
```python
sc.pl.embedding(adata, basis=basis,
                    color=['coarse_grain_merged','euploid_or_t21'],
               size=5, wspace=0.5, ncols=6, components='all',cmap='RdPu',frameon=False,legend_fontsize=14)
```

Example 13 (python):
```python
# umap per celltype to highlight
group_col = 'region'
groups_to_plot = ['Heart',
                  'Great vessels', 
                 'Left atria',
                 'Left ventricle',
                 'Pericardial sac + RA + LA + IAS',
                 'Right ventricle',
                  # 'Pericardium',
                 ]

plt.rcParams['axes.grid'] = False
fig,axs=plt.subplots(nrows=2,ncols=4,figsize=(15,7))
# Flatten the axs array for easier iteration in case of multiple rows and columns
axs_flat = axs.flatten()

for i,grp in enumerate(groups_to_plot):
    sc.pl.embedding(adata, basis=basis,
                    color=[group_col], groups=grp, 
               size=10, frameon=False,
               show=False, ax=axs_flat[i], legend_loc=None, legend_fontsize=5,
               title=grp
              )
# Turn off any remaining subplots
for j in range(len(groups_to_plot), len(axs_flat)):
    fig.delaxes(axs_flat[j])

# fig.suptitle('Main Title for the Figure')
fig.tight_layout()
```

Example 14 (python):
```python
sc.pl.embedding(adata[adata.obs['euploid_or_t21']=='tri21'], 
                basis=basis,
                    color=['conf_score_fine'], 
               size=10, wspace=0.4, ncols=6, components='all',frameon=False,legend_fontsize=30, title='tri21, celltypist confidence score')
```

Example 15 (python):
```python
adata.obs['euploid_or_t21'].cat.categories
```

Example 16 (unknown):
```unknown
Index(['euploid', 'tri21'], dtype='object')
```

Example 17 (python):
```python
%%time
scVI_latent_rep = 'scVI_latent_correcting-donor-region_n-layers-2'
# compute the k-nearest-neighbor graph that is used in both clustering and umap algorithms
sc.pp.neighbors(adata, use_rep=scVI_latent_rep, n_neighbors=20, key_added=f"{scVI_latent_rep}-forMilo")
# Construct neighbourhoods
milo.make_nhoods(adata, neighbors_key=f"{scVI_latent_rep}-forMilo", prop=0.1)
# Count cells in neighbourhoods
# Milo leverages the variation in cell numbers between replicates for the same experimental condition to test for differential abundance. 
milo.count_nhoods(adata, sample_col="sangerID")
# Test for differential abundance between conditions
milo.DA_nhoods(adata, design="~ euploid_or_t21")
## Check results
milo_results = adata.uns["nhood_adata"].obs
milo_results
```

Example 18 (unknown):
```unknown
CPU times: user 24.7 s, sys: 423 ms, total: 25.1 s
Wall time: 26.2 s
```

Example 19 (unknown):
```unknown
milo_results['SpatialFDR'].hist(bins=100)
```

Example 20 (python):
```python
nhood_size = np.array(adata.obsm["nhoods"].sum(0)).ravel()
plt.hist(nhood_size, bins=100);
```

Example 21 (python):
```python
sc.pl.embedding(adata, basis=basis,
                    color=['coarse_grain_merged'],
               size=5, wspace=0.5, ncols=6, components='all',cmap='RdPu',frameon=False,legend_fontsize=14)
```

Example 22 (python):
```python
%%time
# Visualize results on UMAP embedding
milopy.utils.build_nhood_graph(adata,basis=f'X_{basis}')
plt.rcParams["figure.figsize"] = [12,10]
milopl.plot_nhood_graph(adata, 
                        alpha=0.1, ## SpatialFDR level
                        min_size=2, ## Size of smallest dot
                        show=False
                       )
plt.savefig(f'{figdir}/Tri21_RNA_Milo_UMAP.pdf', bbox_inches="tight",dpi=300)
```

Example 23 (unknown):
```unknown
CPU times: user 792 ms, sys: 8.5 ms, total: 801 ms
Wall time: 1.22 s
```

Example 24 (python):
```python
milopy.utils.annotate_nhoods(adata, anno_col='fine_grain_merged')
```

Example 25 (python):
```python
adata_copy = adata.copy()
milopy.utils.annotate_nhoods(adata_copy, anno_col='euploid_or_t21')
df = adata_copy.uns['nhood_adata'].obs.copy()
df[df['SpatialFDR']<0.1]
```

Example 26 (python):
```python
plt.hist(adata.uns['nhood_adata'].obs["nhood_annotation_frac"],bins=100);
plt.xlabel("celltype fraction")
```

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

Example 28 (python):
```python
# We can see that for the majority of neighbourhoods, almost all cells have the same neighbourhood. We can rename neighbourhoods where less than 60% of the cells have the top label as "Mixed"
adata.uns['nhood_adata'].obs.loc[adata.uns['nhood_adata'].obs["nhood_annotation_frac"] < 0.6, "nhood_annotation"] = "Mixed"
```

Example 29 (python):
```python
# save
# modify nhood names
adata.uns['nhood_adata'].obs_names = 'nhood'+adata.uns['nhood_adata'].obs_names
adata.write('/nfs/team205/heart/anndata_objects/Foetal/trisomy21/Euploid_T21Hearts_Aug2024_sel_subsampled-per-cellstate-donor.h5ad')
adata.uns['nhood_adata'].obs.head()
```

Example 30 (python):
```python
adata = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/trisomy21/Euploid_T21Hearts_Aug2024_sel_subsampled-per-cellstate-donor.h5ad')
adata.uns['nhood_adata'].obs.head()
```

Example 31 (python):
```python
# from CamU
adata = sc.read_h5ad('/rfs/project/rfs-iCNyzSAaucw/kk837/data_objects/Foetal/trisomy21/Euploid_T21Hearts_Aug2024_sel_subsampled-per-cellstate-donor.h5ad')
adata.uns['nhood_adata'].obs.head()
```

Example 32 (unknown):
```unknown
figdir = '/rfs/project/rfs-iCNyzSAaucw/kk837/notebooks/Foetal/fetal_ms_figures'
```

Example 33 (python):
```python
# map cell type name to short versions
df = pd.read_csv('/rfs/project/rfs-iCNyzSAaucw/kk837/notebooks/Foetal/finegrain_name_mapping.csv').set_index('Full_name')
mapping = df['Short_name_finalised'].to_dict()
adata.uns['nhood_adata'].obs.replace({'nhood_annotation':mapping},inplace=True)
```

Example 34 (python):
```python
# drop cell types which have lower than certain number of significant nhoods
sig_nhoods_thresh = 3

ad = adata.uns['nhood_adata'].copy()

data = ad.obs[['SpatialFDR','nhood_annotation','logFC']]
data = data[data['nhood_annotation']!='Mixed']
# drop cell types which has lower than 3 nhood
val_counts = data['nhood_annotation'].value_counts()
data = data[data['nhood_annotation'].isin(list(val_counts.index[val_counts>=3]))]
# highlight significant nhoods
mask_sig = data['SpatialFDR'] < 0.1
mask_eup = data['logFC'] < 0
mask_t21 = data['logFC'] > 0
nhoods_eup = list(data.index[mask_sig&mask_eup])
nhoods_t21 = list(data.index[mask_sig&mask_t21])
data['nhood_type'] = 'mixture'
data.loc[nhoods_eup,'nhood_type'] = 'euploid'
data.loc[nhoods_t21,'nhood_type'] = 'tri21'

# select cell types based on the number of significant nhoods
data['significant'] = data['SpatialFDR'] < 0.1 # to count number of significant nhood
sig_counts = data[['nhood_annotation','significant']].groupby('nhood_annotation').sum()
mask = sig_counts['significant'] >= sig_nhoods_thresh
cells_sel = sig_counts.index[mask]
data = data[data['nhood_annotation'].isin(cells_sel)]

# reorder cell types based on the number of significant nhoods
data['nhood_annotation'] = data['nhood_annotation'].astype('category')
data['nhood_annotation'] = data['nhood_annotation'].cat.reorder_categories(
    list(data.groupby('nhood_annotation').sum().sort_values('significant',ascending=False).index)
)

# Set the size of the plot
plt.figure(figsize=(3, 6))

# Create the violin plot
sns.violinplot(data=data, x='logFC', y='nhood_annotation', inner='box', scale='width', color='lightgray')

# Add jittered points without dodge to align them on the same y-axis location
sns.stripplot(data=data, x='logFC', y='nhood_annotation', jitter=True, size=5,
              hue='nhood_type', palette={'euploid': 'blue', 
                                         'mixture': 'gray',
                                         'tri21': 'red'
                                        }, 
              marker='o', alpha=0.7)

# Remove the legend created by hue
plt.legend([],[], frameon=False)
# Display the plot
# plt.title(comp)
plt.savefig(f'{figdir}/Tri21_RNA_Milo_beeswarm_sort_by_n-significant_celltype-selected.pdf', bbox_inches="tight",dpi=300)
plt.show()
```

Example 35 (python):
```python
/tmp/ipykernel_293470/1713106557.py:31: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.
  list(data.groupby('nhood_annotation').sum().sort_values('significant',ascending=False).index)
/tmp/ipykernel_293470/1713106557.py:38: FutureWarning: 

The `scale` parameter has been renamed and will be removed in v0.15.0. Pass `density_norm='width'` for the same effect.
  sns.violinplot(data=data, x='logFC', y='nhood_annotation', inner='box', scale='width', color='lightgray')
```

Example 36 (python):
```python
# all the cell types

ad = adata.uns['nhood_adata'].copy()

data = ad.obs[['SpatialFDR','nhood_annotation','logFC']]
data = data[data['nhood_annotation']!='Mixed']
# drop cell types which has lower than 3 nhood
val_counts = data['nhood_annotation'].value_counts()
data = data[data['nhood_annotation'].isin(list(val_counts.index[val_counts>=3]))]
# highlight significant nhoods
mask_sig = data['SpatialFDR'] < 0.1
mask_eup = data['logFC'] < 0
mask_t21 = data['logFC'] > 0
nhoods_eup = list(data.index[mask_sig&mask_eup])
nhoods_t21 = list(data.index[mask_sig&mask_t21])
data['nhood_type'] = 'mixture'
data.loc[nhoods_eup,'nhood_type'] = 'euploid'
data.loc[nhoods_t21,'nhood_type'] = 'tri21'

# reorder based on number of significant nhoods
# reorder cell types
data['nhood_annotation'] = data['nhood_annotation'].astype('category')
data['significant'] = data['SpatialFDR'] < 0.1 # to count number of significant nhood
data['nhood_annotation'] = data['nhood_annotation'].cat.reorder_categories(
    list(data.groupby('nhood_annotation').sum().sort_values('significant',ascending=False).index)
)

# Set the size of the plot
plt.figure(figsize=(3, 10))

# Create the violin plot
sns.violinplot(data=data, x='logFC', y='nhood_annotation', inner='box', scale='width', color='lightgray')

# Add jittered points without dodge to align them on the same y-axis location
sns.stripplot(data=data, x='logFC', y='nhood_annotation', jitter=True, size=7,
              hue='nhood_type', palette={'euploid': 'blue', 
                                         'mixture': 'gray',
                                         'tri21': 'red'
                                        }, 
              marker='o', alpha=0.7)

# Remove the legend created by hue
plt.legend([],[], frameon=False)
# Display the plot
# plt.title(comp)
# plt.savefig(f'{figdir}/Tri21_RNA_Milo_beeswarm_sort_by_n-significant.pdf', bbox_inches="tight",dpi=300)
plt.show()
```

Example 37 (python):
```python
/tmp/ipykernel_293470/2917940356.py:25: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.
  list(data.groupby('nhood_annotation').sum().sort_values('significant',ascending=False).index)
/tmp/ipykernel_293470/2917940356.py:32: FutureWarning: 

The `scale` parameter has been renamed and will be removed in v0.15.0. Pass `density_norm='width'` for the same effect.
  sns.violinplot(data=data, x='logFC', y='nhood_annotation', inner='box', scale='width', color='lightgray')
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Fig5_Trisomy21/no03-1_combine-and-subsample_euploid-trisomy21_data.html

**Contents:**
- Import modules¶
- Read in anndata¶
- Euploid¶
- Trisomy21¶
- Select samples¶
- Euploid¶
- Trisomy 21¶
- Concatenate¶
- Subsample cells¶

**Examples:**

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

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 numpy as np
import numpy.random as random
import pandas as pd
import scanpy as sc
```

Example 2 (python):
```python
sc.settings.verbosity = 3             # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.settings.set_figure_params(dpi=80, facecolor='white', color_map = 'RdPu',)
```

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

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

Example 5 (python):
```python
import os
os.getcwd()
```

Example 6 (unknown):
```unknown
'/nfs/team205/kk18/notebooks/Foetal/Trisomy21/RNA'
```

Example 7 (python):
```python
# euploid
eup = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/Feb28ObjectRaw_finegrain_updated.h5ad')
eup
```

Example 8 (unknown):
```unknown
AnnData object with n_obs × n_vars = 297473 × 36601
    obs: 'latent_RT_efficiency', 'latent_cell_probability', 'latent_scale', 'sangerID', 'combinedID', 'donor', 'region', 'age', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'donor_by_library-prep', 'multiplexed', 'SOC | status', 'SOC | log_prob_singleton', 'SOC | log_prob_doublet', 'batch_key', '_scvi_batch', 'FACSgate', 'fine_grain', 'mid_grain', 'coarse_grain', 'sex', 'week', 'trimester', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', '_scvi_labels', 'stress_score', 'hb1_score'
    var: 'gene_name_scRNA-0', 'gene_id'
    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', 'mid_grain_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 9 (python):
```python
# subset multiome
eup = eup[eup.obs['kit_10x']=='Multiome-v1']
eup.shape
```

Example 10 (unknown):
```unknown
(211145, 36601)
```

Example 11 (python):
```python
pd.crosstab(eup.obs['age'],eup.obs['sex'])
```

Example 12 (python):
```python
pd.crosstab(eup.obs['age'],eup.obs['region'])
```

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

Example 14 (python):
```python
pd.crosstab(eup.obs['donor'],eup.obs['age'])
```

Example 15 (python):
```python
t21 = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/trisomy21/T21Hearts_cleaned_July2024.h5ad')
t21
```

Example 16 (unknown):
```unknown
AnnData object with n_obs × n_vars = 110693 × 36601
    obs: 'latent_RT_efficiency', 'latent_cell_probability', 'latent_scale', 'sangerID', 'combinedID', 'donor', 'donor_type', 'region', 'region_finest', 'age', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'flushed', 'scrublet_score', 'scrublet_leiden', 'cluster_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', 'batch_key', 'leiden_scVI_0.1', 'leiden_scVI_0.2', 'leiden_scVI_0.3', 'leiden_scVI_0.5', 'leiden_scVI_0.7', 'leiden_scVI_1', 'celltypist_coarse', 'conf_score_coarse', 'cells_to_remove', 'week', 'trimester', 'sex', 'celltypist_fine', 'conf_score_fine'
    var: 'feature_type', 'ambient_expression-0', 'ambient_expression-1', 'ambient_expression-10', 'ambient_expression-11', 'ambient_expression-12', 'ambient_expression-13', 'ambient_expression-14', 'ambient_expression-2', 'ambient_expression-3', 'ambient_expression-4', 'ambient_expression-5', 'ambient_expression-6', 'ambient_expression-7', 'ambient_expression-8', 'ambient_expression-9', 'mt', 'ribo', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'gene_id'
    uns: 'age_colors', 'cells_to_remove_colors', 'celltypist_coarse_colors', 'celltypist_fine_colors', 'donor_colors', 'leiden', 'leiden_scVI_0.1_colors', 'leiden_scVI_0.2_colors', 'leiden_scVI_0.3_colors', 'leiden_scVI_0.5_colors', 'leiden_scVI_0.7_colors', 'leiden_scVI_1_colors', 'neighbors', 'region_colors', 'sangerID_colors', 'umap'
    obsm: 'X_scVI', 'X_umap', 'latent_gene_encoding'
    obsp: 'connectivities', 'distances'
```

Example 17 (python):
```python
pd.crosstab(t21.obs['age'],t21.obs['sex'])
```

Example 18 (python):
```python
pd.crosstab(t21.obs['age'],t21.obs['region'])
```

Example 19 (python):
```python
display(pd.crosstab(t21.obs['donor'],t21.obs['region']))
display(pd.crosstab(t21.obs['donor'],t21.obs['sex']))
display(pd.crosstab(t21.obs['donor'],t21.obs['age']))
```

Example 20 (unknown):
```unknown
age_sel = ['12W0D','13W0D','14W0D']
region_sel = ['apex','atria','atrial septum','heart without node',
                 'left ventricle',
                 'node',
                 'outflow tract','right ventricle',
                 'right ventricle and papillaries',
                 'whole sample']

mask1 = eup.obs['age'].isin(age_sel)
mask2 = eup.obs['region'].isin(region_sel)
eup = eup[mask1&mask2]
eup
```

Example 21 (unknown):
```unknown
View of AnnData object with n_obs × n_vars = 55210 × 36601
    obs: 'latent_RT_efficiency', 'latent_cell_probability', 'latent_scale', 'sangerID', 'combinedID', 'donor', 'region', 'age', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'donor_by_library-prep', 'multiplexed', 'SOC | status', 'SOC | log_prob_singleton', 'SOC | log_prob_doublet', 'batch_key', '_scvi_batch', 'FACSgate', 'fine_grain', 'mid_grain', 'coarse_grain', 'sex', 'week', 'trimester', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', '_scvi_labels', 'stress_score', 'hb1_score'
    var: 'gene_name_scRNA-0', 'gene_id'
    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 22 (python):
```python
display(pd.crosstab(eup.obs['donor'],eup.obs['region']))
display(pd.crosstab(eup.obs['donor'],eup.obs['sex']))
display(pd.crosstab(eup.obs['donor'],eup.obs['age']))
```

Example 23 (unknown):
```unknown
mask = t21.obs['region'].isin(['Pericardium'])==False
t21 = t21[mask]
t21
```

Example 24 (unknown):
```unknown
View of AnnData object with n_obs × n_vars = 92284 × 36601
    obs: 'latent_RT_efficiency', 'latent_cell_probability', 'latent_scale', 'sangerID', 'combinedID', 'donor', 'donor_type', 'region', 'region_finest', 'age', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'flushed', 'scrublet_score', 'scrublet_leiden', 'cluster_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', 'batch_key', 'leiden_scVI_0.1', 'leiden_scVI_0.2', 'leiden_scVI_0.3', 'leiden_scVI_0.5', 'leiden_scVI_0.7', 'leiden_scVI_1', 'celltypist_coarse', 'conf_score_coarse', 'cells_to_remove', 'week', 'trimester', 'sex', 'celltypist_fine', 'conf_score_fine'
    var: 'feature_type', 'ambient_expression-0', 'ambient_expression-1', 'ambient_expression-10', 'ambient_expression-11', 'ambient_expression-12', 'ambient_expression-13', 'ambient_expression-14', 'ambient_expression-2', 'ambient_expression-3', 'ambient_expression-4', 'ambient_expression-5', 'ambient_expression-6', 'ambient_expression-7', 'ambient_expression-8', 'ambient_expression-9', 'mt', 'ribo', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'gene_id'
    uns: 'age_colors', 'cells_to_remove_colors', 'celltypist_coarse_colors', 'celltypist_fine_colors', 'donor_colors', 'leiden', 'leiden_scVI_0.1_colors', 'leiden_scVI_0.2_colors', 'leiden_scVI_0.3_colors', 'leiden_scVI_0.5_colors', 'leiden_scVI_0.7_colors', 'leiden_scVI_1_colors', 'neighbors', 'region_colors', 'sangerID_colors', 'umap'
    obsm: 'X_scVI', 'X_umap', 'latent_gene_encoding'
    obsp: 'connectivities', 'distances'
```

Example 25 (python):
```python
display(pd.crosstab(t21.obs['donor'],t21.obs['region']))
display(pd.crosstab(t21.obs['donor'],t21.obs['sex']))
display(pd.crosstab(t21.obs['donor'],t21.obs['age']))
```

Example 26 (python):
```python
eup.obs['euploid_or_t21'] = 'euploid'
t21.obs['euploid_or_t21'] = 'tri21'
t21.obs.replace({'cell_or_nuclei':{'Nuclei':'nuclei'}},inplace=True)

adata = eup.concatenate(t21, 
            join='inner', 
            batch_key=None, 
            index_unique=None)
adata
```

Example 27 (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 28 (unknown):
```unknown
AnnData object with n_obs × n_vars = 147494 × 36601
    obs: 'latent_RT_efficiency', 'latent_cell_probability', 'latent_scale', 'sangerID', 'combinedID', 'donor', 'region', 'age', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'donor_by_library-prep', 'multiplexed', 'SOC | status', 'SOC | log_prob_singleton', 'SOC | log_prob_doublet', 'batch_key', '_scvi_batch', 'FACSgate', 'fine_grain', 'mid_grain', 'coarse_grain', 'sex', 'week', 'trimester', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', '_scvi_labels', 'stress_score', 'hb1_score', 'euploid_or_t21', 'donor_type', 'region_finest', 'flushed', 'scrublet_leiden', 'cluster_scrublet_score', 'leiden_scVI_0.1', 'leiden_scVI_0.2', 'leiden_scVI_0.3', 'leiden_scVI_0.5', 'leiden_scVI_0.7', 'leiden_scVI_1', 'celltypist_coarse', 'conf_score_coarse', 'cells_to_remove', 'celltypist_fine', 'conf_score_fine'
    var: 'gene_id', 'gene_name_scRNA-0-0', 'feature_type-1', 'ambient_expression-0-1', 'ambient_expression-1-1', 'ambient_expression-10-1', 'ambient_expression-11-1', 'ambient_expression-12-1', 'ambient_expression-13-1', 'ambient_expression-14-1', 'ambient_expression-2-1', 'ambient_expression-3-1', 'ambient_expression-4-1', 'ambient_expression-5-1', 'ambient_expression-6-1', 'ambient_expression-7-1', 'ambient_expression-8-1', 'ambient_expression-9-1', 'mt-1', 'ribo-1', 'n_cells_by_counts-1', 'mean_counts-1', 'log1p_mean_counts-1', 'pct_dropout_by_counts-1', 'total_counts-1', 'log1p_total_counts-1'
    obsm: 'X_scVI', 'X_umap', 'latent_gene_encoding'
```

Example 29 (python):
```python
# merge coarse grain labels
adata.obs['coarse_grain_merged'] = adata.obs['coarse_grain'].astype('str').copy()
t21_obsnames = adata.obs_names[adata.obs['euploid_or_t21']=='tri21']
adata.obs.loc[t21_obsnames,'coarse_grain_merged'] = adata.obs.loc[t21_obsnames,'celltypist_coarse']
set(adata.obs['coarse_grain_merged'])
```

Example 30 (unknown):
```unknown
{'Cardiomyocytes',
 'Endothelium',
 'Epicardium',
 'Leukocytes',
 'Mesenchymal',
 'Neural'}
```

Example 31 (python):
```python
# merge fine grain labels
adata.obs['fine_grain_merged'] = adata.obs['fine_grain'].astype('str').copy()
t21_obsnames = adata.obs_names[adata.obs['euploid_or_t21']=='tri21']
adata.obs.loc[t21_obsnames,'fine_grain_merged'] = adata.obs.loc[t21_obsnames,'celltypist_fine']
print(f"number of NaN in the merged label: {sum(adata.obs['fine_grain_merged'].isna())}")
with pd.option_context('display.max_rows', None):
    display(pd.crosstab(adata.obs['fine_grain_merged'],adata.obs['euploid_or_t21']))
```

Example 32 (unknown):
```unknown
number of NaN in the merged label: 0
```

Example 33 (python):
```python
# save concatenated object
adata.write('/nfs/team205/heart/anndata_objects/Foetal/trisomy21/Euploid_T21Hearts_Aug2024_sel.h5ad')
```

Example 34 (python):
```python
# read in cell IDs which passed ATAC QC: trisomy21
tri21_atac_cells = utils.txt2list("/nfs/team205/heart/anndata_objects/Foetal/trisomy21/ArchR/project_output/project_cellIDs_removed-low-cellnumer-celltypes.txt")
print(len(tri21_atac_cells))
print(tri21_atac_cells[:5])

# read in cell IDs which passed ATAC QC: euploid
eup_atac_cells = utils.txt2list('/nfs/team205/heart/anndata_objects/Foetal/multiome_ATAC/ArchR/project_output_peak-revised-Feb2024/project_cellIDs.txt')
print(len(eup_atac_cells))
print(eup_atac_cells[:5])

# conbine
atac_cells = eup_atac_cells + tri21_atac_cells
print(len(atac_cells))
```

Example 35 (unknown):
```unknown
89416
['BHF_F_Hea13188320_BHF_F_Hea13187622_GGGTCACTCATAGCCG-1', 'BHF_F_Hea13188320_BHF_F_Hea13187622_TAATCACCATTGGGAG-1', 'BHF_F_Hea13188320_BHF_F_Hea13187622_CGAGCTGGTCTCACTG-1', 'BHF_F_Hea13188320_BHF_F_Hea13187622_AAGGTATAGGCTACTG-1', 'BHF_F_Hea13188320_BHF_F_Hea13187622_GCGCGATTCATAATCG-1']
167022
['BHF_F_Hea11064670_BHF_F_Hea11031823_CAGCCTAAGTCTTGAA-1', 'BHF_F_Hea11064670_BHF_F_Hea11031823_GGTGTTGTCAGGCCTA-1', 'BHF_F_Hea11064670_BHF_F_Hea11031823_CGTGTGTCAAGCCACT-1', 'BHF_F_Hea11064670_BHF_F_Hea11031823_AGCTTTAAGCTAAAGG-1', 'BHF_F_Hea11064670_BHF_F_Hea11031823_GTTTCCTCAGGTCCTG-1']
256438
```

Example 36 (python):
```python
shared_obsnames = list(set(adata.obs_names).intersection(atac_cells))
print(len(shared_obsnames))
```

Example 37 (python):
```python
print(adata.shape)
adata_sub = adata[shared_obsnames]
print(adata_sub.shape)
```

Example 38 (unknown):
```unknown
(147494, 36601)
(122174, 36601)
```

Example 39 (python):
```python
with pd.option_context('display.max_rows', None):
    display(pd.crosstab(adata_sub.obs['fine_grain_merged'],adata_sub.obs['euploid_or_t21']))
```

Example 40 (python):
```python
with pd.option_context('display.max_rows', None):
    display(pd.crosstab(adata_sub.obs['fine_grain_merged'],adata_sub.obs['donor']))
```

Example 41 (python):
```python
# select cell states which has more than 30 cells in both euploid and trisomy21
ctab = pd.crosstab(adata_sub.obs['fine_grain_merged'],adata_sub.obs['euploid_or_t21'])
mask = (ctab>=30).sum(axis=1)==2
celltype_sel = ctab.index[mask]

# for the celltypes which passed above, subset nuclei per donor to be the 3rd largest number in the list of nuclei number per donor
selected_obsnames = []
for celltype in celltype_sel:
    ad = adata_sub[adata_sub.obs['fine_grain_merged']==celltype]
    n_subsample = ad.obs['donor'].value_counts()[2] # get 3rd largest number
    ad = utils.sctk_subsample(ad, fraction=1, groupby='donor', max_n=n_subsample, random_state=0)
    selected_obsnames = selected_obsnames+list(ad.obs_names)
    del ad
adata_sub = adata_sub[selected_obsnames]
adata_sub
```

Example 42 (unknown):
```unknown
View of AnnData object with n_obs × n_vars = 83477 × 36601
    obs: 'latent_RT_efficiency', 'latent_cell_probability', 'latent_scale', 'sangerID', 'combinedID', 'donor', 'region', 'age', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'donor_by_library-prep', 'multiplexed', 'SOC | status', 'SOC | log_prob_singleton', 'SOC | log_prob_doublet', 'batch_key', '_scvi_batch', 'FACSgate', 'fine_grain', 'mid_grain', 'coarse_grain', 'sex', 'week', 'trimester', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', '_scvi_labels', 'stress_score', 'hb1_score', 'euploid_or_t21', 'donor_type', 'region_finest', 'flushed', 'scrublet_leiden', 'cluster_scrublet_score', 'leiden_scVI_0.1', 'leiden_scVI_0.2', 'leiden_scVI_0.3', 'leiden_scVI_0.5', 'leiden_scVI_0.7', 'leiden_scVI_1', 'celltypist_coarse', 'conf_score_coarse', 'cells_to_remove', 'celltypist_fine', 'conf_score_fine', 'coarse_grain_merged', 'fine_grain_merged'
    var: 'gene_id', 'gene_name_scRNA-0-0', 'feature_type-1', 'ambient_expression-0-1', 'ambient_expression-1-1', 'ambient_expression-10-1', 'ambient_expression-11-1', 'ambient_expression-12-1', 'ambient_expression-13-1', 'ambient_expression-14-1', 'ambient_expression-2-1', 'ambient_expression-3-1', 'ambient_expression-4-1', 'ambient_expression-5-1', 'ambient_expression-6-1', 'ambient_expression-7-1', 'ambient_expression-8-1', 'ambient_expression-9-1', 'mt-1', 'ribo-1', 'n_cells_by_counts-1', 'mean_counts-1', 'log1p_mean_counts-1', 'pct_dropout_by_counts-1', 'total_counts-1', 'log1p_total_counts-1'
    obsm: 'X_scVI', 'X_umap', 'latent_gene_encoding'
```

Example 43 (python):
```python
with pd.option_context('display.max_rows', None):
    display(pd.crosstab(adata_sub.obs['fine_grain_merged'],adata_sub.obs['euploid_or_t21']))

with pd.option_context('display.max_rows', None):
    display(pd.crosstab(adata_sub.obs['fine_grain_merged'],adata_sub.obs['donor']))
```

Example 44 (python):
```python
pd.crosstab(adata_sub.obs['donor'],adata_sub.obs['euploid_or_t21'])
```

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

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

Example 47 (python):
```python
# save
adata_sub.write('/nfs/team205/heart/anndata_objects/Foetal/trisomy21/Euploid_T21Hearts_Aug2024_sel_subsampled-per-cellstate-donor.h5ad')
```

Example 48 (python):
```python
# save obs for later use
adata_sub.obs.to_csv('/nfs/team205/heart/anndata_objects/Foetal/trisomy21/Euploid_T21Hearts_Aug2024_sel_subsampled-per-cellstate-donor.obs.csv')
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Fig5_Trisomy21/no13_GSEA_summary_finegrain.html

**Contents:**
- Read in GSEA results¶
- Plot: Euploid nhoods vs reference nhoods¶

**Examples:**

Example 1 (python):
```python
import scanpy as sc
import pandas as pd
import numpy as np
import pandas as pd
import anndata
import scipy
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
import warnings
warnings.simplefilter("ignore", UserWarning)
warnings.simplefilter(action='ignore', category=pd.errors.PerformanceWarning)
```

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

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

Example 4 (python):
```python
gsea_res_df = pd.read_csv('/nfs/team205/kk18/notebooks/Foetal/Trisomy21/RNA/milo_Mar2025_refine-nhoods/GSEA_trisomy21-vs-euploid-nhoods_finegrain.csv',
                                index_col=0)
print(gsea_res_df.shape)
gsea_res_df.head()
```

Example 5 (unknown):
```unknown
(20860, 10)
```

Example 6 (python):
```python
from matplotlib.colors import ListedColormap
```

Example 7 (python):
```python
fdr_thresh = 0.1
geneset = 'MSigDB_Hallmark_2020'
plt_title = f'Trisomy21-nhood vs Euploid-nhood, {geneset}'
significant_only = True
figsize = (8,5.5)

#################################
# prepare dataframes
df = gsea_res_df[gsea_res_df['geneset'] == geneset]
df_nes = df.pivot(columns='Term',index='celltype',values='nes')
df_nes = df_nes.fillna(0)
df_sig = df.pivot(columns='Term',index='celltype',values='fdr')
df_sig = df_sig<fdr_thresh

if significant_only:
    # select pathways which have significant NES at least in one cell type
    pathways_sel = list(df_sig.columns[df_sig.sum()>0])
    df_nes = df_nes[pathways_sel]
    df_sig = df_sig[pathways_sel]

# plotting
plt.rcParams['axes.grid'] = False
# Custom colormap: add grey for masked areas
g = sns.clustermap(df_nes,figsize=figsize,cmap='coolwarm',center=0, # vmax=5,
                  dendrogram_ratio=0.05, cbar_pos=(-0.07, 0.9, 0.02, 0.10)
                  ) # 

# Get the reordered indices for rows and columns
row_order = g.dendrogram_row.reordered_ind
col_order = g.dendrogram_col.reordered_ind

# Reorder the mask based on the clustering result
df_sig_reorder = df_sig.iloc[row_order, :].iloc[:, col_order]

# Iterate over the mask and add asterisks
for i in range(df_nes.shape[0]):
    for j in range(df_nes.shape[1]):
        if df_sig_reorder.iloc[i, j]:
            g.ax_heatmap.text(j + 0.5, i + 0.7, '*', ha='center', va='center', color='black', fontsize=14)
            
# Add a title to the colorbar
colorbar = g.ax_heatmap.collections[0].colorbar
colorbar.set_label('NES', rotation=270, labelpad=15)
plt.title(plt_title)
plt.show()
```

Example 8 (python):
```python
fdr_thresh = 0.1
geneset = 'KEGG_2021_Human'
plt_title = f'Trisomy21-nhood vs Euploid-nhood, {geneset}'
significant_only = True
figsize = (12,6)

#################################
# prepare dataframes
df = gsea_res_df[gsea_res_df['geneset'] == geneset]
df_nes = df.pivot(columns='Term',index='celltype',values='nes')
df_nes = df_nes.fillna(0)
df_sig = df.pivot(columns='Term',index='celltype',values='fdr')
df_sig = df_sig<fdr_thresh

if significant_only:
    # select pathways which have significant NES at least in one cell type
    pathways_sel = list(df_sig.columns[df_sig.sum()>0])
    df_nes = df_nes[pathways_sel]
    df_sig = df_sig[pathways_sel]

# plotting
plt.rcParams['axes.grid'] = False
# Custom colormap: add grey for masked areas
g = sns.clustermap(df_nes,figsize=figsize,cmap='coolwarm',center=0, # vmax=5,
                  dendrogram_ratio=0.05, cbar_pos=(-0.07, 0.9, 0.02, 0.10)
                  ) # 

# Get the reordered indices for rows and columns
row_order = g.dendrogram_row.reordered_ind
col_order = g.dendrogram_col.reordered_ind

# Reorder the mask based on the clustering result
df_sig_reorder = df_sig.iloc[row_order, :].iloc[:, col_order]

# Iterate over the mask and add asterisks
for i in range(df_nes.shape[0]):
    for j in range(df_nes.shape[1]):
        if df_sig_reorder.iloc[i, j]:
            g.ax_heatmap.text(j + 0.5, i + 0.7, '*', ha='center', va='center', color='black', fontsize=14)
            
# Add a title to the colorbar
colorbar = g.ax_heatmap.collections[0].colorbar
colorbar.set_label('NES', rotation=270, labelpad=15)
plt.title(plt_title)
plt.show()
```

Example 9 (python):
```python
fdr_thresh = 0.1
genesets = ['MSigDB_Hallmark_2020','KEGG_2021_Human']
remove_keywords = ['infection','disease','cancer','cardiomyopathy','shigellosis',
                   'malaria','hepatitis','bacterial','diabetes','atherosclerosis','sclerosis',
                   'carcinogenesis'
                  ]
plt_title = f'Trisomy21-nhood vs Euploid-nhood, {geneset}'
significant_only = True
figsize = (12,6.5)
savepath = f'{figdir}/Tri21_Milo-reassigned-DEGs_T21-vs-EUPnhoods_GSEA_MSigDB-KEGG.pdf'

#################################
# prepare dataframes
df = gsea_res_df.copy()
df = df[df['geneset'].isin(genesets)]
df.replace({'geneset':{
    'MSigDB_Hallmark_2020':'MH',
    'KEGG_2021_Human':'KEGG',
    'GO_Biological_Process_2023':'GOBP',
    'Reactome_2022':'RT'
}},inplace=True)
df['Term'] = df['geneset'].astype('str')+': '+df['Term'].astype('str')
df_nes = df.pivot(columns='Term',index='celltype',values='nes')
df_nes = df_nes.fillna(0)
df_sig = df.pivot(columns='Term',index='celltype',values='fdr')
df_sig = df_sig<fdr_thresh

if significant_only:
    # select pathways which have significant NES at least in one cell type
    pathways_sel = list(df_sig.columns[df_sig.sum()>0])
    df_nes = df_nes[pathways_sel]
    df_sig = df_sig[pathways_sel]

if remove_keywords!=None:
    pathways_sel = [x for x in df_nes.columns if all([xx.lower() not in x.lower() for xx in remove_keywords])]
    df_nes = df_nes[pathways_sel]
    df_sig = df_sig[pathways_sel]

# plotting
plt.rcParams['axes.grid'] = False
# Custom colormap: add grey for masked areas
g = sns.clustermap(df_nes,figsize=figsize,cmap='coolwarm',center=0, # vmax=5,
                  dendrogram_ratio=0.05, cbar_pos=(-0.07, 0.9, 0.02, 0.10)
                  ) # 

# Get the reordered indices for rows and columns
row_order = g.dendrogram_row.reordered_ind
col_order = g.dendrogram_col.reordered_ind

# Reorder the mask based on the clustering result
df_sig_reorder = df_sig.iloc[row_order, :].iloc[:, col_order]

# Iterate over the mask and add asterisks
for i in range(df_nes.shape[0]):
    for j in range(df_nes.shape[1]):
        if df_sig_reorder.iloc[i, j]:
            g.ax_heatmap.text(j + 0.5, i + 0.7, '*', ha='center', va='center', color='black', fontsize=14)
            
# Add a title to the colorbar
colorbar = g.ax_heatmap.collections[0].colorbar
colorbar.set_label('NES', rotation=270, labelpad=15)
plt.title(plt_title)
if savepath!=None:
    plt.savefig(savepath, bbox_inches="tight",dpi=300)
plt.show()
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Fig5_Trisomy21/no11_deg_finegrain.html

**Contents:**
- Read in adata¶
- DEG test for each celltype, logFC_direction: euploid vs reference¶

42412 rows × 7 columns

**Examples:**

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

import warnings
warnings.simplefilter("ignore", UserWarning)
warnings.simplefilter(action='ignore', category=pd.errors.PerformanceWarning)
```

Example 2 (python):
```python
# Configure rpy2 to Use the Correct R Installation:
import os
os.environ['R_HOME'] = '/software/cellgen/team205/kk18/envs/generic_env/lib/R'
# check
from rpy2.rinterface_lib import openrlib
print(openrlib.R_HOME)
```

Example 3 (unknown):
```unknown
/software/cellgen/team205/kk18/envs/generic_env/lib/R
```

Example 4 (python):
```python
import rpy2.rinterface_lib.callbacks
import anndata2ri
import logging

from rpy2.robjects import pandas2ri
from rpy2.robjects import r

sc.settings.verbosity = 0
rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)

pandas2ri.activate()
anndata2ri.activate()

%load_ext rpy2.ipython
```

Example 5 (unknown):
```unknown
/tmp/ipykernel_41476/3191876078.py:12: DeprecationWarning: The global conversion available with activate() is deprecated and will be removed in the next major release. Use a local converter.
  anndata2ri.activate()
```

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

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

Example 8 (unknown):
```unknown
%%R
library(edgeR)
library(tidyverse)
```

Example 9 (unknown):
```unknown
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
```

Example 10 (unknown):
```unknown
Loading required package: limma
```

Example 11 (python):
```python
# pseudobulked per nhood
# nhood-by-gene
adata = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/trisomy21/Euploid_T21Hearts_Aug2024_sel_subsampled-per-cellstate-donor_cells-in-sig-nhoods-sfdr0p1_pseudobulked-nhood-reassigned-sum_finegrain.h5ad')
adata
```

Example 12 (unknown):
```unknown
AnnData object with n_obs × n_vars = 577 × 36601
    obs: 'index_cell', 'logFC', 'PValue', 'FDR', 'SpatialFDR', 'nhood_annotation', 'deg_group'
    var: 'gene_id', 'gene_name_scRNA-0-0', 'feature_type-1', 'mt-1', 'ribo-1', 'n_cells_by_counts-1', 'mean_counts-1', 'log1p_mean_counts-1', 'pct_dropout_by_counts-1', 'total_counts-1', 'log1p_total_counts-1', 'gene_name_GenomicFeatures', 'gene_biotype_GenomicFeatures', 'gene_seq_start_GenomicFeatures', 'gene_seq_end_GenomicFeatures', 'seq_name_GenomicFeatures', 'seq_strand_GenomicFeatures', 'seq_coord_system_GenomicFeatures', 'description_GenomicFeatures', 'gene_id_version_GenomicFeatures', 'canonical_transcript_GenomicFeatures', 'gene_name'
```

Example 13 (python):
```python
celltype_list = [x for x in adata.obs['nhood_annotation'].unique() if x!='Mixed']
celltype_list
```

Example 14 (unknown):
```unknown
['AtrialCardiomyocytesRight',
 'AtrialCardiomyocytesCycling',
 'VentricularConductionSystemProximal',
 'AtrialCardiomyocytesLeft',
 'AtrioventricularNodePacemakerCells',
 'CoronaryCapillaryEndothelialCells',
 'CoronaryPericytes',
 'CoronaryVesselAdventitialFibroblasts',
 'DuctusArteriosusSmoothMuscleCells',
 'EndocardialCells',
 'EndocardialCushionCells',
 'SubEpicardialFibroblasts',
 'GreatVesselAdventitialFibroblasts',
 'GreatVesselArterialEndothelialCells',
 'GreatVesselSmoothMuscleCells',
 'MacrophagesATF3pos',
 'MacrophagesLYVE1pos',
 'Megakaryocytes',
 'MesothelialEpicardialCells',
 'MyocardialInterstitialFibroblasts',
 'ValveInterstitialCells',
 'Myofibroblasts',
 'ParasympatheticNeurons',
 'VentricularCardiomyocytesCycling',
 'VentricularCardiomyocytesRightTrabeculated',
 'VentricularCardiomyocytesLeftCompact',
 'VentricularCardiomyocytesLeftTrabeculated',
 'VentricularCardiomyocytesRightCompact']
```

Example 15 (python):
```python
# select celltypes which has n_sample>=3
ctab = pd.crosstab(adata.obs['nhood_annotation'],adata.obs['deg_group'])[['euploid','trisomy21']]
celltypes_sel = list(ctab.index[(ctab>=3).sum(axis=1)==2])
celltypes_sel = [x for x in celltypes_sel if x!='Mixed']
adata_sub = adata[adata.obs['nhood_annotation'].isin(celltypes_sel)]
# subset groups
adata_sub = adata_sub[adata_sub.obs['deg_group'].isin(['euploid','trisomy21'])]
pd.crosstab(adata_sub.obs['nhood_annotation'],adata_sub.obs['deg_group'])
```

Example 16 (python):
```python
%%R -i adata_sub
class(adata_sub)
```

Example 17 (unknown):
```unknown
[1] "SingleCellExperiment"
attr(,"package")
[1] "SingleCellExperiment"
```

Example 18 (python):
```python
%%R
# https://github.com/compbiomed/singleCellTK/blob/40527afd2e58e14451aae655c0d8fae4a7a0f54b/R/subsetSCE.R#L140
subsetSCECols <- function(inSCE, index = NULL, bool = NULL, colData = NULL) {

  if(is.null(index) & is.null(bool) & is.null(colData)) {
    stop("At least one of 'index', 'bool', or 'colData' must be supplied.")
  }
  final.ix <- rep(TRUE, ncol(inSCE))

  # Parse index containing integers
  if(!is.null(index)) {
    if(min(index) < 1 | max(index) > ncol(inSCE)) {
      stop("'index' must contain integers between 1 and the number of columns ",
           "in 'inSCE': ", ncol(inSCE))
    }
    final.ix[!seq(ncol(inSCE)) %in% index] <- FALSE
  }

  # Parse Boolean vector
  if(!is.null(bool)) {
    if(length(bool) != ncol(inSCE) | !is.logical(bool)) {
      stop("'bool' must be a logical vector the same length as the number of ",
           "colmns in 'inSCE': ", ncol(inSCE))
    }
    final.ix[!bool] <- FALSE
  }

  # Parse expressions for colData variables
  if(!is.null(colData)) {
    for(i in seq_along(colData)) {
      temp <- eval(parse(text = as.character(colData[i])),
                   envir = as.data.frame(colData(inSCE)))
      if(length(temp) != ncol(inSCE) | !is.logical(temp)) {
        stop("The expression ", colData[i], " did not produce a boolean ",
             "vector the same length as the number of columns in 'inSCE'. ",
             "Please ensure that the spelling of the variable you are ",
             "trying to use matches one of the column names in ",
             "'colData(inSCE)' and that your expression is valid.")
      }
      final.ix[!temp] <- FALSE
    }
  }

  inSCE <- inSCE[,final.ix]
  return(inSCE)
  }

# https://www.sc-best-practices.org/conditions/differential_gene_expression.html
fit_model <- function(adata_){
    # create an edgeR object with counts and grouping factor
    y <- DGEList(assay(adata_, "X"), group = colData(adata_)$deg_group)
    # filter out genes with low counts
    print("Dimensions before subsetting:")
    print(dim(y))
    print("")
    keep <- filterByExpr(y)
    y <- y[keep, , keep.lib.sizes=FALSE]
    print("Dimensions after subsetting:")
    print(dim(y))
    print("")
    # normalize
    y <- calcNormFactors(y)
    # create a vector that is a concatenation of condition and cell type that we will later use with contrasts
    group <- colData(adata_)$deg_group
    # sex <- colData(adata_)$sex
    # replicate <- colData(adata_)$donor
    # create a design matrix: here we have multiple donors so also consider that in the design matrix
    design <- model.matrix(~ 0 + group)
    # estimate dispersion
    y <- estimateDisp(y, design = design)
    # fit the model
    fit <- glmQLFit(y, design)
    return(list("fit"=fit, "design"=design, "y"=y))
}
```

Example 19 (python):
```python
%%time
%%R -o tt_all 
celltype_list <- unique(colData(adata_sub)$nhood_annotation)
for(celltype in celltype_list){
    print(celltype)
    # subset
    sce <- subsetSCECols(adata_sub, colData = "nhood_annotation == celltype")
    # fit model
    outs <-fit_model(sce)
    # get output
    fit <- outs$fit
    y <- outs$y
    # deg
    myContrast <- makeContrasts('grouptrisomy21-groupeuploid', levels = y$design)
    qlf <- glmQLFTest(fit, contrast=myContrast)
    # get all of the DE genes and calculate Benjamini-Hochberg adjusted FDR
    tt <- topTags(qlf, n = Inf)
    tt <- tt$table
    tt <- rownames_to_column(tt, var = "gene_name")
    tt$celltype <- celltype
    # rbind
    if(celltype==celltype_list[1]){
        tt_all = tt
    }else{
        tt_all = rbind(tt_all,tt)
    }
    print('')
}
```

Example 20 (unknown):
```unknown
[1] "AtrialCardiomyocytesRight"
[1] "Dimensions before subsetting:"
[1] 36601    70
[1] ""
[1] "Dimensions after subsetting:"
[1] 3761   70
[1] ""
[1] ""
[1] "AtrialCardiomyocytesLeft"
[1] "Dimensions before subsetting:"
[1] 36601    36
[1] ""
[1] "Dimensions after subsetting:"
[1] 3740   36
[1] ""
[1] ""
[1] "CoronaryCapillaryEndothelialCells"
[1] "Dimensions before subsetting:"
[1] 36601    11
[1] ""
[1] "Dimensions after subsetting:"
[1] 2684   11
[1] ""
[1] ""
[1] "EndocardialCells"
[1] "Dimensions before subsetting:"
[1] 36601    12
[1] ""
[1] "Dimensions after subsetting:"
[1] 5171   12
[1] ""
[1] ""
[1] "SubEpicardialFibroblasts"
[1] "Dimensions before subsetting:"
[1] 36601    15
[1] ""
[1] "Dimensions after subsetting:"
[1] 2628   15
[1] ""
[1] ""
[1] "GreatVesselAdventitialFibroblasts"
[1] "Dimensions before subsetting:"
[1] 36601    16
[1] ""
[1] "Dimensions after subsetting:"
[1] 3758   16
[1] ""
[1] ""
[1] "GreatVesselSmoothMuscleCells"
[1] "Dimensions before subsetting:"
[1] 36601    22
[1] ""
[1] "Dimensions after subsetting:"
[1] 2607   22
[1] ""
[1] ""
[1] "MacrophagesATF3pos"
[1] "Dimensions before subsetting:"
[1] 36601    12
[1] ""
[1] "Dimensions after subsetting:"
[1] 1570   12
[1] ""
[1] ""
[1] "MacrophagesLYVE1pos"
[1] "Dimensions before subsetting:"
[1] 36601    44
[1] ""
[1] "Dimensions after subsetting:"
[1] 2509   44
[1] ""
[1] ""
[1] "MyocardialInterstitialFibroblasts"
[1] "Dimensions before subsetting:"
[1] 36601    43
[1] ""
[1] "Dimensions after subsetting:"
[1] 2267   43
[1] ""
[1] ""
[1] "VentricularCardiomyocytesCycling"
[1] "Dimensions before subsetting:"
[1] 36601    44
[1] ""
[1] "Dimensions after subsetting:"
[1] 3294   44
[1] ""
[1] ""
[1] "VentricularCardiomyocytesLeftCompact"
[1] "Dimensions before subsetting:"
[1] 36601    46
[1] ""
[1] "Dimensions after subsetting:"
[1] 2269   46
[1] ""
[1] ""
[1] "VentricularCardiomyocytesLeftTrabeculated"
[1] "Dimensions before subsetting:"
[1] 36601     8
[1] ""
[1] "Dimensions after subsetting:"
[1] 2930    8
[1] ""
[1] ""
[1] "VentricularCardiomyocytesRightCompact"
[1] "Dimensions before subsetting:"
[1] 36601    20
[1] ""
[1] "Dimensions after subsetting:"
[1] 3224   20
[1] ""
[1] ""
CPU times: user 16.5 s, sys: 244 ms, total: 16.7 s
Wall time: 17.9 s
```

Example 21 (unknown):
```unknown
# save
tt_all.to_csv('./milo_nhood-reassigned_pseudobulk_edgeR_res_tri21-nhoods_vs_euploid-nhoods_finegrain.csv')
```

---

## Notebook

**URL:** http://127.0.0.1:9020/notebooks_html/Fig5_Trisomy21/no09_pseudobulk_finegrain.html

**Contents:**
- Read in adata¶
- Annotate nhoods¶
- Pseudobulk: nhood-by-gene adata¶

**Examples:**

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

import seaborn as sns

import warnings
warnings.simplefilter("ignore", UserWarning)
```

Example 2 (python):
```python
# Configure rpy2 to Use the Correct R Installation:
import os
os.environ['R_HOME'] = '/software/cellgen/team205/kk18/envs/generic_env/lib/R'
# check
from rpy2.rinterface_lib import openrlib
print(openrlib.R_HOME)
```

Example 3 (unknown):
```unknown
/software/cellgen/team205/kk18/envs/generic_env/lib/R
```

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

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

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

Example 7 (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
# atlas
adata = sc.read_h5ad('/nfs/team205/heart/anndata_objects/Foetal/trisomy21/Euploid_T21Hearts_Aug2024_sel_subsampled-per-cellstate-donor_cells-in-sig-nhoods-sfdr0p1.h5ad')
adata
```

Example 9 (python):
```python
AnnData object with n_obs × n_vars = 14846 × 36601
    obs: 'latent_RT_efficiency', 'latent_cell_probability', 'latent_scale', 'sangerID', 'combinedID', 'donor', 'region', 'age', 'facility', 'cell_or_nuclei', 'modality', 'kit_10x', 'scrublet_score', 'doublet_pval', 'doublet_bh_pval', 'n_genes', 'n_counts', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'HB_score', 'donor_by_library-prep', 'multiplexed', 'SOC | status', 'SOC | log_prob_singleton', 'SOC | log_prob_doublet', 'batch_key', '_scvi_batch', 'FACSgate', 'fine_grain', 'mid_grain', 'coarse_grain', 'sex', 'week', 'trimester', 'heart_or_greatvessels', 'cycling', 'S_score', 'G2M_score', 'phase', '_scvi_labels', 'stress_score', 'hb1_score', 'euploid_or_t21', 'donor_type', 'region_finest', 'flushed', 'scrublet_leiden', 'cluster_scrublet_score', 'leiden_scVI_0.1', 'leiden_scVI_0.2', 'leiden_scVI_0.3', 'leiden_scVI_0.5', 'leiden_scVI_0.7', 'leiden_scVI_1', 'celltypist_coarse', 'conf_score_coarse', 'cells_to_remove', 'celltypist_fine', 'conf_score_fine', 'coarse_grain_merged', 'fine_grain_merged', 'nhood_ixs_random', 'nhood_ixs_refined', 'nhood_kth_distance'
    var: 'gene_id', 'gene_name_scRNA-0-0', 'feature_type-1', 'mt-1', 'ribo-1', 'n_cells_by_counts-1', 'mean_counts-1', 'log1p_mean_counts-1', 'pct_dropout_by_counts-1', 'total_counts-1', 'log1p_total_counts-1', 'gene_name_GenomicFeatures', 'gene_biotype_GenomicFeatures', 'gene_seq_start_GenomicFeatures', 'gene_seq_end_GenomicFeatures', 'seq_name_GenomicFeatures', 'seq_strand_GenomicFeatures', 'seq_coord_system_GenomicFeatures', 'description_GenomicFeatures', 'gene_id_version_GenomicFeatures', 'canonical_transcript_GenomicFeatures', 'gene_name'
    uns: 'coarse_grain_merged_colors', 'correcting-donor-region_n-layers-2', 'donor_colors', 'euploid_or_t21_colors', 'fine_grain_merged_colors', 'mid_grain_colors', 'nhood_adata', 'nhood_neighbors_key', 'region_colors', 'scVI_latent_correcting-donor-region_n-layers-2-forMilo', 'umap'
    obsm: 'X_umap_correcting-donor-region_n-layers-2', 'latent_gene_encoding', 'nhoods', 'nhoods_reassigned', 'scVI_latent_correcting-donor-region_n-layers-2'
    obsp: 'correcting-donor-region_n-layers-2_connectivities', 'correcting-donor-region_n-layers-2_distances', 'scVI_latent_correcting-donor-region_n-layers-2-forMilo_connectivities', 'scVI_latent_correcting-donor-region_n-layers-2-forMilo_distances'
```

Example 10 (python):
```python
milopy.utils.annotate_nhoods(adata, anno_col='fine_grain_merged')
```

Example 11 (python):
```python
plt.hist(adata.uns['nhood_adata'].obs["nhood_annotation_frac"],bins=100);
plt.xlabel("celltype fraction")
```

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

Example 13 (python):
```python
# We can see that for the majority of neighbourhoods, almost all cells have the same neighbourhood. We can rename neighbourhoods where less than 60% of the cells have the top label as "Mixed"
adata.uns['nhood_adata'].obs["nhood_annotation"] = adata.uns['nhood_adata'].obs["nhood_annotation"].astype('str')
adata.uns['nhood_adata'].obs.loc[adata.uns['nhood_adata'].obs["nhood_annotation_frac"] < 0.6, "nhood_annotation"] = "Mixed"
```

Example 14 (python):
```python
adata.uns['nhood_adata'].obs['deg_group'].value_counts()
```

Example 15 (unknown):
```unknown
deg_group
trisomy21    339
euploid      238
Name: count, dtype: int64
```

Example 16 (python):
```python
%%time
# first take sum
# with raw gene counts
# np.dot of gene-by-cell and cell-by-nhood (boolean) will make gene-by-nhood with summed counts
X = adata.X.todense().T.copy() # gene-by-cell
Y = adata.obsm['nhoods_reassigned'].todense().copy()
Z = np.matmul(X,Y) # gene-by-nhood

# prepare anndata
adata_nhood = anndata.AnnData(
    X=np.array(Z).T,
    obs=adata.uns['nhood_adata'].obs[['index_cell','logFC','PValue','FDR','SpatialFDR','nhood_annotation','deg_group']].copy(),
    var=adata.var.copy()
)
adata_nhood
```

Example 17 (unknown):
```unknown
CPU times: user 18.1 s, sys: 3.45 s, total: 21.6 s
Wall time: 22.6 s
```

Example 18 (unknown):
```unknown
AnnData object with n_obs × n_vars = 577 × 36601
    obs: 'index_cell', 'logFC', 'PValue', 'FDR', 'SpatialFDR', 'nhood_annotation', 'deg_group'
    var: 'gene_id', 'gene_name_scRNA-0-0', 'feature_type-1', 'mt-1', 'ribo-1', 'n_cells_by_counts-1', 'mean_counts-1', 'log1p_mean_counts-1', 'pct_dropout_by_counts-1', 'total_counts-1', 'log1p_total_counts-1', 'gene_name_GenomicFeatures', 'gene_biotype_GenomicFeatures', 'gene_seq_start_GenomicFeatures', 'gene_seq_end_GenomicFeatures', 'seq_name_GenomicFeatures', 'seq_strand_GenomicFeatures', 'seq_coord_system_GenomicFeatures', 'description_GenomicFeatures', 'gene_id_version_GenomicFeatures', 'canonical_transcript_GenomicFeatures', 'gene_name'
```

Example 19 (python):
```python
# save
adata_nhood.write('/nfs/team205/heart/anndata_objects/Foetal/trisomy21/Euploid_T21Hearts_Aug2024_sel_subsampled-per-cellstate-donor_cells-in-sig-nhoods-sfdr0p1_pseudobulked-nhood-reassigned-sum_finegrain.h5ad')
```

---
