Tutorial for MERFISH WMB dataset
Need additional packages: scanpy seaborn
Load the packages
[2]:
%reload_ext autoreload
%autoreload 2
import os
import time
import scanpy as sc
import pandas as pd
import numpy as np
import anndata as ad
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.patches import Patch
from matplotlib.lines import Line2D
from sklearn.metrics import adjusted_rand_score, adjusted_mutual_info_score, f1_score
from Harmonics import *
import warnings
warnings.filterwarnings("ignore")
sc.settings.verbosity = 0
sc.settings.set_figure_params(dpi=30, dpi_save=500)
from matplotlib import rcParams
rcParams["figure.dpi"] = 30
rcParams["savefig.dpi"] = 500
[ ]:
data_dir = '../../../Data/Spatial/Transcriptomics/MERFISH_WMB_Zhang2023/animal3_sagittal/processed/'
save_dir = f'../../results/MERFISH_WMB_Zhang2023/animal3_sagittal/Harmonics/'
if not os.path.exists(save_dir):
os.makedirs(save_dir)
Define the function to change p values to corresponding star representation, used to show the results of additional tests implemented in Harmonics.
[4]:
def p2stars(p):
if p < 0.001:
return '***'
elif p < 0.01:
return '**'
elif p < 0.05:
return '*'
else:
return ''
Load dataset
[9]:
slice_name_list = np.loadtxt(data_dir + f"slice_name_list.txt", dtype=str, delimiter=" ").tolist()
adata_list = []
for slice_name in slice_name_list:
adata = ad.read_h5ad(data_dir + slice_name + '.h5ad')
adata_list.append(adata)
Run Harmonics
Instantiate Harmonics
[5]:
model = Harmonics_Model(adata_list,
slice_name_list,
cond_list=None, # default
cond_name_list=None, # default
concat_label='slice_name', # default
proportion_label=None, # default
seed=1234, # default
parallel=True, # recommand to set to True for large dataset and False for small dataset
verbose=True, # default
)
Dataset comprises 23 slices, 2081549 cells/spots in total.
Preprocess the data (Generating the connection graph and calculating neighborhood cell type destribution for cells)
[6]:
model.preprocess(ct_key='subclass_transfer',
spatial_key='spatial', # default
method='joint', # default
n_step=3, # default
n_neighbors=20, # default
cut_percentage=99, # default
)
Generating Delaunay neighbor graph...
100%|██████████| 23/23 [00:27<00:00, 1.19s/it]
All done!
Performing graph completion...
100%|██████████| 23/23 [04:08<00:00, 10.82s/it]
All done!
The cell types of interest are:
ABC NN
ACB-BST-FS D1 Gaba
AD Serpinb7 Glut
ADP-MPO Trp73 Glut
AHN Onecut3 Gaba
AHN-RCH-LHA Otp Fezf1 Glut
AHN-SBPV-PVHd Pdrm12 Gaba
APN C1ql2 Glut
APN C1ql4 Glut
ARH-PVi Six6 Dopa-Gaba
ARH-PVp Tbx3 Gaba
ARH-PVp Tbx3 Glut
AV Col27a1 Glut
AVPV-MEPO-SFO Tbr1 Glut
Astro-CB NN
Astro-NT NN
Astro-OLF NN
Astro-TE NN
Astroependymal NN
B-PB Nr4a2 Glut
BAM NN
BST Tac2 Gaba
BST-MPN Six3 Nrgn Gaba
BST-SI-AAA Six3 Slc22a3 Gaba
BST-po Iigp1 Glut
Bergmann NN
CA1-ProS Glut
CA2-FC-IG Glut
CA3 Glut
CB Granule Glut
CB PLI Gly-Gaba
CBN Dmbx1 Gaba
CBN Neurod2 Pvalb Glut
CBX Golgi Gly-Gaba
CBX MLI Cdh22 Gaba
CBX MLI Megf11 Gaba
CBX Purkinje Gaba
CEA-AAA-BST Six3 Sp9 Gaba
CEA-BST Ebf1 Pdyn Gaba
CEA-BST Gal Avp Gaba
CEA-BST Rai14 Pdyn Crh Gaba
CEA-BST Six3 Cyp26b1 Gaba
CHOR NN
CLA-EPd-CTX Car3 Glut
CM-IAD-CL-PCN Sema5b Glut
COAa-PAA-MEA Barhl2 Glut
COAp Grxcr2 Glut
CS-PRNr-DR En1 Sox2 Gaba
CS-PRNr-PCG Tmem163 Otp Gaba
CS-RPO Meis2 Gaba
CT SUB Glut
CU-ECU Pax2 Gly-Gaba
CU-ECU-SPVI Foxb1 Glut
CUN Evx2 Lhx2 Glut
CUN-PPN Evx2 Meis2 Glut
DC NN
DCO Il22 Gly-Gaba
DCO UBC Glut
DG Glut
DG-PIR Ex IMN
DMH Hmx2 Gaba
DMH Hmx2 Glut
DMH Nkx2-4 Glut
DMH Prdm13 Gaba
DMH-LHA Gsx1 Gaba
DMH-LHA Vgll2 Glut
DMX VII Tbx20 Chol
DTN-LDT-IPN Otp Pax3 Gaba
ENTmv-PA-COAp Glut
Endo NN
Ependymal NN
GPe-SI Sox6 Cyp26b1 Gaba
GPi Tbr1 Cngb3 Gaba-Glut
GRN-IRN-MDRNd Ikzf1 Gly-Gaba
HB Calcb Chol
HPF CR Glut
HY Gnrh1 Glut
Hypendymal NN
IA Mgp Gaba
IC Six3 En2 Gaba
IC Tfap2d Maf Glut
IF-RL-CLI-PAG Foxa1 Glut
IO Fgl2 Glut
IPN Otp Crisp1 Gaba
IPN-LDT Vsx2 Nkx6-1 Glut
IRN Dmbx1 Pax2 Gly-Gaba
IRN Vip Glut
IT AON-TT-DP Glut
IT EP-CLA Glut
L2 IT ENT-po Glut
L2 IT PPP-APr Glut
L2/3 IT CTX Glut
L2/3 IT ENT Glut
L2/3 IT PIR-ENTl Glut
L2/3 IT PPP Glut
L2/3 IT RSP Glut
L4 RSP-ACA Glut
L4/5 IT CTX Glut
L5 ET CTX Glut
L5 IT CTX Glut
L5 NP CTX Glut
L5 PPP Glut
L5/6 IT TPE-ENT Glut
L6 CT CTX Glut
L6 IT CTX Glut
L6b CTX Glut
L6b EPd Glut
L6b/CT ENT Glut
LA-BLA-BMA-PA Glut
LDT Fgf7 Gaba
LDT Vsx2 Nkx6-1 Nfib Glut
LDT-DTN Gata3 Nfix Gaba
LDT-PCG St18 Gaba
LDT-PCG Vsx2 Lhx4 Glut
LDT-PCG-CS Gata3 Lhx1 Gaba
LGv-SPFp-SPFm Nkx2-2 Tcf7l2 Gaba
LGv-ZI Otx2 Gaba
LH Pou4f1 Sox1 Glut
LHA Barhl2 Glut
LHA Pmch Glut
LHA-AHN-PVH Otp Trh Glut
LHA-MEA Otp Glut
LSX Nkx2-1 Gaba
LSX Otx2 Gaba
LSX Prdm12 Slit2 Gaba
LSX Prdm12 Zeb2 Gaba
LSX Sall3 Lmo1 Gaba
LSX Sall3 Pax6 Gaba
Lamp5 Gaba
Lamp5 Lhx6 Gaba
Lymphoid NN
MARN-GRN Pyy Glut
MARN-PPY Ngfr Gly-Gaba
MB-MY Tph2 Glut-Sero
MB-ant-ve Dmrta2 Glut
MDRN Hoxb5 Ebf2 Gly-Gaba
MDRNd Bves Glut
MDRNd Prox1 Pax6 Gly-Gaba
MDRNv Crp Glut
MDRNv Lhx4 Qrfprl Glut
MEA Otp Foxp2 Glut
MEA Slc17a7 Glut
MEA-BST Lhx6 Nfib Gaba
MEA-BST Lhx6 Nr2e1 Gaba
MEA-BST Lhx6 Sp9 Gaba
MEA-BST Otp Zic2 Glut
MEA-BST Sox6 Gaba
MEA-COA-BMA Ccdc42 Glut
MEV Ppp1r1c Glut
MG-POL-SGN Nts Glut
MH Tac2 Glut
MM Foxb1 Glut
MM-ant Foxb1 Glut
MPN-MPO-LPO Lhx6 Zfhx3 Gaba
MPN-MPO-PVpo Hmx2 Glut
MPO-ADP Lhx8 Gaba
MRN Pou3f1 C1ql4 Glut
MRN-PAG Nkx6-1 Glut
MRN-PPN-CUN Pax8 Gaba
MRN-VTN-PPN Pax5 Cdh23 Gaba
MS-SF Bsx Glut
MV Nkx6-1 Gly-Gaba
MV Nr4a2 Gly-Gaba
MV Pax6 Gly-Gaba
MV Xdh Gly-Gaba
MV-SPIV Phox2b Ebf3 Lbx1 Glut
MV-SPIV Slc6a2 Glut
MV-SPIV Zic4 Neurod2 Glut
MV-SPIV-PRP Dmbx1 Gly-Gaba
MY Lhx1 Gly-Gaba
MY Prox1 Lmo7 Gly-Gaba
Microglia NN
Monocytes NN
ND-INC Foxd2 Glut
NDB-SI-MA-STRv Lhx8 Gaba
NDB-SI-ant Prdm12 Gaba
NI-RPO Gata3 Nr4a2 Gaba
NLL Gata3 Gly-Gaba
NLL-SOC Spp1 Glut
NLL-po Pax7 Gaba
NLOT Rho Glut
NP PPP Glut
NP SUB Glut
NTS Aldh1a2 Glut
NTS Dbh Glut
NTS Mbnl3 Glut
NTS Phox2b Glut
NTS-MDRNd Prox1 Zic1 Gly-Gaba
NTS-PARN Neurod2 Gly-Gaba
OB Dopa-Gaba
OB Eomes Ms4a15 Glut
OB Meis2 Thsd7b Gaba
OB Trdn Gaba
OB-STR-CTX Inh IMN
OB-in Frmd7 Gaba
OB-mi Frmd7 Gaba
OB-out Frmd7 Gaba
OEC NN
OPC NN
OT D3 Folh1 Gaba
Oligo NN
PAG Pou4f1 Bnc2 Glut
PAG Pou4f1 Ebf2 Glut
PAG Pou4f2 Glut
PAG Pou4f2 Mesi2 Glut
PAG Pou4f3 Glut
PAG Tcf24 Glut
PAG Ucn Glut
PAG-MRN Pou3f1 Glut
PAG-MRN Rln3 Gaba
PAG-MRN Tfap2b Glut
PAG-MRN-RN Foxa2 Gaba
PAG-ND-PCG Onecut1 Gaba
PAG-PPN Pax5 Sox21 Gaba
PAG-RN Nkx2-2 Otx1 Gaba
PAG-SC Neurod2 Meis2 Glut
PAG-SC Pou4f1 Zic1 Glut
PAL-STR Gaba-Chol
PARN-MDRNd-NTS Gbx2 Gly-Gaba
PAS-MV Ebf2 Gly-Gaba
PB Evx2 Glut
PB Lmx1a Glut
PB Pax5 Glut
PB Sst Gly-Gaba
PB-NTS Phox2b Ebf3 Lmx1b Glut
PB-PSV Phox2b Glut
PB-SUT Tlx3 Lhx2 Glut
PBG Mtnr1a Glut-Chol
PCG-PRNr Vsx2 Nkx6-1 Glut
PDTg Otp Olig3 Gaba
PDTg Otp Shroom3 Gaba
PDTg-PCG Pax6 Gaba
PF Fzd5 Glut
PG-TRN-LRN Fat2 Glut
PGRN-PARN-MDRN Hoxb5 Glut
PGRNd Dmbx1 Glut
PH Pitx2 Glut
PH-LHA Foxb1 Glut
PH-SUM Foxa1 Glut
PH-an Pitx2 Glut
PH-ant-LHA Otp Bsx Glut
PMd-LHA Foxb1 Glut
PMv-TMv Pitx2 Glut
POR Gata3 Gly-Gaba
POR Spp1 Gly-Gaba
PPN-CUN-PCG Otp En1 Gaba
PPY-PGRNl Vip Glyc-Gaba
PRC-PAG Pax6 Glut
PRC-PAG Tcf7l2 Irx2 Glut
PRNc Otp Gly-Gaba
PRNc Prox1 Brs3 Gly-Gaba
PRNc-NI-SG-RPO Vsx2 Nr4a2 Glut
PRNc-PARN Tlx1 Glut
PRNr Otp Nfib Glut
PRP Gata3 Slc6a5 Gly-Gaba
PRP Otp Gly-Gaba
PRP-NI-PRNc-GRN Otp Glut
PRT Mecom Gaba
PRT Tcf7l2 Gaba
PSV Lmx1a Trpv6 Glut
PSV Pax2 Gly-Gaba
PSV Pvalb Lhx2 Glut
PVH-SO-PVa Otp Glut
PVHd Gsc Gaba
PVHd-DMH Lhx6 Gaba
PVHd-SBPV Six3 Prox1 Gaba
PVR Six3 Sox3 Gaba
PVT-PT Ntrk1 Glut
PVpo-VMPO-MPN Hmx2 Gaba
Peri NN
Pineal Crx Glut
Pvalb Gaba
Pvalb chandelier Gaba
RE-Xi Nox4 Glut
RHP-COA Ndnf Gaba
RN Spp1 Glut
RO-RPA Pkd2l1 Gaba
RPA Pax6 Hoxb5 Gly-Gaba
RT-ZI Gnb3 Gaba
SBPV-PVa Six6 Satb2 Gaba
SC Bnc2 Glut
SC Lef1 Otx2 Gaba
SC Otx2 Gcnt4 Gaba
SC Tnnt1 Gli3 Gaba
SC-PAG Lef1 Emx2 Gaba
SCH Six6 Cdc14a Gaba
SCdg-PAG Tfap2b Glut
SCig Foxb1 Glut
SCig Foxb1 Otx2 Glut
SCig Tfap2b Chrnb3 Glut
SCig-an-PPT Foxb1 Glut
SCiw Pitx2 Glut
SCm-PAG Cdh23 Gaba
SCop Pou4f2 Neurod2 Glut
SCop Sln Glut
SCs Dmbx1 Gaba
SCs Lef1 Gli3 Gaba
SCs Pax7 Nfia Gaba
SCsg Gabrr2 Gaba
SCsg Pde5a Glut
SI-MA-ACB Ebf1 Bnc2 Gaba
SI-MA-LPO-LHA Skor1 Glut
SI-MPO-LPO Lhx8 Gaba
SMC NN
SNc-VTA-RAmb Foxa1 Dopa
SNr Six3 Gaba
SNr-VTA Pax5 Npas1 Gaba
SPA-SPFm-SPFp-POL-PIL-PoT Sp9 Glut
SPVC Ccdc172 Glut
SPVC Mafa Glut
SPVC Nmu Glut
SPVI-SPVC Sall3 Lhx1 Gly-Gaba
SPVI-SPVC Sall3 Nfib Gly-Gaba
SPVI-SPVC Tlx3 Ebf3 Glut
SPVO Mafa Meis2 Glut
STN-PSTN Pitx2 Glut
STR D1 Gaba
STR D1 Sema5a Gaba
STR D2 Gaba
STR Lhx8 Gaba
STR Prox1 Lhx6 Gaba
STR-PAL Chst9 Gaba
SUB-ProS Glut
Sncg Gaba
Sst Chodl Gaba
Sst Gaba
TH Prkcd Grin2c Glut
TMd-DMH Foxd2 Gaba
TMv-PMv Tbx3 Hist-Gaba
TRS-BAC Sln Glut
TU-ARH Otp Six6 Gaba
Tanycyte NN
VCO Mafa Meis2 Glut
VLMC NN
VMH Fezf1 Glut
VMH Nr5a1 Glut
Vip Gaba
ZI Pax6 Gaba
Generating one-hot matrix...
100%|██████████| 23/23 [00:08<00:00, 2.56it/s]
All done!
Dataset comprises 338 cell types.
Calculating cell type distribution for microenvironments...
Microenvironments comprise 42.09 cells/spots on average.
Minimum: 20, Maximum: 102
Perform overclustered initialization on the cell type distributions of cell neighborhoods. Resulting in Qmax niches. The distributions of niches are also computed. We start from 100 niches on this dataset due to its large capacity and we anticipate to find fine structures of the brain.
[7]:
model.initialize_clusters(dim_reduction=True, # default
explained_var=None, # default
n_components=None, # default
n_components_max=100, # default
standardize=True, # default
method='kmeans', # default
Qmax=100,
)
Performing dimension reduction...
Returning 100 principal components.
Initializing niches...
100 initial niches defined.
Perform hierarchical distribution matching to reduce the niche number to no less than Qmin. This step results in niche assignment under a sequence of different niche numbers (usually from Qmax to Qmin).
[ ]:
model.hier_dist_match(assign_metric='jsd', # default
weighted_merge=True, # default
max_iters=100, # default
tol=1e-4, # default
test_kmeans=False, # default
Qmin=2, # default
)
Starting from 100 cell niches...
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99]
56%|█████▌ | 56/100 [15:58<12:33, 17.12s/it]
Distribution of cell niches (centers) converge at iteration 57.
100 cell niches left.
Merging cell niche 64 and cell niche 98...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
1%| | 1/100 [00:32<53:43, 32.56s/it]
Distribution of cell niches (centers) converge at iteration 2.
99 cell niches left.
Merging cell niche 67 and cell niche 85...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
2%|▏ | 2/100 [00:49<40:28, 24.78s/it]
Distribution of cell niches (centers) converge at iteration 3.
98 cell niches left.
Merging cell niche 67 and cell niche 43...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
1%| | 1/100 [00:33<54:32, 33.06s/it]
Distribution of cell niches (centers) converge at iteration 2.
97 cell niches left.
Merging cell niche 64 and cell niche 53...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 44, 45, 46, 47, 48, 49, 50, 51, 52, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
2%|▏ | 2/100 [00:51<42:07, 25.79s/it]
Distribution of cell niches (centers) converge at iteration 3.
96 cell niches left.
Merging cell niche 64 and cell niche 41...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 50, 51, 52, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
4%|▍ | 4/100 [01:22<33:11, 20.75s/it]
Distribution of cell niches (centers) converge at iteration 5.
95 cell niches left.
Merging cell niche 64 and cell niche 70...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 50, 51, 52, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
1%| | 1/100 [00:31<52:36, 31.88s/it]
Distribution of cell niches (centers) converge at iteration 2.
94 cell niches left.
Merging cell niche 33 and cell niche 50...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 51, 52, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
27%|██▋ | 27/100 [07:38<20:39, 16.98s/it]
Distribution of cell niches (centers) converge at iteration 28.
93 cell niches left.
Merging cell niche 93 and cell niche 33...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 51, 52, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
1%| | 1/100 [00:33<55:53, 33.87s/it]
Distribution of cell niches (centers) converge at iteration 2.
92 cell niches left.
Merging cell niche 23 and cell niche 64...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 51, 52, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 65, 66, 67, 68, 69, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
2%|▏ | 2/100 [00:48<39:32, 24.21s/it]
Distribution of cell niches (centers) converge at iteration 3.
91 cell niches left.
Merging cell niche 77 and cell niche 23...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 51, 52, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 65, 66, 67, 68, 69, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
1%| | 1/100 [00:31<52:39, 31.91s/it]
Distribution of cell niches (centers) converge at iteration 2.
90 cell niches left.
Merging cell niche 9 and cell niche 67...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 51, 52, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 65, 66, 68, 69, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
2%|▏ | 2/100 [00:46<38:17, 23.45s/it]
Distribution of cell niches (centers) converge at iteration 3.
89 cell niches left.
Merging cell niche 65 and cell niche 77...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 51, 52, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 65, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
2%|▏ | 2/100 [00:47<38:48, 23.76s/it]
Distribution of cell niches (centers) converge at iteration 3.
88 cell niches left.
Merging cell niche 61 and cell niche 65...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 51, 52, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
3%|▎ | 3/100 [01:03<33:57, 21.01s/it]
Distribution of cell niches (centers) converge at iteration 4.
87 cell niches left.
Merging cell niche 21 and cell niche 61...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 51, 52, 54, 55, 56, 57, 58, 59, 60, 62, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
1%| | 1/100 [00:32<52:57, 32.10s/it]
Distribution of cell niches (centers) converge at iteration 2.
86 cell niches left.
Merging cell niche 93 and cell niche 60...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 51, 52, 54, 55, 56, 57, 58, 59, 62, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
1%| | 1/100 [00:31<51:14, 31.06s/it]
Distribution of cell niches (centers) converge at iteration 2.
85 cell niches left.
Merging cell niche 48 and cell niche 21...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 51, 52, 54, 55, 56, 57, 58, 59, 62, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
1%| | 1/100 [00:30<50:14, 30.45s/it]
Distribution of cell niches (centers) converge at iteration 2.
84 cell niches left.
Merging cell niche 42 and cell niche 5...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 51, 52, 54, 55, 56, 57, 58, 59, 62, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 99]
4%|▍ | 4/100 [01:16<30:40, 19.17s/it]
Distribution of cell niches (centers) converge at iteration 5.
83 cell niches left.
Merging cell niche 48 and cell niche 91...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 51, 52, 54, 55, 56, 57, 58, 59, 62, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 82, 83, 84, 86, 87, 88, 89, 90, 92, 93, 94, 95, 96, 97, 99]
1%| | 1/100 [00:31<51:16, 31.08s/it]
Distribution of cell niches (centers) converge at iteration 2.
82 cell niches left.
Merging cell niche 9 and cell niche 88...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 51, 52, 54, 55, 56, 57, 58, 59, 62, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 82, 83, 84, 86, 87, 89, 90, 92, 93, 94, 95, 96, 97, 99]
13%|█▎ | 13/100 [03:34<23:54, 16.49s/it]
Distribution of cell niches (centers) converge at iteration 14.
81 cell niches left.
Merging cell niche 9 and cell niche 97...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 48, 49, 51, 52, 54, 55, 56, 57, 58, 59, 62, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 82, 83, 84, 86, 87, 89, 90, 92, 93, 94, 95, 96, 99]
1%| | 1/100 [00:30<50:02, 30.33s/it]
Distribution of cell niches (centers) converge at iteration 2.
80 cell niches left.
Merging cell niche 80 and cell niche 48...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 49, 51, 52, 54, 55, 56, 57, 58, 59, 62, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 82, 83, 84, 86, 87, 89, 90, 92, 93, 94, 95, 96, 99]
1%| | 1/100 [00:29<49:20, 29.90s/it]
Distribution of cell niches (centers) converge at iteration 2.
79 cell niches left.
Merging cell niche 80 and cell niche 82...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 49, 51, 52, 54, 55, 56, 57, 58, 59, 62, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 83, 84, 86, 87, 89, 90, 92, 93, 94, 95, 96, 99]
2%|▏ | 2/100 [00:44<36:21, 22.26s/it]
Distribution of cell niches (centers) converge at iteration 3.
78 cell niches left.
Merging cell niche 80 and cell niche 16...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 17, 18, 19, 20, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37, 38, 39, 40, 42, 44, 45, 46, 47, 49, 51, 52, 54, 55, 56, 57, 58, 59, 62, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 83, 84, 86, 87, 89, 90, 92, 93, 94, 95, 96, 99]
2%|▏ | 2/100 [00:44<36:29, 22.35s/it]
Distribution of cell niches (centers) converge at iteration 3.
77 cell niches left.
Merging cell niche 80 and cell niche 37...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 17, 18, 19, 20, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 38, 39, 40, 42, 44, 45, 46, 47, 49, 51, 52, 54, 55, 56, 57, 58, 59, 62, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 83, 84, 86, 87, 89, 90, 92, 93, 94, 95, 96, 99]
2%|▏ | 2/100 [00:43<35:31, 21.75s/it]
Distribution of cell niches (centers) converge at iteration 3.
76 cell niches left.
Merging cell niche 80 and cell niche 92...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 17, 18, 19, 20, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 38, 39, 40, 42, 44, 45, 46, 47, 49, 51, 52, 54, 55, 56, 57, 58, 59, 62, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 83, 84, 86, 87, 89, 90, 93, 94, 95, 96, 99]
4%|▍ | 4/100 [01:13<29:17, 18.30s/it]
Distribution of cell niches (centers) converge at iteration 5.
75 cell niches left.
Merging cell niche 35 and cell niche 9...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 10, 11, 12, 13, 14, 15, 17, 18, 19, 20, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 38, 39, 40, 42, 44, 45, 46, 47, 49, 51, 52, 54, 55, 56, 57, 58, 59, 62, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 83, 84, 86, 87, 89, 90, 93, 94, 95, 96, 99]
9%|▉ | 9/100 [02:26<24:37, 16.23s/it]
Distribution of cell niches (centers) converge at iteration 10.
74 cell niches left.
Merging cell niche 35 and cell niche 62...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 10, 11, 12, 13, 14, 15, 17, 18, 19, 20, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 38, 39, 40, 42, 44, 45, 46, 47, 49, 51, 52, 54, 55, 56, 57, 58, 59, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 83, 84, 86, 87, 89, 90, 93, 94, 95, 96, 99]
8%|▊ | 8/100 [02:08<24:42, 16.11s/it]
Distribution of cell niches (centers) converge at iteration 9.
73 cell niches left.
Merging cell niche 80 and cell niche 19...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 10, 11, 12, 13, 14, 15, 17, 18, 20, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32, 34, 35, 36, 38, 39, 40, 42, 44, 45, 46, 47, 49, 51, 52, 54, 55, 56, 57, 58, 59, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 83, 84, 86, 87, 89, 90, 93, 94, 95, 96, 99]
10%|█ | 10/100 [02:36<23:28, 15.65s/it]
Distribution of cell niches (centers) converge at iteration 11.
72 cell niches left.
Merging cell niche 93 and cell niche 26...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 10, 11, 12, 13, 14, 15, 17, 18, 20, 22, 24, 25, 27, 28, 29, 30, 31, 32, 34, 35, 36, 38, 39, 40, 42, 44, 45, 46, 47, 49, 51, 52, 54, 55, 56, 57, 58, 59, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 80, 81, 83, 84, 86, 87, 89, 90, 93, 94, 95, 96, 99]
2%|▏ | 2/100 [00:42<34:18, 21.00s/it]
Distribution of cell niches (centers) converge at iteration 3.
71 cell niches left.
Merging cell niche 56 and cell niche 80...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 10, 11, 12, 13, 14, 15, 17, 18, 20, 22, 24, 25, 27, 28, 29, 30, 31, 32, 34, 35, 36, 38, 39, 40, 42, 44, 45, 46, 47, 49, 51, 52, 54, 55, 56, 57, 58, 59, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 81, 83, 84, 86, 87, 89, 90, 93, 94, 95, 96, 99]
5%|▌ | 5/100 [01:23<26:30, 16.74s/it]
Distribution of cell niches (centers) converge at iteration 6.
70 cell niches left.
Merging cell niche 55 and cell niche 56...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 10, 11, 12, 13, 14, 15, 17, 18, 20, 22, 24, 25, 27, 28, 29, 30, 31, 32, 34, 35, 36, 38, 39, 40, 42, 44, 45, 46, 47, 49, 51, 52, 54, 55, 57, 58, 59, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 81, 83, 84, 86, 87, 89, 90, 93, 94, 95, 96, 99]
2%|▏ | 2/100 [00:41<33:44, 20.65s/it]
Distribution of cell niches (centers) converge at iteration 3.
69 cell niches left.
Merging cell niche 69 and cell niche 42...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 10, 11, 12, 13, 14, 15, 17, 18, 20, 22, 24, 25, 27, 28, 29, 30, 31, 32, 34, 35, 36, 38, 39, 40, 44, 45, 46, 47, 49, 51, 52, 54, 55, 57, 58, 59, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 81, 83, 84, 86, 87, 89, 90, 93, 94, 95, 96, 99]
2%|▏ | 2/100 [00:41<33:35, 20.56s/it]
Distribution of cell niches (centers) converge at iteration 3.
68 cell niches left.
Merging cell niche 63 and cell niche 55...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 10, 11, 12, 13, 14, 15, 17, 18, 20, 22, 24, 25, 27, 28, 29, 30, 31, 32, 34, 35, 36, 38, 39, 40, 44, 45, 46, 47, 49, 51, 52, 54, 57, 58, 59, 63, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 81, 83, 84, 86, 87, 89, 90, 93, 94, 95, 96, 99]
11%|█ | 11/100 [02:47<22:37, 15.26s/it]
Distribution of cell niches (centers) converge at iteration 12.
67 cell niches left.
Merging cell niche 95 and cell niche 63...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 10, 11, 12, 13, 14, 15, 17, 18, 20, 22, 24, 25, 27, 28, 29, 30, 31, 32, 34, 35, 36, 38, 39, 40, 44, 45, 46, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 81, 83, 84, 86, 87, 89, 90, 93, 94, 95, 96, 99]
21%|██ | 21/100 [05:04<19:05, 14.49s/it]
Distribution of cell niches (centers) converge at iteration 22.
66 cell niches left.
Merging cell niche 2 and cell niche 24...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 2, 3, 4, 6, 7, 8, 10, 11, 12, 13, 14, 15, 17, 18, 20, 22, 25, 27, 28, 29, 30, 31, 32, 34, 35, 36, 38, 39, 40, 44, 45, 46, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 81, 83, 84, 86, 87, 89, 90, 93, 94, 95, 96, 99]
3%|▎ | 3/100 [00:55<29:58, 18.54s/it]
Distribution of cell niches (centers) converge at iteration 4.
65 cell niches left.
Merging cell niche 46 and cell niche 2...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 8, 10, 11, 12, 13, 14, 15, 17, 18, 20, 22, 25, 27, 28, 29, 30, 31, 32, 34, 35, 36, 38, 39, 40, 44, 45, 46, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 81, 83, 84, 86, 87, 89, 90, 93, 94, 95, 96, 99]
1%| | 1/100 [00:26<44:22, 26.90s/it]
Distribution of cell niches (centers) converge at iteration 2.
64 cell niches left.
Merging cell niche 7 and cell niche 35...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 8, 10, 11, 12, 13, 14, 15, 17, 18, 20, 22, 25, 27, 28, 29, 30, 31, 32, 34, 36, 38, 39, 40, 44, 45, 46, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 81, 83, 84, 86, 87, 89, 90, 93, 94, 95, 96, 99]
4%|▍ | 4/100 [01:06<26:27, 16.53s/it]
Distribution of cell niches (centers) converge at iteration 5.
63 cell niches left.
Merging cell niche 7 and cell niche 17...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 8, 10, 11, 12, 13, 14, 15, 18, 20, 22, 25, 27, 28, 29, 30, 31, 32, 34, 36, 38, 39, 40, 44, 45, 46, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 75, 76, 78, 79, 81, 83, 84, 86, 87, 89, 90, 93, 94, 95, 96, 99]
4%|▍ | 4/100 [01:05<26:09, 16.34s/it]
Distribution of cell niches (centers) converge at iteration 5.
62 cell niches left.
Merging cell niche 46 and cell niche 78...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 8, 10, 11, 12, 13, 14, 15, 18, 20, 22, 25, 27, 28, 29, 30, 31, 32, 34, 36, 38, 39, 40, 44, 45, 46, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 75, 76, 79, 81, 83, 84, 86, 87, 89, 90, 93, 94, 95, 96, 99]
10%|█ | 10/100 [02:27<22:06, 14.74s/it]
Distribution of cell niches (centers) converge at iteration 11.
61 cell niches left.
Merging cell niche 93 and cell niche 84...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 8, 10, 11, 12, 13, 14, 15, 18, 20, 22, 25, 27, 28, 29, 30, 31, 32, 34, 36, 38, 39, 40, 44, 45, 46, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 75, 76, 79, 81, 83, 86, 87, 89, 90, 93, 94, 95, 96, 99]
3%|▎ | 3/100 [00:53<28:34, 17.67s/it]
Distribution of cell niches (centers) converge at iteration 4.
60 cell niches left.
Merging cell niche 93 and cell niche 14...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 8, 10, 11, 12, 13, 15, 18, 20, 22, 25, 27, 28, 29, 30, 31, 32, 34, 36, 38, 39, 40, 44, 45, 46, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 75, 76, 79, 81, 83, 86, 87, 89, 90, 93, 94, 95, 96, 99]
12%|█▏ | 12/100 [02:50<20:50, 14.21s/it]
Distribution of cell niches (centers) converge at iteration 13.
59 cell niches left.
Merging cell niche 75 and cell niche 31...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 8, 10, 11, 12, 13, 15, 18, 20, 22, 25, 27, 28, 29, 30, 32, 34, 36, 38, 39, 40, 44, 45, 46, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 75, 76, 79, 81, 83, 86, 87, 89, 90, 93, 94, 95, 96, 99]
9%|▉ | 9/100 [02:09<21:47, 14.36s/it]
Distribution of cell niches (centers) converge at iteration 10.
58 cell niches left.
Merging cell niche 46 and cell niche 36...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 8, 10, 11, 12, 13, 15, 18, 20, 22, 25, 27, 28, 29, 30, 32, 34, 38, 39, 40, 44, 45, 46, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 75, 76, 79, 81, 83, 86, 87, 89, 90, 93, 94, 95, 96, 99]
1%| | 1/100 [00:24<41:06, 24.92s/it]
Distribution of cell niches (centers) converge at iteration 2.
57 cell niches left.
Merging cell niche 8 and cell niche 46...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 8, 10, 11, 12, 13, 15, 18, 20, 22, 25, 27, 28, 29, 30, 32, 34, 38, 39, 40, 44, 45, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 75, 76, 79, 81, 83, 86, 87, 89, 90, 93, 94, 95, 96, 99]
12%|█▏ | 12/100 [02:45<20:17, 13.83s/it]
Distribution of cell niches (centers) converge at iteration 13.
56 cell niches left.
Merging cell niche 93 and cell niche 75...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 8, 10, 11, 12, 13, 15, 18, 20, 22, 25, 27, 28, 29, 30, 32, 34, 38, 39, 40, 44, 45, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 76, 79, 81, 83, 86, 87, 89, 90, 93, 94, 95, 96, 99]
1%| | 1/100 [00:26<42:57, 26.03s/it]
Distribution of cell niches (centers) converge at iteration 2.
55 cell niches left.
Merging cell niche 68 and cell niche 8...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 11, 12, 13, 15, 18, 20, 22, 25, 27, 28, 29, 30, 32, 34, 38, 39, 40, 44, 45, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 76, 79, 81, 83, 86, 87, 89, 90, 93, 94, 95, 96, 99]
23%|██▎ | 23/100 [05:02<16:54, 13.17s/it]
Distribution of cell niches (centers) converge at iteration 24.
54 cell niches left.
Merging cell niche 68 and cell niche 87...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 11, 12, 13, 15, 18, 20, 22, 25, 27, 28, 29, 30, 32, 34, 38, 39, 40, 44, 45, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 76, 79, 81, 83, 86, 89, 90, 93, 94, 95, 96, 99]
12%|█▏ | 12/100 [02:41<19:44, 13.47s/it]
Distribution of cell niches (centers) converge at iteration 13.
53 cell niches left.
Merging cell niche 27 and cell niche 89...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 11, 12, 13, 15, 18, 20, 22, 25, 27, 28, 29, 30, 32, 34, 38, 39, 40, 44, 45, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 76, 79, 81, 83, 86, 90, 93, 94, 95, 96, 99]
2%|▏ | 2/100 [00:35<29:19, 17.95s/it]
Distribution of cell niches (centers) converge at iteration 3.
52 cell niches left.
Merging cell niche 27 and cell niche 95...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 11, 12, 13, 15, 18, 20, 22, 25, 27, 28, 29, 30, 32, 34, 38, 39, 40, 44, 45, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 76, 79, 81, 83, 86, 90, 93, 94, 96, 99]
61%|██████ | 61/100 [12:23<07:55, 12.19s/it]
Distribution of cell niches (centers) converge at iteration 62.
51 cell niches left.
Merging cell niche 52 and cell niche 15...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 11, 12, 13, 18, 20, 22, 25, 27, 28, 29, 30, 32, 34, 38, 39, 40, 44, 45, 47, 49, 51, 52, 54, 57, 58, 59, 66, 68, 69, 71, 72, 73, 74, 76, 79, 81, 83, 86, 90, 93, 94, 96, 99]
26%|██▌ | 26/100 [05:20<15:10, 12.31s/it]
Distribution of cell niches (centers) converge at iteration 27.
50 cell niches left.
Merging cell niche 52 and cell niche 68...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 11, 12, 13, 18, 20, 22, 25, 27, 28, 29, 30, 32, 34, 38, 39, 40, 44, 45, 47, 49, 51, 52, 54, 57, 58, 59, 66, 69, 71, 72, 73, 74, 76, 79, 81, 83, 86, 90, 93, 94, 96, 99]
1%| | 1/100 [00:24<41:05, 24.90s/it]
Distribution of cell niches (centers) converge at iteration 2.
49 cell niches left.
Merging cell niche 27 and cell niche 93...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 11, 12, 13, 18, 20, 22, 25, 27, 28, 29, 30, 32, 34, 38, 39, 40, 44, 45, 47, 49, 51, 52, 54, 57, 58, 59, 66, 69, 71, 72, 73, 74, 76, 79, 81, 83, 86, 90, 94, 96, 99]
1%| | 1/100 [00:22<37:38, 22.82s/it]
Distribution of cell niches (centers) converge at iteration 2.
48 cell niches left.
Merging cell niche 74 and cell niche 52...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 11, 12, 13, 18, 20, 22, 25, 27, 28, 29, 30, 32, 34, 38, 39, 40, 44, 45, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 74, 76, 79, 81, 83, 86, 90, 94, 96, 99]
47%|████▋ | 47/100 [09:14<10:24, 11.79s/it]
Distribution of cell niches (centers) converge at iteration 48.
47 cell niches left.
Merging cell niche 40 and cell niche 81...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 11, 12, 13, 18, 20, 22, 25, 27, 28, 29, 30, 32, 34, 38, 39, 40, 44, 45, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 74, 76, 79, 83, 86, 90, 94, 96, 99]
19%|█▉ | 19/100 [03:46<16:04, 11.91s/it]
Distribution of cell niches (centers) converge at iteration 20.
46 cell niches left.
Merging cell niche 40 and cell niche 74...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 11, 12, 13, 18, 20, 22, 25, 27, 28, 29, 30, 32, 34, 38, 39, 40, 44, 45, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 76, 79, 83, 86, 90, 94, 96, 99]
2%|▏ | 2/100 [00:33<27:17, 16.71s/it]
Distribution of cell niches (centers) converge at iteration 3.
45 cell niches left.
Merging cell niche 76 and cell niche 40...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 11, 12, 13, 18, 20, 22, 25, 27, 28, 29, 30, 32, 34, 38, 39, 44, 45, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 76, 79, 83, 86, 90, 94, 96, 99]
61%|██████ | 61/100 [11:35<07:24, 11.41s/it]
Distribution of cell niches (centers) converge at iteration 62.
44 cell niches left.
Merging cell niche 83 and cell niche 34...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 11, 12, 13, 18, 20, 22, 25, 27, 28, 29, 30, 32, 38, 39, 44, 45, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 76, 79, 83, 86, 90, 94, 96, 99]
28%|██▊ | 28/100 [05:18<13:39, 11.38s/it]
Distribution of cell niches (centers) converge at iteration 29.
43 cell niches left.
Merging cell niche 83 and cell niche 30...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 11, 12, 13, 18, 20, 22, 25, 27, 28, 29, 32, 38, 39, 44, 45, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 76, 79, 83, 86, 90, 94, 96, 99]
4%|▍ | 4/100 [00:53<21:22, 13.36s/it]
Distribution of cell niches (centers) converge at iteration 5.
42 cell niches left.
Merging cell niche 83 and cell niche 11...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 12, 13, 18, 20, 22, 25, 27, 28, 29, 32, 38, 39, 44, 45, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 76, 79, 83, 86, 90, 94, 96, 99]
2%|▏ | 2/100 [00:32<26:20, 16.13s/it]
Distribution of cell niches (centers) converge at iteration 3.
41 cell niches left.
Merging cell niche 90 and cell niche 83...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 12, 13, 18, 20, 22, 25, 27, 28, 29, 32, 38, 39, 44, 45, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 76, 79, 86, 90, 94, 96, 99]
4%|▍ | 4/100 [00:54<21:53, 13.68s/it]
Distribution of cell niches (centers) converge at iteration 5.
40 cell niches left.
Merging cell niche 90 and cell niche 99...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 12, 13, 18, 20, 22, 25, 27, 28, 29, 32, 38, 39, 44, 45, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 76, 79, 86, 90, 94, 96]
23%|██▎ | 23/100 [04:13<14:10, 11.04s/it]
Distribution of cell niches (centers) converge at iteration 24.
39 cell niches left.
Merging cell niche 90 and cell niche 44...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 12, 13, 18, 20, 22, 25, 27, 28, 29, 32, 38, 39, 45, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 76, 79, 86, 90, 94, 96]
1%| | 1/100 [00:21<35:07, 21.28s/it]
Distribution of cell niches (centers) converge at iteration 2.
38 cell niches left.
Merging cell niche 90 and cell niche 39...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 12, 13, 18, 20, 22, 25, 27, 28, 29, 32, 38, 45, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 76, 79, 86, 90, 94, 96]
1%| | 1/100 [00:20<34:38, 21.00s/it]
Distribution of cell niches (centers) converge at iteration 2.
37 cell niches left.
Merging cell niche 76 and cell niche 90...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 12, 13, 18, 20, 22, 25, 27, 28, 29, 32, 38, 45, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 76, 79, 86, 94, 96]
2%|▏ | 2/100 [00:30<24:36, 15.06s/it]
Distribution of cell niches (centers) converge at iteration 3.
36 cell niches left.
Merging cell niche 76 and cell niche 45...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 4, 6, 7, 10, 12, 13, 18, 20, 22, 25, 27, 28, 29, 32, 38, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 76, 79, 86, 94, 96]
45%|████▌ | 45/100 [07:41<09:24, 10.27s/it]
Distribution of cell niches (centers) converge at iteration 46.
35 cell niches left.
Merging cell niche 25 and cell niche 4...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 6, 7, 10, 12, 13, 18, 20, 22, 25, 27, 28, 29, 32, 38, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 76, 79, 86, 94, 96]
2%|▏ | 2/100 [00:29<23:59, 14.69s/it]
Distribution of cell niches (centers) converge at iteration 3.
34 cell niches left.
Merging cell niche 25 and cell niche 22...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 6, 7, 10, 12, 13, 18, 20, 25, 27, 28, 29, 32, 38, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 76, 79, 86, 94, 96]
2%|▏ | 2/100 [00:29<23:44, 14.53s/it]
Distribution of cell niches (centers) converge at iteration 3.
33 cell niches left.
Merging cell niche 25 and cell niche 76...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 6, 7, 10, 12, 13, 18, 20, 25, 27, 28, 29, 32, 38, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 72, 73, 79, 86, 94, 96]
4%|▍ | 4/100 [00:49<19:42, 12.32s/it]
Distribution of cell niches (centers) converge at iteration 5.
32 cell niches left.
Merging cell niche 25 and cell niche 72...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 6, 7, 10, 12, 13, 18, 20, 25, 27, 28, 29, 32, 38, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 73, 79, 86, 94, 96]
29%|██▉ | 29/100 [05:00<12:14, 10.35s/it]
Distribution of cell niches (centers) converge at iteration 30.
31 cell niches left.
Merging cell niche 25 and cell niche 12...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 6, 7, 10, 13, 18, 20, 25, 27, 28, 29, 32, 38, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 73, 79, 86, 94, 96]
2%|▏ | 2/100 [00:28<23:30, 14.39s/it]
Distribution of cell niches (centers) converge at iteration 3.
30 cell niches left.
Merging cell niche 47 and cell niche 25...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 6, 7, 10, 13, 18, 20, 27, 28, 29, 32, 38, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 73, 79, 86, 94, 96]
2%|▏ | 2/100 [00:28<23:01, 14.10s/it]
Distribution of cell niches (centers) converge at iteration 3.
29 cell niches left.
Merging cell niche 47 and cell niche 96...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 6, 7, 10, 13, 18, 20, 27, 28, 29, 32, 38, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 73, 79, 86, 94]
40%|████ | 40/100 [06:20<09:30, 9.51s/it]
Distribution of cell niches (centers) converge at iteration 41.
28 cell niches left.
Merging cell niche 47 and cell niche 27...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 6, 7, 10, 13, 18, 20, 28, 29, 32, 38, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 73, 79, 86, 94]
1%| | 1/100 [00:17<29:23, 17.81s/it]
Distribution of cell niches (centers) converge at iteration 2.
27 cell niches left.
Merging cell niche 47 and cell niche 79...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 6, 7, 10, 13, 18, 20, 28, 29, 32, 38, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 73, 86, 94]
44%|████▍ | 44/100 [06:44<08:35, 9.20s/it]
Distribution of cell niches (centers) converge at iteration 45.
26 cell niches left.
Merging cell niche 47 and cell niche 32...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 6, 7, 10, 13, 18, 20, 28, 29, 38, 47, 49, 51, 54, 57, 58, 59, 66, 69, 71, 73, 86, 94]
3%|▎ | 3/100 [00:36<19:35, 12.12s/it]
Distribution of cell niches (centers) converge at iteration 4.
25 cell niches left.
Merging cell niche 47 and cell niche 54...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 6, 7, 10, 13, 18, 20, 28, 29, 38, 47, 49, 51, 57, 58, 59, 66, 69, 71, 73, 86, 94]
13%|█▎ | 13/100 [02:04<13:52, 9.57s/it]
Distribution of cell niches (centers) converge at iteration 14.
24 cell niches left.
Merging cell niche 51 and cell niche 38...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 6, 7, 10, 13, 18, 20, 28, 29, 47, 49, 51, 57, 58, 59, 66, 69, 71, 73, 86, 94]
4%|▍ | 4/100 [00:42<16:59, 10.62s/it]
Distribution of cell niches (centers) converge at iteration 5.
23 cell niches left.
Merging cell niche 51 and cell niche 94...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 6, 7, 10, 13, 18, 20, 28, 29, 47, 49, 51, 57, 58, 59, 66, 69, 71, 73, 86]
2%|▏ | 2/100 [00:25<20:28, 12.53s/it]
Distribution of cell niches (centers) converge at iteration 3.
22 cell niches left.
Merging cell niche 51 and cell niche 86...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 6, 7, 10, 13, 18, 20, 28, 29, 47, 49, 51, 57, 58, 59, 66, 69, 71, 73]
1%| | 1/100 [00:16<27:31, 16.69s/it]
Distribution of cell niches (centers) converge at iteration 2.
21 cell niches left.
Merging cell niche 57 and cell niche 51...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 6, 7, 10, 13, 18, 20, 28, 29, 47, 49, 57, 58, 59, 66, 69, 71, 73]
4%|▍ | 4/100 [00:40<16:18, 10.19s/it]
Distribution of cell niches (centers) converge at iteration 5.
20 cell niches left.
Merging cell niche 57 and cell niche 6...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 7, 10, 13, 18, 20, 28, 29, 47, 49, 57, 58, 59, 66, 69, 71, 73]
30%|███ | 30/100 [04:12<09:48, 8.41s/it]
Distribution of cell niches (centers) converge at iteration 31.
19 cell niches left.
Merging cell niche 47 and cell niche 66...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 7, 10, 13, 18, 20, 28, 29, 47, 49, 57, 58, 59, 69, 71, 73]
1%| | 1/100 [00:15<26:23, 15.99s/it]
Distribution of cell niches (centers) converge at iteration 2.
18 cell niches left.
Merging cell niche 58 and cell niche 47...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 7, 10, 13, 18, 20, 28, 29, 49, 57, 58, 59, 69, 71, 73]
4%|▍ | 4/100 [00:41<16:45, 10.48s/it]
Distribution of cell niches (centers) converge at iteration 5.
17 cell niches left.
Merging cell niche 57 and cell niche 58...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 7, 10, 13, 18, 20, 28, 29, 49, 57, 59, 69, 71, 73]
1%| | 1/100 [00:16<26:40, 16.16s/it]
Distribution of cell niches (centers) converge at iteration 2.
16 cell niches left.
Merging cell niche 29 and cell niche 57...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 7, 10, 13, 18, 20, 28, 29, 49, 59, 69, 71, 73]
5%|▌ | 5/100 [00:45<14:18, 9.04s/it]
Distribution of cell niches (centers) converge at iteration 6.
15 cell niches left.
Merging cell niche 29 and cell niche 28...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 7, 10, 13, 18, 20, 29, 49, 59, 69, 71, 73]
37%|███▋ | 37/100 [04:46<08:08, 7.75s/it]
Distribution of cell niches (centers) converge at iteration 38.
14 cell niches left.
Merging cell niche 29 and cell niche 73...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 7, 10, 13, 18, 20, 29, 49, 59, 69, 71]
2%|▏ | 2/100 [00:21<17:27, 10.69s/it]
Distribution of cell niches (centers) converge at iteration 3.
13 cell niches left.
Merging cell niche 10 and cell niche 18...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 7, 10, 13, 20, 29, 49, 59, 69, 71]
8%|▊ | 8/100 [01:05<12:30, 8.16s/it]
Distribution of cell niches (centers) converge at iteration 9.
12 cell niches left.
Merging cell niche 49 and cell niche 29...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 7, 10, 13, 20, 49, 59, 69, 71]
8%|▊ | 8/100 [01:04<12:16, 8.00s/it]
Distribution of cell niches (centers) converge at iteration 9.
11 cell niches left.
Merging cell niche 3 and cell niche 49...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 7, 10, 13, 20, 59, 69, 71]
14%|█▍ | 14/100 [01:45<10:45, 7.51s/it]
Distribution of cell niches (centers) converge at iteration 15.
10 cell niches left.
Merging cell niche 1 and cell niche 10...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 7, 13, 20, 59, 69, 71]
6%|▌ | 6/100 [00:47<12:21, 7.89s/it]
Distribution of cell niches (centers) converge at iteration 7.
9 cell niches left.
Merging cell niche 3 and cell niche 59...
Done!
Assigning cells to cell niche...
Current state: [0, 1, 3, 7, 13, 20, 69, 71]
4%|▍ | 4/100 [00:32<13:04, 8.17s/it]
Distribution of cell niches (centers) converge at iteration 5.
8 cell niches left.
Merging cell niche 3 and cell niche 0...
Done!
Assigning cells to cell niche...
Current state: [1, 3, 7, 13, 20, 69, 71]
3%|▎ | 3/100 [00:25<13:48, 8.55s/it]
Distribution of cell niches (centers) converge at iteration 4.
7 cell niches left.
Merging cell niche 13 and cell niche 3...
Done!
Assigning cells to cell niche...
Current state: [1, 7, 13, 20, 69, 71]
7%|▋ | 7/100 [00:50<11:11, 7.23s/it]
Distribution of cell niches (centers) converge at iteration 8.
6 cell niches left.
Merging cell niche 69 and cell niche 13...
Done!
Assigning cells to cell niche...
Current state: [1, 7, 20, 69, 71]
3%|▎ | 3/100 [00:24<13:11, 8.16s/it]
Distribution of cell niches (centers) converge at iteration 4.
5 cell niches left.
Merging cell niche 71 and cell niche 1...
Done!
Assigning cells to cell niche...
Current state: [7, 20, 69, 71]
1%| | 1/100 [00:11<19:29, 11.82s/it]
Distribution of cell niches (centers) converge at iteration 2.
4 cell niches left.
Merging cell niche 71 and cell niche 69...
Done!
Assigning cells to cell niche...
Current state: [7, 20, 71]
1%| | 1/100 [00:11<19:12, 11.64s/it]
Distribution of cell niches (centers) converge at iteration 2.
3 cell niches left.
Merging cell niche 20 and cell niche 71...
Done!
Assigning cells to cell niche...
Current state: [7, 20]
1%| | 1/100 [00:11<18:37, 11.29s/it]
Distribution of cell niches (centers) converge at iteration 2.
2 cell niches left.
Niche count no more than 2.
Finished!
Metric-guided solution selection
metric=’jsd_v2’: bootstrap-based minJSD strategy
The results of niche assignments are saved in .obs[niche_key]
[9]:
adata_list, adata_concat = model.select_solution(n_niche=None, # default
niche_key=f'niche_label', # default
auto=True, # default
metric='jsd_v2', # default
threshold=0.1, # default
return_adata=True, # default
plot=True, # default
save=False, # default
fig_size=(24, 6), # default
save_dir=save_dir,
file_name=f'score_vs_nichecount_basic.pdf',
)
Automatically selecting best solution...
100%|██████████| 100/100 [00:12<00:00, 8.00it/s]
100%|██████████| 100/100 [00:12<00:00, 7.71it/s]
100%|██████████| 100/100 [00:11<00:00, 8.40it/s]
100%|██████████| 100/100 [00:13<00:00, 7.68it/s]
100%|██████████| 100/100 [00:11<00:00, 8.53it/s]
100%|██████████| 100/100 [00:12<00:00, 8.26it/s]
100%|██████████| 100/100 [00:11<00:00, 8.45it/s]
100%|██████████| 100/100 [00:12<00:00, 7.91it/s]
100%|██████████| 100/100 [00:11<00:00, 8.48it/s]
100%|██████████| 100/100 [00:12<00:00, 8.21it/s]
100%|██████████| 100/100 [00:11<00:00, 8.74it/s]
100%|██████████| 100/100 [00:11<00:00, 8.73it/s]
100%|██████████| 100/100 [00:12<00:00, 8.30it/s]
100%|██████████| 100/100 [00:12<00:00, 8.27it/s]
100%|██████████| 100/100 [00:11<00:00, 8.70it/s]
100%|██████████| 100/100 [00:11<00:00, 8.48it/s]
100%|██████████| 100/100 [00:11<00:00, 8.50it/s]
100%|██████████| 100/100 [00:11<00:00, 8.51it/s]
100%|██████████| 100/100 [00:11<00:00, 8.56it/s]
100%|██████████| 100/100 [00:12<00:00, 8.31it/s]
100%|██████████| 100/100 [00:11<00:00, 8.55it/s]
100%|██████████| 100/100 [00:12<00:00, 7.97it/s]
100%|██████████| 100/100 [00:13<00:00, 7.33it/s]
100%|██████████| 100/100 [00:12<00:00, 8.24it/s]
100%|██████████| 100/100 [00:11<00:00, 8.63it/s]
100%|██████████| 100/100 [00:11<00:00, 8.35it/s]
100%|██████████| 100/100 [00:11<00:00, 8.93it/s]
100%|██████████| 100/100 [00:11<00:00, 8.54it/s]
100%|██████████| 100/100 [00:11<00:00, 9.07it/s]
100%|██████████| 100/100 [00:11<00:00, 8.42it/s]
100%|██████████| 100/100 [00:11<00:00, 8.63it/s]
100%|██████████| 100/100 [00:11<00:00, 9.03it/s]
100%|██████████| 100/100 [00:11<00:00, 8.58it/s]
100%|██████████| 100/100 [00:12<00:00, 8.08it/s]
100%|██████████| 100/100 [00:11<00:00, 8.90it/s]
100%|██████████| 100/100 [00:11<00:00, 8.33it/s]
100%|██████████| 100/100 [00:10<00:00, 9.26it/s]
100%|██████████| 100/100 [00:11<00:00, 8.72it/s]
100%|██████████| 100/100 [00:11<00:00, 9.06it/s]
100%|██████████| 100/100 [00:11<00:00, 8.37it/s]
100%|██████████| 100/100 [00:13<00:00, 7.18it/s]
100%|██████████| 100/100 [00:14<00:00, 6.98it/s]
100%|██████████| 100/100 [00:14<00:00, 6.99it/s]
100%|██████████| 100/100 [00:14<00:00, 7.06it/s]
100%|██████████| 100/100 [00:11<00:00, 8.68it/s]
100%|██████████| 100/100 [00:11<00:00, 8.68it/s]
100%|██████████| 100/100 [00:13<00:00, 7.60it/s]
100%|██████████| 100/100 [00:12<00:00, 8.27it/s]
100%|██████████| 100/100 [00:11<00:00, 8.77it/s]
100%|██████████| 100/100 [00:11<00:00, 8.77it/s]
100%|██████████| 100/100 [00:11<00:00, 8.79it/s]
100%|██████████| 100/100 [00:11<00:00, 8.88it/s]
100%|██████████| 100/100 [00:10<00:00, 9.63it/s]
100%|██████████| 100/100 [00:10<00:00, 9.12it/s]
100%|██████████| 100/100 [00:11<00:00, 8.98it/s]
100%|██████████| 100/100 [00:11<00:00, 8.91it/s]
100%|██████████| 100/100 [00:10<00:00, 9.37it/s]
100%|██████████| 100/100 [00:11<00:00, 8.62it/s]
100%|██████████| 100/100 [00:11<00:00, 8.90it/s]
100%|██████████| 100/100 [00:13<00:00, 7.20it/s]
100%|██████████| 100/100 [00:14<00:00, 7.12it/s]
100%|██████████| 100/100 [00:13<00:00, 7.35it/s]
100%|██████████| 100/100 [00:12<00:00, 8.05it/s]
100%|██████████| 100/100 [00:13<00:00, 7.60it/s]
100%|██████████| 100/100 [00:12<00:00, 8.31it/s]
100%|██████████| 100/100 [00:12<00:00, 8.12it/s]
100%|██████████| 100/100 [00:11<00:00, 8.34it/s]
100%|██████████| 100/100 [00:12<00:00, 8.13it/s]
100%|██████████| 100/100 [00:12<00:00, 8.15it/s]
100%|██████████| 100/100 [00:11<00:00, 8.66it/s]
100%|██████████| 100/100 [00:12<00:00, 8.07it/s]
100%|██████████| 100/100 [00:11<00:00, 8.45it/s]
100%|██████████| 100/100 [00:11<00:00, 8.71it/s]
100%|██████████| 100/100 [00:11<00:00, 9.04it/s]
100%|██████████| 100/100 [00:10<00:00, 9.37it/s]
100%|██████████| 100/100 [00:11<00:00, 8.78it/s]
100%|██████████| 100/100 [00:11<00:00, 8.95it/s]
100%|██████████| 100/100 [00:10<00:00, 9.14it/s]
100%|██████████| 100/100 [00:10<00:00, 9.11it/s]
100%|██████████| 100/100 [00:10<00:00, 9.42it/s]
100%|██████████| 100/100 [00:11<00:00, 9.03it/s]
100%|██████████| 100/100 [00:14<00:00, 7.08it/s]
100%|██████████| 100/100 [00:14<00:00, 6.70it/s]
100%|██████████| 100/100 [00:14<00:00, 6.89it/s]
100%|██████████| 100/100 [00:14<00:00, 7.04it/s]
100%|██████████| 100/100 [00:14<00:00, 6.71it/s]
100%|██████████| 100/100 [00:13<00:00, 7.50it/s]
100%|██████████| 100/100 [00:11<00:00, 8.91it/s]
100%|██████████| 100/100 [00:11<00:00, 8.33it/s]
100%|██████████| 100/100 [00:11<00:00, 8.53it/s]
100%|██████████| 100/100 [00:11<00:00, 8.89it/s]
100%|██████████| 100/100 [00:11<00:00, 8.67it/s]
100%|██████████| 100/100 [00:11<00:00, 9.02it/s]
100%|██████████| 100/100 [00:11<00:00, 8.89it/s]
100%|██████████| 100/100 [00:10<00:00, 9.53it/s]
100%|██████████| 100/100 [00:11<00:00, 8.87it/s]
100%|██████████| 100/100 [00:11<00:00, 8.46it/s]
100%|██████████| 100/100 [00:14<00:00, 7.09it/s]
100%|██████████| 100/100 [00:14<00:00, 6.88it/s]
Suggested range of niche count is [61, 62].
Recommended number of niches are [62]
Selecting 62 niches as the best solution.
Done!
Save and reload the results
[ ]:
adata_concat.write_h5ad(save_dir + f'Harmonics_result_0.h5ad')
[10]:
slice_name_list = np.loadtxt(data_dir + f"slice_name_list.txt", dtype=str, delimiter=" ").tolist()
adata_new = ad.read_h5ad(save_dir + f'Harmonics_result_0.h5ad')
# adata_new = adata_concat.copy()
adata_new
[10]:
AnnData object with n_obs × n_vars = 2081549 × 1122
obs: 'donor_id', 'cell_type_ontology_term_id', 'tissue_ontology_term_id', 'subclass_transfer', 'major_brain_region', 'ccf_region_name', 'brain_section_label', 'cell_type', 'tissue', 'development_stage', 'slice_name', 'celltype_idx', 'n_neighbors', 'niche_label_jsd', 'niche_label_jsd_v2', 'niche_label_fmi', 'niche_label_ari', 'niche_label_nmi', 'niche_label_asw', 'niche_label_js_asw', 'niche_label_fisher', 'niche_label_chi', 'niche_label_dbi', 'niche_label_dass_min', 'niche_label_dass_mean', 'niche_label_dafisher', 'niche_label_dachi', 'niche_label_0.09', 'niche_label_0.11', 'niche_label_100', 'niche_label_99', 'niche_label_98', 'niche_label_97', 'niche_label_96', 'niche_label_95', 'niche_label_94', 'niche_label_93', 'niche_label_92', 'niche_label_91', 'niche_label_90', 'niche_label_89', 'niche_label_88', 'niche_label_87', 'niche_label_86', 'niche_label_85', 'niche_label_84', 'niche_label_83', 'niche_label_82', 'niche_label_81', 'niche_label_80', 'niche_label_79', 'niche_label_78', 'niche_label_77', 'niche_label_76', 'niche_label_75', 'niche_label_74', 'niche_label_73', 'niche_label_72', 'niche_label_71', 'niche_label_70', 'niche_label_69', 'niche_label_68', 'niche_label_67', 'niche_label_66', 'niche_label_65', 'niche_label_64', 'niche_label_63', 'niche_label_62', 'niche_label_61', 'niche_label_60', 'niche_label_59', 'niche_label_58', 'niche_label_57', 'niche_label_56', 'niche_label_55', 'niche_label_54', 'niche_label_53', 'niche_label_52', 'niche_label_51', 'niche_label_50', 'niche_label_49', 'niche_label_48', 'niche_label_47', 'niche_label_46', 'niche_label_45', 'niche_label_44', 'niche_label_43', 'niche_label_42', 'niche_label_41', 'niche_label_40', 'niche_label_39', 'niche_label_38', 'niche_label_37', 'niche_label_36', 'niche_label_35', 'niche_label_34', 'niche_label_33', 'niche_label_32', 'niche_label_31', 'niche_label_30', 'niche_label_29', 'niche_label_28', 'niche_label_27', 'niche_label_26', 'niche_label_25', 'niche_label_24', 'niche_label_23', 'niche_label_22', 'niche_label_21', 'niche_label_20', 'niche_label_19', 'niche_label_18', 'niche_label_17', 'niche_label_16', 'niche_label_15', 'niche_label_14', 'niche_label_13', 'niche_label_12', 'niche_label_11', 'niche_label_10', 'niche_label_9', 'niche_label_8', 'niche_label_7', 'niche_label_6', 'niche_label_5', 'niche_label_4', 'niche_label_3', 'niche_label_2', 'niche_label'
uns: 'ct2idx', 'idx2ct', 'niche_cell_count', 'niche_dist', 'niche_label_summary', 'score_dict'
obsm: 'micro_dist', 'onehot', 'spatial', 'spatial_ccf'
Plot the results
[11]:
niches = adata_new.uns['niche_label_summary']
niche_colors = ['#ffff00', '#1ce6ff', '#ff34ff', '#ff4a46', '#008941', '#006fa6', '#a30059', '#ffdbe5', '#7a4900', '#0000a6',
'#63ffac', '#b79762', '#004d43', '#8fb0ff', '#997d87', '#5a0007', '#809693', '#6a3a4c', '#1b4400', '#4fc601',
'#3b5dff', '#4a3b53', '#ff2f80', '#61615a', '#001e09', '#6b7900', '#00c2a0', '#ffaa92', '#ff90c9', '#b903aa',
'#d16100', '#ddefff', '#000035', '#7b4f4b', '#a1c299', '#300018', '#0aa6d8', '#013349', '#00846f', '#372101',
'#ffb500', '#c2ffed', '#a079bf', '#cc0744', '#c0b9b2', '#c2ff99', '#ba0900', '#00489c', "#E1D656", '#6f0062',
'#0cbd66', '#eec3ff', '#456d75', '#b77b68', '#7a87a1', '#788d66', '#885578', '#fad09f', '#ff8a9a', '#d157a0',
'#bec459', '#456648', ]
niche_color_dict = {niches[i]: niche_colors[i] for i in range(len(niches))}
mbr = sorted(set(adata_new.obs['major_brain_region']))
mbr_colors = ['#1f77b4', '#ff7f0e', '#279e68', '#d62728', '#aa40fc', '#8c564b', '#e377c2', '#b5bd61', '#17becf', '#aec7e8',
'#ffbb78', '#98df8a', '#ff9896', '#c5b0d5', '#c49c94']
mbr_color_dict = {mbr[i]: mbr_colors[i] for i in range(len(mbr))}
[12]:
for i in range(len(slice_name_list)):
print(slice_name_list[i])
adata = adata_new[adata_new.obs['slice_name'] == slice_name_list[i], :].copy()
# adata_list[i].obs['niche_label'] = adata.obs['niche_label'].copy()
fig, axes = plt.subplots(1, 2, figsize=(24, 6))
# idx = 0
# niche = str(idx)
# adata.obs['highlight'] = [label if label == niche else 'others' for label in adata.obs['niche_label']]
# highlight_dict = {niche: niche_colors[idx], 'others': 'lightgray'}
# sc.pl.embedding(adata, basis='spatial', color='highlight', palette=highlight_dict,
# ax=axes[0], s=10, show=False, frameon=False, title=f'Niche {niche}', legend_fontsize=16)
sc.pl.embedding(adata, basis='spatial', color='niche_label', palette=niche_color_dict,
ax=axes[0], s=10, show=False, frameon=False, title='Cell Niche', legend_fontsize=16)
sc.pl.embedding(adata, basis='spatial', color='major_brain_region', palette=mbr_color_dict,
ax=axes[1], s=10, show=False, frameon=False, title='Major Region Annotation', legend_fontsize=16)
x_mean = (max(adata.obsm['spatial'][:, 0]) + min(adata.obsm['spatial'][:, 0])) / 2
y_mean = (max(adata.obsm['spatial'][:, 1]) + min(adata.obsm['spatial'][:, 1])) / 2
if x_mean-7800 > min(adata.obsm['spatial'][:, 0]) or x_mean+7800 < max(adata.obsm['spatial'][:, 0]):
raise ValueError('Please adjust the x_lim !')
if y_mean-3600 > min(adata.obsm['spatial'][:, 1]) or y_mean+3600 < max(adata.obsm['spatial'][:, 1]):
raise ValueError('Please adjust the y_lim !')
for ax in axes:
ax.set_xlim(x_mean-7800, x_mean+7800)
ax.set_ylim(y_mean-3600, y_mean+3600)
niche_legend_elements = [
Line2D([0], [0], marker='o', color='w', label=label,
markerfacecolor=color, markersize=10)
for label, color in niche_color_dict.items()
]
axes[0].legend(handles=niche_legend_elements, loc=(1.05, -0.1), frameon=False, ncol=3)
axes[0].axis('off')
mbr_legend_elements = [
Line2D([0], [0], marker='o', color='w', label=label,
markerfacecolor=color, markersize=10)
for label, color in mbr_color_dict.items()
]
axes[1].legend(handles=mbr_legend_elements, loc=(1.05, 0.1), frameon=False, ncol=1)
axes[1].axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.001
C57BL6J-3.002
C57BL6J-3.003
C57BL6J-3.004
C57BL6J-3.005
C57BL6J-3.006
C57BL6J-3.007
C57BL6J-3.008
C57BL6J-3.009
C57BL6J-3.010
C57BL6J-3.011
C57BL6J-3.012
C57BL6J-3.013
C57BL6J-3.015
C57BL6J-3.016
C57BL6J-3.017
C57BL6J-3.019
C57BL6J-3.020
C57BL6J-3.021
C57BL6J-3.022
C57BL6J-3.023
C57BL6J-3.024
C57BL6J-3.025
Match niches to major brain regions
[13]:
import scipy.cluster.hierarchy as sch
import scipy.spatial.distance as ssd
# indices = np.where(adata_new.obs['major_brain_region'] != 'n/a')[0]
region = adata_new.obs['major_brain_region']
niche_labels = adata_new.obs['niche_label']
cross_table = pd.crosstab(region, niche_labels)
normalized_matrix = cross_table.div(cross_table.sum(axis=0), axis=1)
col_dist = ssd.pdist(normalized_matrix.values.T, metric='correlation')
col_link = sch.linkage(col_dist, method='average')
col_link_opt = sch.optimal_leaf_ordering(col_link, ssd.squareform(col_dist))
col_indices = sch.dendrogram(col_link_opt, no_plot=True)['leaves'][::-1]
row_labels = ['Midbrain', 'Pons', 'Medulla', 'n/a', 'Cerebellum', 'Ventricular_systems', 'Fiber_tracts',
'Thalamus', 'Hypothalamus', 'Pallidum', 'Striatum', 'Olfactory', 'Cortical_subplate', 'Isocortex', 'Hippocampus', ]
sorted_matrix = normalized_matrix.loc[row_labels].iloc[:, col_indices]
fig, ax = plt.subplots(figsize=(20, 8))
sns.heatmap(
sorted_matrix,
ax=ax,
cmap='viridis',
cbar_kws={'label': 'Percentage'},
linewidths=0.,
linecolor='gray',
)
ax.grid(False)
ax.set_title('Main Region-Cell Niche Distribution', fontsize=18, pad=16)
ax.set_xlabel('Cell Niches', fontsize=16, labelpad=12)
ax.set_ylabel('Major Brain Regions', fontsize=16, labelpad=12)
ax.set_xticks(np.arange(len(sorted_matrix.columns)) + 0.5)
ax.set_xticklabels(sorted_matrix.columns, rotation=90, ha='center', fontsize=16)
ax.set_yticks(np.arange(len(sorted_matrix.index)) + 0.5)
ax.set_yticklabels(sorted_matrix.index, rotation=0, ha='right', fontsize=16)
plt.tight_layout()
plt.show()
Cell type enrichment analysis
[14]:
ct_df = ct_enrichment_test(adata_new.uns['niche_dist'],
adata_new.uns['niche_cell_count'],
adata_new.uns['idx2ct'],
adata_new.uns['niche_label_summary'],
method='fisher',
alpha=0.05,
fdr_method='fdr_by',
log2fc_threshold=1,
prop_threshold=0.01,
verbose=True,
)
ct_df.head()
62 niches and 338 cell types in total.
[14]:
| niche_idx | niche | celltype_idx | celltype | oddsratio | p-value | q-value | log2fc | prop | enrichment | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 0 | 0 | ABC NN | 0.072610 | 1.054303e-91 | 1.509361e-89 | -3.778787 | 0.000267 | False |
| 1 | 0 | 0 | 1 | ACB-BST-FS D1 Gaba | 0.526762 | 6.277995e-03 | 1.277208e-01 | -0.924496 | 0.000216 | False |
| 2 | 0 | 0 | 2 | AD Serpinb7 Glut | 0.000000 | 1.498499e-05 | 4.148422e-04 | -20.552456 | 0.000000 | False |
| 3 | 0 | 0 | 3 | ADP-MPO Trp73 Glut | 0.000000 | 9.613958e-08 | 3.101726e-06 | -21.047221 | 0.000000 | False |
| 4 | 0 | 0 | 4 | AHN Onecut3 Gaba | 0.000000 | 2.112220e-14 | 9.267749e-13 | -21.986109 | 0.000000 | False |
[15]:
niche_labels = adata_new.uns['niche_label_summary'].copy()
ct_labels = sorted(adata_new.obs['subclass_transfer'].unique())
cn_dist_count = adata_new.uns['niche_dist'].toarray() * adata_new.uns['niche_cell_count'][:, np.newaxis]
cn_dist_norm = cn_dist_count / np.sum(cn_dist_count, axis=0)
matrix_df_norm = pd.DataFrame(
data=cn_dist_norm.T,
index=ct_labels,
columns=niche_labels,
)
ct_df['stars'] = ct_df['q-value'].apply(p2stars)
stars_df = pd.DataFrame(
'',
index=matrix_df_norm.index,
columns=matrix_df_norm.columns
)
for _, row in ct_df[ct_df['enrichment']].iterrows():
niche = row['niche']
ct = row['celltype']
if (ct in stars_df.index) and (niche in stars_df.columns):
stars_df.loc[ct, niche] = row['stars']
fig, ax = plt.subplots(1, 1, figsize=(25, 80))
sns_heatmap = sns.heatmap(
matrix_df_norm,
cmap='Blues',
# cbar_kws={'label': 'Cell type proportion'},
linewidths=0.5,
linecolor='gray',
# square=True,
ax=ax,
cbar_kws={
'label': 'Cell type proportion',
'shrink': 0.4,
'aspect': 50
}
)
for i, ct in enumerate(matrix_df_norm.index):
for j, niche in enumerate(matrix_df_norm.columns):
star = stars_df.iloc[i, j]
if star:
if matrix_df_norm.iloc[i, j] > np.max(matrix_df_norm.values) * 0.7:
color='white'
else:
color='black'
ax.text(j + 0.5, i + 0.6, star, ha='center', va='center', color=color, fontsize=16, fontweight='bold')
n_rows, n_cols = matrix_df_norm.shape
ax.plot([0, n_cols], [n_rows, n_rows], color='gray', linewidth=0.5, clip_on=False)
ax.plot([n_cols, n_cols], [0, n_rows], color='gray', linewidth=0.5, clip_on=False)
ax.set_xticks(np.arange(len(matrix_df_norm.columns)) + 0.5)
ax.set_xticklabels(matrix_df_norm.columns, rotation=90, ha='center', fontsize=16)
ax.set_yticks(np.arange(len(matrix_df_norm.index)) + 0.5)
ax.set_yticklabels(matrix_df_norm.index, rotation=0, ha='right', fontsize=16)
ax.set_ylabel('Cell Type', fontsize=16)
ax.set_xlabel('Niche', fontsize=16)
ax.set_title('Row Normalized Cell Type Proportions', fontsize=16)
ax.collections[0].colorbar.ax.yaxis.label.set_size(16)
ax.collections[0].colorbar.ax.tick_params(labelsize=14)
ax.grid(False)
plt.tight_layout()
plt.show()
Annotate each niche with top enriched cell types
[16]:
slice_niche_counts = pd.crosstab(
adata_new.obs['slice_name'],
adata_new.obs['niche_label']
)
slice_niche_counts
[16]:
| niche_label | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| slice_name | |||||||||||||||||||||
| C57BL6J-3.001 | 2710 | 6586 | 5215 | 1949 | 2080 | 9741 | 1 | 6260 | 1465 | 73 | ... | 2963 | 3746 | 1278 | 0 | 1043 | 6221 | 2844 | 3011 | 0 | 434 |
| C57BL6J-3.002 | 4301 | 1874 | 8815 | 1619 | 1811 | 14138 | 44 | 12427 | 1346 | 187 | ... | 1961 | 3533 | 1429 | 0 | 1691 | 9680 | 856 | 4053 | 0 | 493 |
| C57BL6J-3.003 | 6066 | 1862 | 8549 | 1890 | 1831 | 18150 | 259 | 5607 | 1181 | 252 | ... | 1624 | 4632 | 1761 | 522 | 2089 | 10678 | 546 | 3364 | 0 | 773 |
| C57BL6J-3.004 | 6612 | 1817 | 2435 | 2395 | 1615 | 17075 | 553 | 4027 | 1185 | 283 | ... | 116 | 4722 | 2022 | 5002 | 1915 | 12159 | 331 | 3914 | 32 | 359 |
| C57BL6J-3.005 | 5302 | 2891 | 1587 | 2177 | 1312 | 15446 | 397 | 3928 | 962 | 282 | ... | 13 | 4940 | 2167 | 3877 | 1440 | 11180 | 744 | 4108 | 253 | 146 |
| C57BL6J-3.006 | 4681 | 3033 | 1610 | 1946 | 910 | 13688 | 43 | 3065 | 356 | 272 | ... | 67 | 3943 | 2194 | 3038 | 2543 | 10162 | 565 | 3834 | 448 | 0 |
| C57BL6J-3.007 | 3764 | 4126 | 3332 | 2252 | 574 | 10650 | 38 | 3981 | 16 | 394 | ... | 27 | 4484 | 1379 | 964 | 1946 | 8409 | 414 | 4188 | 759 | 0 |
| C57BL6J-3.008 | 3475 | 4139 | 4706 | 2643 | 8 | 7128 | 278 | 3545 | 0 | 473 | ... | 3 | 4534 | 1221 | 630 | 1621 | 7459 | 816 | 5137 | 2869 | 0 |
| C57BL6J-3.009 | 3338 | 4532 | 6276 | 2576 | 2 | 7032 | 693 | 3507 | 0 | 956 | ... | 8 | 4521 | 984 | 428 | 1122 | 4975 | 1736 | 5324 | 2215 | 0 |
| C57BL6J-3.010 | 3492 | 4405 | 6214 | 3534 | 0 | 6499 | 685 | 3599 | 0 | 1270 | ... | 0 | 4591 | 879 | 386 | 1444 | 4055 | 3388 | 4485 | 2347 | 273 |
| C57BL6J-3.011 | 4299 | 4990 | 0 | 3920 | 0 | 15939 | 856 | 3569 | 0 | 1512 | ... | 0 | 3897 | 0 | 198 | 2656 | 2332 | 0 | 2323 | 1791 | 415 |
| C57BL6J-3.012 | 3869 | 4784 | 0 | 4384 | 0 | 13451 | 553 | 3264 | 0 | 1669 | ... | 0 | 1985 | 0 | 164 | 1585 | 490 | 0 | 1457 | 10 | 441 |
| C57BL6J-3.013 | 3847 | 5351 | 0 | 5312 | 0 | 13686 | 180 | 3617 | 0 | 1994 | ... | 0 | 34 | 0 | 184 | 54 | 92 | 0 | 1944 | 0 | 210 |
| C57BL6J-3.015 | 5375 | 5096 | 0 | 4826 | 0 | 18618 | 0 | 2920 | 0 | 2403 | ... | 0 | 0 | 0 | 175 | 0 | 7 | 0 | 1372 | 0 | 0 |
| C57BL6J-3.016 | 3674 | 5570 | 0 | 5118 | 0 | 4154 | 0 | 3062 | 0 | 1044 | ... | 0 | 0 | 0 | 92 | 0 | 1 | 0 | 286 | 0 | 0 |
| C57BL6J-3.017 | 3672 | 5496 | 0 | 4749 | 0 | 6199 | 0 | 3132 | 0 | 710 | ... | 0 | 0 | 0 | 26 | 0 | 0 | 0 | 0 | 0 | 0 |
| C57BL6J-3.019 | 3528 | 6210 | 0 | 3004 | 0 | 6141 | 0 | 3280 | 0 | 192 | ... | 0 | 0 | 0 | 181 | 0 | 0 | 0 | 0 | 0 | 0 |
| C57BL6J-3.020 | 2626 | 6133 | 0 | 2390 | 0 | 4361 | 0 | 3343 | 0 | 42 | ... | 0 | 0 | 0 | 131 | 0 | 0 | 0 | 0 | 0 | 0 |
| C57BL6J-3.021 | 2106 | 5818 | 0 | 2588 | 0 | 4690 | 0 | 3457 | 0 | 14 | ... | 0 | 0 | 0 | 51 | 0 | 0 | 0 | 0 | 0 | 0 |
| C57BL6J-3.022 | 1048 | 5960 | 0 | 2811 | 0 | 2061 | 0 | 4541 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| C57BL6J-3.023 | 517 | 5641 | 0 | 2453 | 0 | 96 | 0 | 4969 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| C57BL6J-3.024 | 348 | 7352 | 0 | 2106 | 0 | 4 | 0 | 622 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| C57BL6J-3.025 | 64 | 3076 | 0 | 2324 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
23 rows × 62 columns
[17]:
n_niche = len(sorted(set(ct_df['niche'])))
n_top = 5
top_enriched = pd.DataFrame(
index=[str(i) for i in range(n_niche)],
columns=[f'Top{i+1}' for i in range(n_top)],
dtype=object
)
for i in range(n_niche):
niche = str(i)
sub = ct_df[(ct_df['niche'] == niche) & (ct_df['enrichment'] == True)]
# celltypes_in_sub = sub['celltype'].unique()
# adata = adata_new[(adata_new.obs['niche_label'] == niche) & (adata_new.obs['subclass_transfer'].isin(celltypes_in_sub))]
# top = adata.obs['subclass_transfer'].value_counts().head(n_top).index.tolist()
sub_sorted = sub.sort_values('log2fc', ascending=False)
top = sub_sorted['celltype'].tolist()[:n_top]
top += ['/'] * (n_top - len(top))
top_enriched.loc[niche] = top
# plot
counts = slice_niche_counts[niche]
if counts[4] > 100:
idx=4
else:
idx = counts.values.argmax()
print(slice_name_list[idx])
print(top)
adata = adata_new[adata_new.obs['slice_name'] == slice_name_list[idx], :].copy()
adata.obs['highlight'] = [label if label == niche else 'others' for label in adata.obs['niche_label']]
highlight_dict = {niche: niche_colors[i], 'others': 'lightgray'}
fig, axes = plt.subplots(1, 2, figsize=(22, 6))
sc.pl.embedding(adata, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[0], s=10, show=False, frameon=False, title=f'Niche {niche}', legend_fontsize=16)
sc.pl.embedding(adata, basis='spatial', color='major_brain_region', palette=mbr_color_dict,
ax=axes[1], s=10, show=False, frameon=False, title='Major Region Annotation', legend_fontsize=16)
x_mean = (max(adata.obsm['spatial'][:, 0]) + min(adata.obsm['spatial'][:, 0])) / 2
y_mean = (max(adata.obsm['spatial'][:, 1]) + min(adata.obsm['spatial'][:, 1])) / 2
if x_mean-7800 > min(adata.obsm['spatial'][:, 0]) or x_mean+7800 < max(adata.obsm['spatial'][:, 0]):
raise ValueError('Please adjust the x_lim !')
if y_mean-3600 > min(adata.obsm['spatial'][:, 1]) or y_mean+3600 < max(adata.obsm['spatial'][:, 1]):
raise ValueError('Please adjust the y_lim !')
for ax in axes:
ax.set_xlim(x_mean-7800, x_mean+7800)
ax.set_ylim(y_mean-3600, y_mean+3600)
niche_legend_elements = [
Line2D([0], [0], marker='o', color='w', label=label,
markerfacecolor=color, markersize=10)
for label, color in highlight_dict.items()
]
axes[0].legend(handles=niche_legend_elements, loc=(1.05, 0.5), frameon=False, ncol=1)
axes[0].axis('off')
mbr_legend_elements = [
Line2D([0], [0], marker='o', color='w', label=label,
markerfacecolor=color, markersize=10)
for label, color in mbr_color_dict.items()
]
axes[1].legend(handles=mbr_legend_elements, loc=(1.05, 0.1), frameon=False, ncol=1)
axes[1].axis('off')
plt.tight_layout()
plt.show()
# top_enriched.to_csv(save_dir + 'top_enriched_ct.csv', index=True, header=True, encoding='utf-8')
top_enriched.head()
C57BL6J-3.005
['CBX Purkinje Gaba', 'Bergmann NN', 'CB PLI Gly-Gaba', 'CBX MLI Megf11 Gaba', 'CBX MLI Cdh22 Gaba']
C57BL6J-3.005
['L4/5 IT CTX Glut', 'Pvalb Gaba', 'Vip Gaba', 'Sst Gaba', 'L5 IT CTX Glut']
C57BL6J-3.005
['OB-out Frmd7 Gaba', 'OB Dopa-Gaba', 'OB Meis2 Thsd7b Gaba', 'OB Eomes Ms4a15 Glut', 'Astro-OLF NN']
C57BL6J-3.005
['HPF CR Glut', 'ABC NN', 'Lamp5 Gaba', 'VLMC NN', 'Astro-TE NN']
C57BL6J-3.005
['SCsg Pde5a Glut', 'SC Otx2 Gcnt4 Gaba', 'SC Lef1 Otx2 Gaba', 'SCs Pax7 Nfia Gaba', 'SCsg Gabrr2 Gaba']
C57BL6J-3.005
['CB Granule Glut', 'Astro-CB NN', '/', '/', '/']
C57BL6J-3.005
['MEA-BST Lhx6 Sp9 Gaba', 'MEA-BST Sox6 Gaba', 'MEA-BST Lhx6 Nr2e1 Gaba', 'MEA-BST Otp Zic2 Glut', 'BST Tac2 Gaba']
C57BL6J-3.005
['L5 ET CTX Glut', 'L5 IT CTX Glut', 'L5 NP CTX Glut', 'L5/6 IT TPE-ENT Glut', 'Pvalb Gaba']
C57BL6J-3.005
['DMH Hmx2 Glut', 'DMH-LHA Gsx1 Gaba', 'DMH Nkx2-4 Glut', 'PH-ant-LHA Otp Bsx Glut', 'DMH-LHA Vgll2 Glut']
C57BL6J-3.005
['CA3 Glut', 'CA2-FC-IG Glut', 'RHP-COA Ndnf Gaba', 'Lamp5 Lhx6 Gaba', 'Sst Gaba']
C57BL6J-3.005
['CBX MLI Cdh22 Gaba', 'CBX MLI Megf11 Gaba', 'ABC NN', 'VLMC NN', 'SMC NN']
C57BL6J-3.005
['PB Lmx1a Glut', 'PB Evx2 Glut', 'PB Pax5 Glut', 'PB-NTS Phox2b Ebf3 Lmx1b Glut', 'CUN-PPN Evx2 Meis2 Glut']
C57BL6J-3.005
['Ependymal NN', 'TU-ARH Otp Six6 Gaba', 'MEA-BST Sox6 Gaba', 'RT-ZI Gnb3 Gaba', 'OB-STR-CTX Inh IMN']
C57BL6J-3.005
['L6 CT CTX Glut', 'L6 IT CTX Glut', 'L6b CTX Glut', 'CLA-EPd-CTX Car3 Glut', 'L5 NP CTX Glut']
C57BL6J-3.005
['OB Trdn Gaba', 'OB-in Frmd7 Gaba', 'Astro-OLF NN', 'OB-STR-CTX Inh IMN', 'OB Eomes Ms4a15 Glut']
C57BL6J-3.005
['SCig Tfap2b Chrnb3 Glut', 'SCiw Pitx2 Glut', 'SCs Lef1 Gli3 Gaba', 'SCig Foxb1 Glut', 'SC Tnnt1 Gli3 Gaba']
C57BL6J-3.020
['L6b EPd Glut', 'IT EP-CLA Glut', 'CLA-EPd-CTX Car3 Glut', 'L6b CTX Glut', 'L6b/CT ENT Glut']
C57BL6J-3.005
['CEA-BST Six3 Cyp26b1 Gaba', 'CEA-BST Ebf1 Pdyn Gaba', 'CEA-BST Rai14 Pdyn Crh Gaba', 'ACB-BST-FS D1 Gaba', 'CEA-AAA-BST Six3 Sp9 Gaba']
C57BL6J-3.005
['DG Glut', 'DG-PIR Ex IMN', '/', '/', '/']
C57BL6J-3.005
['L2/3 IT RSP Glut', 'L4 RSP-ACA Glut', 'Vip Gaba', 'L5 IT CTX Glut', 'Lamp5 Gaba']
C57BL6J-3.005
['SMC NN', 'VLMC NN', 'BAM NN', 'Astro-NT NN', 'Endo NN']
C57BL6J-3.005
['CU-ECU-SPVI Foxb1 Glut', 'SPVI-SPVC Sall3 Nfib Gly-Gaba', 'SPVI-SPVC Sall3 Lhx1 Gly-Gaba', 'SPVI-SPVC Tlx3 Ebf3 Glut', 'PARN-MDRNd-NTS Gbx2 Gly-Gaba']
C57BL6J-3.002
['LSX Prdm12 Zeb2 Gaba', 'LSX Prdm12 Slit2 Gaba', 'LSX Nkx2-1 Gaba', 'LSX Sall3 Lmo1 Gaba', 'LSX Otx2 Gaba']
C57BL6J-3.005
['CHOR NN', 'Ependymal NN', 'MEA-BST Sox6 Gaba', 'BAM NN', 'VLMC NN']
C57BL6J-3.005
['RT-ZI Gnb3 Gaba', 'ZI Pax6 Gaba', 'Astro-NT NN', 'Oligo NN', '/']
C57BL6J-3.013
['MEA Otp Foxp2 Glut', 'MEA-COA-BMA Ccdc42 Glut', 'COAa-PAA-MEA Barhl2 Glut', 'NLOT Rho Glut', 'MEA-BST Sox6 Gaba']
C57BL6J-3.005
['APN C1ql2 Glut', 'APN C1ql4 Glut', 'PRT Tcf7l2 Gaba', 'MRN Pou3f1 C1ql4 Glut', 'PRC-PAG Pax6 Glut']
C57BL6J-3.005
['MPN-MPO-PVpo Hmx2 Glut', 'PVpo-VMPO-MPN Hmx2 Gaba', 'AVPV-MEPO-SFO Tbr1 Glut', 'MPO-ADP Lhx8 Gaba', 'SI-MPO-LPO Lhx8 Gaba']
C57BL6J-3.005
['NP SUB Glut', 'SUB-ProS Glut', 'ENTmv-PA-COAp Glut', 'CT SUB Glut', 'Sst Gaba']
C57BL6J-3.005
['NDB-SI-MA-STRv Lhx8 Gaba', 'SI-MA-LPO-LHA Skor1 Glut', 'SI-MPO-LPO Lhx8 Gaba', 'PAL-STR Gaba-Chol', 'MPO-ADP Lhx8 Gaba']
C57BL6J-3.005
['Oligo NN', 'L6b CTX Glut', '/', '/', '/']
C57BL6J-3.005
['OT D3 Folh1 Gaba', 'STR-PAL Chst9 Gaba', 'STR D1 Sema5a Gaba', 'STR D1 Gaba', '/']
C57BL6J-3.005
['CA1-ProS Glut', 'CA2-FC-IG Glut', 'SUB-ProS Glut', 'CA3 Glut', '/']
C57BL6J-3.017
['L6b/CT ENT Glut', 'NP PPP Glut', 'L5/6 IT TPE-ENT Glut', 'ENTmv-PA-COAp Glut', 'L6b EPd Glut']
C57BL6J-3.005
['STN-PSTN Pitx2 Glut', 'PMd-LHA Foxb1 Glut', 'PH-SUM Foxa1 Glut', 'IF-RL-CLI-PAG Foxa1 Glut', 'PH-LHA Foxb1 Glut']
C57BL6J-3.005
['IC Tfap2d Maf Glut', 'IC Six3 En2 Gaba', 'SCsg Gabrr2 Gaba', 'Astro-NT NN', '/']
C57BL6J-3.005
['IT AON-TT-DP Glut', 'L2/3 IT PIR-ENTl Glut', 'Pvalb Gaba', '/', '/']
C57BL6J-3.017
['LA-BLA-BMA-PA Glut', 'COAp Grxcr2 Glut', 'MEA Slc17a7 Glut', 'Vip Gaba', 'Sst Gaba']
C57BL6J-3.019
['L2/3 IT PIR-ENTl Glut', 'L2/3 IT ENT Glut', 'Vip Gaba', '/', '/']
C57BL6J-3.015
['L2/3 IT PPP Glut', 'L2 IT ENT-po Glut', 'L2/3 IT ENT Glut', 'L2 IT PPP-APr Glut', 'L2/3 IT RSP Glut']
C57BL6J-3.005
['VMH Fezf1 Glut', 'VMH Nr5a1 Glut', 'TU-ARH Otp Six6 Gaba', 'DMH Hmx2 Gaba', 'AHN Onecut3 Gaba']
C57BL6J-3.005
['STR D2 Gaba', 'STR D1 Gaba', 'STR D1 Sema5a Gaba', '/', '/']
C57BL6J-3.005
['L2/3 IT CTX Glut', 'Vip Gaba', 'Lamp5 Gaba', 'Pvalb Gaba', 'Sst Gaba']
C57BL6J-3.005
['STR-PAL Chst9 Gaba', 'STR D1 Sema5a Gaba', 'NDB-SI-ant Prdm12 Gaba', 'ACB-BST-FS D1 Gaba', 'CEA-AAA-BST Six3 Sp9 Gaba']
C57BL6J-3.005
['ABC NN', 'VLMC NN', 'CT SUB Glut', 'BAM NN', 'MEA-COA-BMA Ccdc42 Glut']
C57BL6J-3.005
['MH Tac2 Glut', 'LH Pou4f1 Sox1 Glut', 'PF Fzd5 Glut', 'SPA-SPFm-SPFp-POL-PIL-PoT Sp9 Glut', 'Astroependymal NN']
C57BL6J-3.005
['NTS-PARN Neurod2 Gly-Gaba', 'NTS Phox2b Glut', 'DMX VII Tbx20 Chol', 'PARN-MDRNd-NTS Gbx2 Gly-Gaba', 'NTS Dbh Glut']
C57BL6J-3.005
['ZI Pax6 Gaba', 'LGv-ZI Otx2 Gaba', 'LGv-SPFp-SPFm Nkx2-2 Tcf7l2 Gaba', 'PVHd-DMH Lhx6 Gaba', 'Astro-NT NN']
C57BL6J-3.005
['SCig-an-PPT Foxb1 Glut', 'PRT Mecom Gaba', 'SCs Dmbx1 Gaba', 'PRT Tcf7l2 Gaba', 'SCs Pax7 Nfia Gaba']
C57BL6J-3.005
['AHN-RCH-LHA Otp Fezf1 Glut', 'AHN-SBPV-PVHd Pdrm12 Gaba', 'AHN Onecut3 Gaba', 'PVHd-SBPV Six3 Prox1 Gaba', 'LHA Barhl2 Glut']
C57BL6J-3.005
['OB-mi Frmd7 Gaba', 'OB Eomes Ms4a15 Glut', 'Astro-OLF NN', 'OB Meis2 Thsd7b Gaba', 'OB Dopa-Gaba']
C57BL6J-3.005
['SCdg-PAG Tfap2b Glut', 'PAG Pou4f2 Glut', 'PAG Pou4f2 Mesi2 Glut', 'PAG Pou4f3 Glut', 'PAG Pou4f1 Ebf2 Glut']
C57BL6J-3.001
['DTN-LDT-IPN Otp Pax3 Gaba', 'LDT-PCG St18 Gaba', 'LDT-DTN Gata3 Nfix Gaba', 'LDT-PCG-CS Gata3 Lhx1 Gaba', 'NI-RPO Gata3 Nr4a2 Gaba']
C57BL6J-3.005
['TH Prkcd Grin2c Glut', 'CM-IAD-CL-PCN Sema5b Glut', 'AV Col27a1 Glut', 'Astro-NT NN', '/']
C57BL6J-3.005
['CUN Evx2 Lhx2 Glut', 'PAG-MRN Pou3f1 Glut', 'PAG-PPN Pax5 Sox21 Gaba', 'PAG-RN Nkx2-2 Otx1 Gaba', 'PAG-MRN-RN Foxa2 Gaba']
C57BL6J-3.005
['OB-STR-CTX Inh IMN', 'Astro-OLF NN', 'Ependymal NN', 'Astro-TE NN', '/']
C57BL6J-3.005
['PG-TRN-LRN Fat2 Glut', 'PSV Pvalb Lhx2 Glut', 'Astro-NT NN', 'Oligo NN', '/']
C57BL6J-3.005
['MY Lhx1 Gly-Gaba', 'PRNc Otp Gly-Gaba', 'PGRN-PARN-MDRN Hoxb5 Glut', 'Astro-NT NN', 'Oligo NN']
C57BL6J-3.005
['OEC NN', 'OB-out Frmd7 Gaba', 'CT SUB Glut', 'Astro-OLF NN', 'BAM NN']
C57BL6J-3.005
['MM Foxb1 Glut', 'CBN Dmbx1 Gaba', 'IO Fgl2 Glut', 'CBN Neurod2 Pvalb Glut', 'SNc-VTA-RAmb Foxa1 Dopa']
C57BL6J-3.005
['SPVC Ccdc172 Glut', 'SPVC Nmu Glut', 'SPVC Mafa Glut', 'SPVI-SPVC Sall3 Lhx1 Gly-Gaba', 'SPVI-SPVC Tlx3 Ebf3 Glut']
C57BL6J-3.005
['MEA-BST Lhx6 Nfib Gaba', 'MEA-BST Otp Zic2 Glut', 'MEA-BST Lhx6 Sp9 Gaba', 'MPN-MPO-LPO Lhx6 Zfhx3 Gaba', 'BST-MPN Six3 Nrgn Gaba']
[17]:
| Top1 | Top2 | Top3 | Top4 | Top5 | |
|---|---|---|---|---|---|
| 0 | CBX Purkinje Gaba | Bergmann NN | CB PLI Gly-Gaba | CBX MLI Megf11 Gaba | CBX MLI Cdh22 Gaba |
| 1 | L4/5 IT CTX Glut | Pvalb Gaba | Vip Gaba | Sst Gaba | L5 IT CTX Glut |
| 2 | OB-out Frmd7 Gaba | OB Dopa-Gaba | OB Meis2 Thsd7b Gaba | OB Eomes Ms4a15 Glut | Astro-OLF NN |
| 3 | HPF CR Glut | ABC NN | Lamp5 Gaba | VLMC NN | Astro-TE NN |
| 4 | SCsg Pde5a Glut | SC Otx2 Gcnt4 Gaba | SC Lef1 Otx2 Gaba | SCs Pax7 Nfia Gaba | SCsg Gabrr2 Gaba |
Niche annotation
[18]:
region_dict = {
'0': 'Cerebellar cortex, Purkinje layer',
'1': 'Isocortex, layer 4',
'2': 'MOB, glomerular layer',
'3': 'Hippocampal region/Cerebral cortex, Layer 1',
'4': 'Superior colliculus, sensory related',
'5': 'Cerebellar cortex, granular layer',
'6': 'Medial amygdalar nucleus/Bed nuclei of the stria terminalis',
'7': 'Isocortex, layer 5',
'8': 'Dorsomedial nucleus of the hypothalamus',
'9': 'CA3',
'10': 'Cerebellar cortex, molecular layer',
'11': 'Pons, parabrachial nucleus',
'12': 'Ventricular systems, ependymal enriched',
'13': 'Isocortex, layer 6',
'14': 'MOB, granule layer',
'15': 'Superior colliculus, motor related',
'16': 'Cortical subplate 1',
'17': 'Central amygdalar nucleus/Bed nuclei of the stria terminalis',
'18': 'Dentate gyrus',
'19': 'Retrosplenial area, layer 2-4',
'20': 'SMC enriched',
'21': 'Medulla 1',
'22': 'Lateral septal complex',
'23': 'Ventricular systems, choroid plexus',
'24': 'Reticular nucleus of the thalamus',
'25': 'Medial amygdalar nucleus/Cortical amygdalar area',
'26': 'Anterior pretectal nucleus',
'27': 'Medial preoptic area/Medial preoptic nucleus',
'28': 'Subiculum',
'29': 'Pallidum, substantia innominata',
'30': 'Fiber tracts',
'31': 'Olfactory tubercle',
'32': 'CA1',
'33': 'Entorhinal area',
'34': 'Posterior hypothalamic nucleus/Dorsal premammillary nucleus',
'35': 'Inferior colliculus',
'36': 'Anterior olfactory nucleus',
'37': 'Cortical subplate 2',
'38': 'Piriform area',
'39': 'Presubiculum/Postsubiculum/Entorhinal area',
'40': 'Ventromedial hypothalamic nucleus',
'41': 'Striatum, dorsal region/ventral region',
'42': 'Isocortex, Layer 2/3',
'43': 'Striatum-Pallidum interface',
'44': 'ABC enriched',
'45': 'Epithalamus/Parafascicular nucleus',
'46': 'Medulla 2',
'47': 'Zona incerta',
'48': 'Posterior pretectal nucleus',
'49': 'Anterior hypothalamic nucleus',
'50': 'MOB, external plexiform layer',
'51': 'Periaqueductal gray',
'52': 'Pons',
'53': 'Thalamus',
'54': 'Midbrain reticular nucleus/Periaqueductal gray',
'55': 'Ventricular systems, OB-STR-CTX Inh IMN enriched',
'56': 'Pontine gray/Principal sensory nucleus of the trigeminal',
'57': 'Midbrain/Pons/Medulla',
'58': 'MOB, olfactory nerve layer',
'59': 'Mammillary body/Cerebellar nuclei/Midbrain/Medulla/Pons',
'60': 'Medulla 3',
'61': 'Bed nuclei of the stria terminalis'
}
Examine the niche delineation of some ROIs
Cerebellar cortex
Show relative niches and corresponding annotations
[19]:
template = 'C57BL6J-3.008'
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
cc_niches = ['10', '0', '5']
print(template)
fig, axes = plt.subplots(1, 4, figsize=(24, 6))
sc.pl.embedding(adata, basis='spatial', color='niche_label', palette=niche_color_dict,
ax=axes[0], s=50, show=False, frameon=False, legend_loc=None)
axes[0].set_title('Cerebellar Cortex', fontsize=20)
for i, niche in enumerate(cc_niches):
adata.obs['highlight'] = [label if label == niche else 'others' for label in adata.obs['niche_label']]
highlight_dict = {niche: niche_colors[int(niche)], 'others': 'lightgray'}
sc.pl.embedding(adata, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[i+1], s=50, show=False, frameon=False, legend_loc=None)
axes[i+1].set_title(region_dict[niche], fontsize=20)
for ax in axes:
ax.set_xlim(max(adata.obsm['spatial'][:, 0])-4000, max(adata.obsm['spatial'][:, 0])-900)
ax.set_ylim(max(adata.obsm['spatial'][:, 1])-4000, max(adata.obsm['spatial'][:, 1])-900)
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.008
Extract data of corresponding niches
[20]:
adata_cc = adata_new[adata_new.obs['niche_label'].isin(cc_niches), :].copy()
adata_cc
[20]:
AnnData object with n_obs × n_vars = 337904 × 1122
obs: 'donor_id', 'cell_type_ontology_term_id', 'tissue_ontology_term_id', 'subclass_transfer', 'major_brain_region', 'ccf_region_name', 'brain_section_label', 'cell_type', 'tissue', 'development_stage', 'slice_name', 'celltype_idx', 'n_neighbors', 'niche_label_jsd', 'niche_label_jsd_v2', 'niche_label_fmi', 'niche_label_ari', 'niche_label_nmi', 'niche_label_asw', 'niche_label_js_asw', 'niche_label_fisher', 'niche_label_chi', 'niche_label_dbi', 'niche_label_dass_min', 'niche_label_dass_mean', 'niche_label_dafisher', 'niche_label_dachi', 'niche_label_0.09', 'niche_label_0.11', 'niche_label_100', 'niche_label_99', 'niche_label_98', 'niche_label_97', 'niche_label_96', 'niche_label_95', 'niche_label_94', 'niche_label_93', 'niche_label_92', 'niche_label_91', 'niche_label_90', 'niche_label_89', 'niche_label_88', 'niche_label_87', 'niche_label_86', 'niche_label_85', 'niche_label_84', 'niche_label_83', 'niche_label_82', 'niche_label_81', 'niche_label_80', 'niche_label_79', 'niche_label_78', 'niche_label_77', 'niche_label_76', 'niche_label_75', 'niche_label_74', 'niche_label_73', 'niche_label_72', 'niche_label_71', 'niche_label_70', 'niche_label_69', 'niche_label_68', 'niche_label_67', 'niche_label_66', 'niche_label_65', 'niche_label_64', 'niche_label_63', 'niche_label_62', 'niche_label_61', 'niche_label_60', 'niche_label_59', 'niche_label_58', 'niche_label_57', 'niche_label_56', 'niche_label_55', 'niche_label_54', 'niche_label_53', 'niche_label_52', 'niche_label_51', 'niche_label_50', 'niche_label_49', 'niche_label_48', 'niche_label_47', 'niche_label_46', 'niche_label_45', 'niche_label_44', 'niche_label_43', 'niche_label_42', 'niche_label_41', 'niche_label_40', 'niche_label_39', 'niche_label_38', 'niche_label_37', 'niche_label_36', 'niche_label_35', 'niche_label_34', 'niche_label_33', 'niche_label_32', 'niche_label_31', 'niche_label_30', 'niche_label_29', 'niche_label_28', 'niche_label_27', 'niche_label_26', 'niche_label_25', 'niche_label_24', 'niche_label_23', 'niche_label_22', 'niche_label_21', 'niche_label_20', 'niche_label_19', 'niche_label_18', 'niche_label_17', 'niche_label_16', 'niche_label_15', 'niche_label_14', 'niche_label_13', 'niche_label_12', 'niche_label_11', 'niche_label_10', 'niche_label_9', 'niche_label_8', 'niche_label_7', 'niche_label_6', 'niche_label_5', 'niche_label_4', 'niche_label_3', 'niche_label_2', 'niche_label'
uns: 'ct2idx', 'idx2ct', 'niche_cell_count', 'niche_dist', 'niche_label_summary', 'score_dict'
obsm: 'micro_dist', 'onehot', 'spatial', 'spatial_ccf'
Filter out cell types that comprise less than 500 cell in these niches
[21]:
min_cells = 500
ct_counts = adata_cc.obs['subclass_transfer'].value_counts()
valid_ct = ct_counts[ct_counts >= min_cells].index
adata_cc = adata_cc[adata_cc.obs['subclass_transfer'].isin(valid_ct)].copy()
adata_cc.shape, len(list(set(adata_cc.obs['subclass_transfer'])))
[21]:
((335561, 1122), 19)
Compute the number of cells in each niche
[22]:
ct_onehot = label2onehot(adata_cc.obs['subclass_transfer'])
indices = [cc_niches.index(str(label)) for label in adata_cc.obs['niche_label']]
cn_onehot = index2onehot(indices)
cell_count_niche = cn_onehot.T.sum(axis=1)
cc_dist = (cn_onehot.T @ ct_onehot) / cell_count_niche
cell_count_niche = np.array(cell_count_niche).flatten()
cc_ct = sorted(set(adata_cc.obs['subclass_transfer']))
idx2ct = {str(i): cc_ct[i] for i in range(len(cc_ct))}
cell_count_niche
[22]:
array([ 49620., 77698., 208243.])
Run cell type enrichment analysis
[23]:
cc_df = ct_enrichment_test(cc_dist,
cell_count_niche,
idx2ct,
cc_niches,
method='fisher',
alpha=0.05,
fdr_method='fdr_by',
log2fc_threshold=1,
prop_threshold=0.01,
verbose=True,
)
cc_df.head()
3 niches and 19 cell types in total.
[23]:
| niche_idx | niche | celltype_idx | celltype | oddsratio | p-value | q-value | log2fc | prop | enrichment | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 10 | 0 | ABC NN | 118.754800 | 0.000000e+00 | 0.000000e+00 | 6.863601 | 0.019549 | True |
| 1 | 0 | 10 | 1 | Astro-CB NN | 0.124247 | 6.177931e-164 | 4.794325e-163 | -2.991976 | 0.001693 | False |
| 2 | 0 | 10 | 2 | Astro-NT NN | 0.382726 | 1.338256e-34 | 7.356331e-34 | -1.379438 | 0.002660 | False |
| 3 | 0 | 10 | 3 | BAM NN | 12.437840 | 1.870073e-199 | 1.591696e-198 | 3.618020 | 0.008545 | False |
| 4 | 0 | 10 | 4 | Bergmann NN | 0.099564 | 0.000000e+00 | 0.000000e+00 | -3.245458 | 0.006530 | False |
Plot the results
[24]:
matrix_df = pd.DataFrame(
data=cc_dist.toarray(),
index=cc_niches,
columns=cc_ct,
)
cn_dist_count = cc_dist.toarray() * cell_count_niche[:, np.newaxis]
cn_dist_norm = cn_dist_count / np.sum(cn_dist_count, axis=0)
matrix_df_norm = pd.DataFrame(
data=cn_dist_norm,
index=cc_niches,
columns=cc_ct,
)
cc_df['stars'] = cc_df['q-value'].apply(p2stars)
stars_df = pd.DataFrame(
'',
index=matrix_df.index,
columns=matrix_df.columns
)
for _, row in cc_df[cc_df['enrichment']].iterrows():
niche = row['niche']
ct = row['celltype']
if (niche in stars_df.index) and (ct in stars_df.columns):
stars_df.loc[niche, ct] = row['stars']
fig, axes = plt.subplots(1, 1, figsize=(15, 2))
sns_heatmap_0 = sns.heatmap(
matrix_df,
cmap='Blues',
# cbar_kws={'label': 'Cell type proportion'},
linewidths=0.5,
linecolor='gray',
# square=True,
ax=axes
)
for i, niche in enumerate(matrix_df.index):
for j, ct in enumerate(matrix_df.columns):
star = stars_df.iloc[i, j]
if star:
if matrix_df.iloc[i, j] > np.max(matrix_df.values) * 0.7:
color='white'
else:
color='black'
axes.text(j + 0.5, i + 0.6, star, ha='center', va='center', color=color, fontsize=20, fontweight='bold')
n_rows, n_cols = matrix_df.shape
axes.plot([0, n_cols], [n_rows, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.plot([n_cols, n_cols], [0, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.set_xticklabels(axes.get_xticklabels(), rotation=90, ha='center', fontsize=20)
axes.set_yticklabels(axes.get_yticklabels(), rotation=0, ha='right', fontsize=20)
axes.set_ylabel('Niche', fontsize=20)
axes.set_xlabel('Cell Type', fontsize=20)
axes.set_title('Cell Type Proportions', fontsize=20)
axes.collections[0].colorbar.ax.yaxis.label.set_size(20)
axes.collections[0].colorbar.ax.tick_params(labelsize=16)
axes.grid(False)
plt.tight_layout()
plt.show()
fig, axes = plt.subplots(1, 1, figsize=(15, 2))
sns_heatmap_1 = sns.heatmap(
matrix_df_norm,
cmap='Blues',
# cbar_kws={'label': 'Cell type proportion'},
linewidths=0.5,
linecolor='gray',
# square=True,
ax=axes
)
for i, niche in enumerate(matrix_df.index):
for j, ct in enumerate(matrix_df.columns):
star = stars_df.iloc[i, j]
if star:
if matrix_df_norm.iloc[i, j] > np.max(matrix_df_norm.values) * 0.7:
color='white'
else:
color='black'
axes.text(j + 0.5, i + 0.6, star, ha='center', va='center', color=color, fontsize=20, fontweight='bold')
n_rows, n_cols = matrix_df.shape
axes.plot([0, n_cols], [n_rows, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.plot([n_cols, n_cols], [0, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.set_xticklabels(axes.get_xticklabels(), rotation=90, ha='center', fontsize=20)
axes.set_yticklabels(axes.get_yticklabels(), rotation=0, ha='right', fontsize=20)
axes.set_ylabel('Niche', fontsize=20)
axes.set_xlabel('Cell Type', fontsize=20)
axes.set_title('Column Normalized Cell Type Proportions', fontsize=20)
axes.collections[0].colorbar.ax.yaxis.label.set_size(20)
axes.collections[0].colorbar.ax.tick_params(labelsize=16)
axes.grid(False)
plt.tight_layout()
plt.show()
Display enriched cell type for each niche
[25]:
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
enriched_ct = ['CBX MLI Cdh22 Gaba', 'Bergmann NN', 'CBX Purkinje Gaba', 'CB Granule Glut']
print(template)
fig, axes = plt.subplots(1, 4, figsize=(24, 6))
for i, ct in enumerate(enriched_ct):
adata.obs['highlight'] = [label if label == ct else 'others' for label in adata.obs['subclass_transfer']]
highlight_dict = {ct: 'brown', 'others': 'lightgray'}
adata_tmp = adata[adata.obs['highlight'].isin(['others'])].copy()
sc.pl.embedding(adata_tmp, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[i], s=50, show=False, frameon=False, legend_loc=None)
adata_tmp = adata[~(adata.obs['highlight'].isin(['others']))].copy()
sc.pl.embedding(adata_tmp, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[i], s=50, show=False, frameon=False, legend_loc=None)
axes[i].set_title(ct, fontsize=20)
for ax in axes:
ax.set_xlim(max(adata.obsm['spatial'][:, 0])-4000, max(adata.obsm['spatial'][:, 0])-900)
ax.set_ylim(max(adata.obsm['spatial'][:, 1])-4000, max(adata.obsm['spatial'][:, 1])-900)
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.008
Preprocess data and find differentially expressed genes for each niche
[26]:
# adata_cc = adata_new[adata_new.obs['niche_label'].isin(cc_niches), :].copy()
sc.pp.filter_genes(adata_cc, min_cells=100)
sc.pp.normalize_total(adata_cc, target_sum=1e3)
sc.pp.log1p(adata_cc)
sc.tl.rank_genes_groups(adata_cc, groupby='niche_label', method='wilcoxon', pts=True)
adata_cc
[26]:
AnnData object with n_obs × n_vars = 335561 × 1122
obs: 'donor_id', 'cell_type_ontology_term_id', 'tissue_ontology_term_id', 'subclass_transfer', 'major_brain_region', 'ccf_region_name', 'brain_section_label', 'cell_type', 'tissue', 'development_stage', 'slice_name', 'celltype_idx', 'n_neighbors', 'niche_label_jsd', 'niche_label_jsd_v2', 'niche_label_fmi', 'niche_label_ari', 'niche_label_nmi', 'niche_label_asw', 'niche_label_js_asw', 'niche_label_fisher', 'niche_label_chi', 'niche_label_dbi', 'niche_label_dass_min', 'niche_label_dass_mean', 'niche_label_dafisher', 'niche_label_dachi', 'niche_label_0.09', 'niche_label_0.11', 'niche_label_100', 'niche_label_99', 'niche_label_98', 'niche_label_97', 'niche_label_96', 'niche_label_95', 'niche_label_94', 'niche_label_93', 'niche_label_92', 'niche_label_91', 'niche_label_90', 'niche_label_89', 'niche_label_88', 'niche_label_87', 'niche_label_86', 'niche_label_85', 'niche_label_84', 'niche_label_83', 'niche_label_82', 'niche_label_81', 'niche_label_80', 'niche_label_79', 'niche_label_78', 'niche_label_77', 'niche_label_76', 'niche_label_75', 'niche_label_74', 'niche_label_73', 'niche_label_72', 'niche_label_71', 'niche_label_70', 'niche_label_69', 'niche_label_68', 'niche_label_67', 'niche_label_66', 'niche_label_65', 'niche_label_64', 'niche_label_63', 'niche_label_62', 'niche_label_61', 'niche_label_60', 'niche_label_59', 'niche_label_58', 'niche_label_57', 'niche_label_56', 'niche_label_55', 'niche_label_54', 'niche_label_53', 'niche_label_52', 'niche_label_51', 'niche_label_50', 'niche_label_49', 'niche_label_48', 'niche_label_47', 'niche_label_46', 'niche_label_45', 'niche_label_44', 'niche_label_43', 'niche_label_42', 'niche_label_41', 'niche_label_40', 'niche_label_39', 'niche_label_38', 'niche_label_37', 'niche_label_36', 'niche_label_35', 'niche_label_34', 'niche_label_33', 'niche_label_32', 'niche_label_31', 'niche_label_30', 'niche_label_29', 'niche_label_28', 'niche_label_27', 'niche_label_26', 'niche_label_25', 'niche_label_24', 'niche_label_23', 'niche_label_22', 'niche_label_21', 'niche_label_20', 'niche_label_19', 'niche_label_18', 'niche_label_17', 'niche_label_16', 'niche_label_15', 'niche_label_14', 'niche_label_13', 'niche_label_12', 'niche_label_11', 'niche_label_10', 'niche_label_9', 'niche_label_8', 'niche_label_7', 'niche_label_6', 'niche_label_5', 'niche_label_4', 'niche_label_3', 'niche_label_2', 'niche_label'
var: 'n_cells'
uns: 'ct2idx', 'idx2ct', 'niche_cell_count', 'niche_dist', 'niche_label_summary', 'score_dict', 'log1p', 'rank_genes_groups'
obsm: 'micro_dist', 'onehot', 'spatial', 'spatial_ccf'
[27]:
n_top_genes = 5
result = adata_cc.uns['rank_genes_groups']
groups = result['names'].dtype.names
selected_genes = []
for group in groups:
genes = pd.Index(result["names"][group], name="gene")
pts_series = result["pts"][group]
pts_aligned = pts_series.reindex(genes)
df = pd.DataFrame({
'gene': result['names'][group],
'pvals_adj': result['pvals_adj'][group],
'logfoldchanges': result['logfoldchanges'][group],
'pts': pts_aligned.values,
})
df_filtered = df[(df['pvals_adj'] < 0.05) & (df['logfoldchanges'] > 1) & (df['pts'] > 0.3)]
topn = df_filtered.sort_values('logfoldchanges', ascending=False).head(n_top_genes)
selected_genes.extend(topn['gene'].values.tolist())
Plot the results
[28]:
result = sc.pl.heatmap(adata_cc, var_names=selected_genes, groupby='niche_label', cmap='viridis', use_raw=False, show=False, standard_scale='var')
fig = result['heatmap_ax'].get_figure()
fig.set_size_inches(6, 9)
ax = result['heatmap_ax']
ax.tick_params(axis='x', labelsize=16)
ax = result['groupby_ax']
ax.tick_params(axis='y', labelsize=16)
plt.tight_layout()
plt.show()
[29]:
result = sc.pl.dotplot(adata_cc, var_names=selected_genes, groupby='niche_label', use_raw=False, show=False, standard_scale="var", dot_max=0.6)
fig = result['mainplot_ax'].get_figure()
fig.set_size_inches(8, 3)
ax = result['mainplot_ax']
ax.tick_params(axis='x', labelsize=16)
ax.tick_params(axis='y', labelsize=16)
ax = result['size_legend_ax']
ax.tick_params(axis='x', labelsize=12)
ax = result['color_legend_ax']
ax.tick_params(axis='x', labelsize=12)
plt.tight_layout()
plt.show()
[30]:
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
print(template)
x = adata.obsm['spatial'][:, 0]
y = adata.obsm['spatial'][:, 1]
x_min, x_max = max(adata.obsm['spatial'][:, 0])-4000, max(adata.obsm['spatial'][:, 0])-900
y_min, y_max = max(adata.obsm['spatial'][:, 1])-4000, max(adata.obsm['spatial'][:, 1])-900
mask = (x >= x_min) & (x <= x_max) & (y >= y_min) & (y <= y_max)
adata_crop = adata[mask].copy()
sc.pp.normalize_total(adata_crop, target_sum=1e3)
sc.pp.log1p(adata_crop)
fig, axes = plt.subplots(3, 5, figsize=(24, 14))
for i, marker in enumerate(selected_genes):
row = i // 5
col = i % 5
sc.pl.embedding(adata_crop, basis='spatial', color=marker,
ax=axes[row, col], s=30, show=False, frameon=False,
cmap='viridis', title=marker)
axes[row, col].set_title(marker, fontsize=20)
for row in axes:
for ax in row:
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.008
Hippocampal formation
Show relative niches and corresponding annotations
[31]:
template = 'C57BL6J-3.008'
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
hf_niches = ['28', '32', '9', '18',]
print(template)
fig, axes = plt.subplots(2, 3, figsize=(20, 10))
sc.pl.embedding(adata, basis='spatial', color='niche_label', palette=niche_color_dict,
ax=axes[0, 0], s=70, show=False, frameon=False, legend_loc=None)
axes[0, 0].set_title('Hippocampal Formation', fontsize=20)
for i, niche in enumerate(hf_niches):
row_idx = (i+1) // 3
col_idx = (i+1) % 3
adata.obs['highlight'] = [label if label == niche else 'others' for label in adata.obs['niche_label']]
highlight_dict = {niche: niche_colors[int(niche)], 'others': 'lightgray'}
sc.pl.embedding(adata, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[row_idx, col_idx], s=70, show=False, frameon=False, legend_loc=None)
axes[row_idx, col_idx].set_title(region_dict[niche], fontsize=20)
for row in axes:
for ax in row:
ax.set_xlim(max(adata.obsm['spatial'][:, 0])-8150, max(adata.obsm['spatial'][:, 0])-5700)
ax.set_ylim(max(adata.obsm['spatial'][:, 1])-2300, max(adata.obsm['spatial'][:, 1])-900)
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.008
Extract data of corresponding niches
[32]:
adata_hf = adata_new[adata_new.obs['niche_label'].isin(hf_niches), :].copy()
adata_hf
[32]:
AnnData object with n_obs × n_vars = 72676 × 1122
obs: 'donor_id', 'cell_type_ontology_term_id', 'tissue_ontology_term_id', 'subclass_transfer', 'major_brain_region', 'ccf_region_name', 'brain_section_label', 'cell_type', 'tissue', 'development_stage', 'slice_name', 'celltype_idx', 'n_neighbors', 'niche_label_jsd', 'niche_label_jsd_v2', 'niche_label_fmi', 'niche_label_ari', 'niche_label_nmi', 'niche_label_asw', 'niche_label_js_asw', 'niche_label_fisher', 'niche_label_chi', 'niche_label_dbi', 'niche_label_dass_min', 'niche_label_dass_mean', 'niche_label_dafisher', 'niche_label_dachi', 'niche_label_0.09', 'niche_label_0.11', 'niche_label_100', 'niche_label_99', 'niche_label_98', 'niche_label_97', 'niche_label_96', 'niche_label_95', 'niche_label_94', 'niche_label_93', 'niche_label_92', 'niche_label_91', 'niche_label_90', 'niche_label_89', 'niche_label_88', 'niche_label_87', 'niche_label_86', 'niche_label_85', 'niche_label_84', 'niche_label_83', 'niche_label_82', 'niche_label_81', 'niche_label_80', 'niche_label_79', 'niche_label_78', 'niche_label_77', 'niche_label_76', 'niche_label_75', 'niche_label_74', 'niche_label_73', 'niche_label_72', 'niche_label_71', 'niche_label_70', 'niche_label_69', 'niche_label_68', 'niche_label_67', 'niche_label_66', 'niche_label_65', 'niche_label_64', 'niche_label_63', 'niche_label_62', 'niche_label_61', 'niche_label_60', 'niche_label_59', 'niche_label_58', 'niche_label_57', 'niche_label_56', 'niche_label_55', 'niche_label_54', 'niche_label_53', 'niche_label_52', 'niche_label_51', 'niche_label_50', 'niche_label_49', 'niche_label_48', 'niche_label_47', 'niche_label_46', 'niche_label_45', 'niche_label_44', 'niche_label_43', 'niche_label_42', 'niche_label_41', 'niche_label_40', 'niche_label_39', 'niche_label_38', 'niche_label_37', 'niche_label_36', 'niche_label_35', 'niche_label_34', 'niche_label_33', 'niche_label_32', 'niche_label_31', 'niche_label_30', 'niche_label_29', 'niche_label_28', 'niche_label_27', 'niche_label_26', 'niche_label_25', 'niche_label_24', 'niche_label_23', 'niche_label_22', 'niche_label_21', 'niche_label_20', 'niche_label_19', 'niche_label_18', 'niche_label_17', 'niche_label_16', 'niche_label_15', 'niche_label_14', 'niche_label_13', 'niche_label_12', 'niche_label_11', 'niche_label_10', 'niche_label_9', 'niche_label_8', 'niche_label_7', 'niche_label_6', 'niche_label_5', 'niche_label_4', 'niche_label_3', 'niche_label_2', 'niche_label'
uns: 'ct2idx', 'idx2ct', 'niche_cell_count', 'niche_dist', 'niche_label_summary', 'score_dict'
obsm: 'micro_dist', 'onehot', 'spatial', 'spatial_ccf'
Filter out cell types that comprise less than 500 cell in these niches
[33]:
min_cells = 500
ct_counts = adata_hf.obs['subclass_transfer'].value_counts()
valid_ct = ct_counts[ct_counts >= min_cells].index
adata_hf = adata_hf[adata_hf.obs['subclass_transfer'].isin(valid_ct)].copy()
adata_hf.shape, len(list(set(adata_hf.obs['subclass_transfer'])))
[33]:
((69208, 1122), 17)
Compute the number of cells in each niche
[34]:
ct_onehot = label2onehot(adata_hf.obs['subclass_transfer'])
indices = [hf_niches.index(str(label)) for label in adata_hf.obs['niche_label']]
cn_onehot = index2onehot(indices)
cell_count_niche = cn_onehot.T.sum(axis=1)
hf_dist = (cn_onehot.T @ ct_onehot) / cell_count_niche
cell_count_niche = np.array(cell_count_niche).flatten()
hf_ct = sorted(set(adata_hf.obs['subclass_transfer']))
idx2ct = {str(i): hf_ct[i] for i in range(len(hf_ct))}
cell_count_niche
[34]:
array([11520., 17966., 13007., 26715.])
Run cell type enrichment analysis
[35]:
hf_df = ct_enrichment_test(hf_dist,
cell_count_niche,
idx2ct,
hf_niches,
method='fisher',
alpha=0.05,
fdr_method='fdr_by',
log2fc_threshold=1,
prop_threshold=0.01,
verbose=True,
)
hf_df.head()
4 niches and 17 cell types in total.
[35]:
| niche_idx | niche | celltype_idx | celltype | oddsratio | p-value | q-value | log2fc | prop | enrichment | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 28 | 0 | Astro-TE NN | 1.448319 | 5.534405e-32 | 3.476840e-31 | 0.472642 | 0.135330 | False |
| 1 | 0 | 28 | 1 | CA1-ProS Glut | 0.105427 | 0.000000e+00 | 0.000000e+00 | -2.933487 | 0.028472 | False |
| 2 | 0 | 28 | 2 | CA2-FC-IG Glut | 0.034992 | 6.445520e-39 | 4.386659e-38 | -4.823074 | 0.000347 | False |
| 3 | 0 | 28 | 3 | CA3 Glut | 0.013635 | 0.000000e+00 | 0.000000e+00 | -6.017838 | 0.001823 | False |
| 4 | 0 | 28 | 4 | CT SUB Glut | 217.203752 | 0.000000e+00 | 0.000000e+00 | 7.658627 | 0.070052 | True |
Plot the results
[36]:
matrix_df = pd.DataFrame(
data=hf_dist.toarray(),
index=hf_niches,
columns=hf_ct,
)
cn_dist_count = hf_dist.toarray() * cell_count_niche[:, np.newaxis]
cn_dist_norm = cn_dist_count / np.sum(cn_dist_count, axis=0)
matrix_df_norm = pd.DataFrame(
data=cn_dist_norm,
index=hf_niches,
columns=hf_ct,
)
hf_df['stars'] = hf_df['q-value'].apply(p2stars)
stars_df = pd.DataFrame(
'',
index=matrix_df.index,
columns=matrix_df.columns
)
for _, row in hf_df[hf_df['enrichment']].iterrows():
niche = row['niche']
ct = row['celltype']
if (niche in stars_df.index) and (ct in stars_df.columns):
stars_df.loc[niche, ct] = row['stars']
fig, axes = plt.subplots(1, 1, figsize=(18, 3))
sns_heatmap_0 = sns.heatmap(
matrix_df,
cmap='Blues',
# cbar_kws={'label': 'Cell type proportion'},
linewidths=0.5,
linecolor='gray',
# square=True,
ax=axes
)
for i, niche in enumerate(matrix_df.index):
for j, ct in enumerate(matrix_df.columns):
star = stars_df.iloc[i, j]
if star:
if matrix_df.iloc[i, j] > np.max(matrix_df.values) * 0.7:
color='white'
else:
color='black'
axes.text(j + 0.5, i + 0.6, star, ha='center', va='center', color=color, fontsize=20, fontweight='bold')
n_rows, n_cols = matrix_df.shape
axes.plot([0, n_cols], [n_rows, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.plot([n_cols, n_cols], [0, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.set_xticklabels(axes.get_xticklabels(), rotation=90, ha='center', fontsize=20)
axes.set_yticklabels(axes.get_yticklabels(), rotation=0, ha='right', fontsize=20)
axes.set_ylabel('Niche', fontsize=20)
axes.set_xlabel('Cell Type', fontsize=20)
axes.set_title('Cell Type Proportions', fontsize=20)
axes.collections[0].colorbar.ax.yaxis.label.set_size(20)
axes.collections[0].colorbar.ax.tick_params(labelsize=16)
axes.grid(False)
plt.tight_layout()
plt.show()
fig, axes = plt.subplots(1, 1, figsize=(18, 3))
sns_heatmap_1 = sns.heatmap(
matrix_df_norm,
cmap='Blues',
# cbar_kws={'label': 'Cell type proportion'},
linewidths=0.5,
linecolor='gray',
# square=True,
ax=axes
)
for i, niche in enumerate(matrix_df.index):
for j, ct in enumerate(matrix_df.columns):
star = stars_df.iloc[i, j]
if star:
if matrix_df_norm.iloc[i, j] > np.max(matrix_df_norm.values) * 0.7:
color='white'
else:
color='black'
axes.text(j + 0.5, i + 0.6, star, ha='center', va='center', color=color, fontsize=20, fontweight='bold')
n_rows, n_cols = matrix_df.shape
axes.plot([0, n_cols], [n_rows, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.plot([n_cols, n_cols], [0, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.set_xticklabels(axes.get_xticklabels(), rotation=90, ha='center', fontsize=20)
axes.set_yticklabels(axes.get_yticklabels(), rotation=0, ha='right', fontsize=20)
axes.set_ylabel('Niche', fontsize=20)
axes.set_xlabel('Cell Type', fontsize=20)
axes.set_title('Column Normalized Cell Type Proportions', fontsize=20)
axes.collections[0].colorbar.ax.yaxis.label.set_size(20)
axes.collections[0].colorbar.ax.tick_params(labelsize=16)
axes.grid(False)
plt.tight_layout()
plt.show()
Display enriched cell type for each niche
[37]:
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
enriched_ct = ['SUB-ProS Glut', 'CA1-ProS Glut', 'CA2-FC-IG Glut', 'CA3 Glut', 'DG Glut']
print(template)
fig, axes = plt.subplots(2, 3, figsize=(20, 10))
for i, ct in enumerate(enriched_ct):
row_idx = i // 3
col_idx = i % 3
adata.obs['highlight'] = [label if label == ct else 'others' for label in adata.obs['subclass_transfer']]
highlight_dict = {ct: 'brown', 'others': 'lightgray'}
adata_tmp = adata[(adata.obs['highlight'].isin(['others']))].copy()
sc.pl.embedding(adata_tmp, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[row_idx, col_idx], s=70, show=False, frameon=False, legend_loc=None)
adata_tmp = adata[~(adata.obs['highlight'].isin(['others']))].copy()
sc.pl.embedding(adata_tmp, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[row_idx, col_idx], s=70, show=False, frameon=False, legend_loc=None)
axes[row_idx, col_idx].set_title(ct, fontsize=20)
for row in axes:
for ax in row:
ax.set_xlim(max(adata.obsm['spatial'][:, 0])-8150, max(adata.obsm['spatial'][:, 0])-5700)
ax.set_ylim(max(adata.obsm['spatial'][:, 1])-2300, max(adata.obsm['spatial'][:, 1])-900)
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.008
Preprocess data and find differentially expressed genes for each niche
[38]:
# adata_hf = adata_new[adata_new.obs['niche_label'].isin(hf_niches), :].copy()
sc.pp.filter_genes(adata_hf, min_cells=100)
sc.pp.normalize_total(adata_hf, target_sum=1e3)
sc.pp.log1p(adata_hf)
sc.tl.rank_genes_groups(adata_hf, groupby='niche_label', method='wilcoxon', pts=True)
adata_hf
[38]:
AnnData object with n_obs × n_vars = 69208 × 1122
obs: 'donor_id', 'cell_type_ontology_term_id', 'tissue_ontology_term_id', 'subclass_transfer', 'major_brain_region', 'ccf_region_name', 'brain_section_label', 'cell_type', 'tissue', 'development_stage', 'slice_name', 'celltype_idx', 'n_neighbors', 'niche_label_jsd', 'niche_label_jsd_v2', 'niche_label_fmi', 'niche_label_ari', 'niche_label_nmi', 'niche_label_asw', 'niche_label_js_asw', 'niche_label_fisher', 'niche_label_chi', 'niche_label_dbi', 'niche_label_dass_min', 'niche_label_dass_mean', 'niche_label_dafisher', 'niche_label_dachi', 'niche_label_0.09', 'niche_label_0.11', 'niche_label_100', 'niche_label_99', 'niche_label_98', 'niche_label_97', 'niche_label_96', 'niche_label_95', 'niche_label_94', 'niche_label_93', 'niche_label_92', 'niche_label_91', 'niche_label_90', 'niche_label_89', 'niche_label_88', 'niche_label_87', 'niche_label_86', 'niche_label_85', 'niche_label_84', 'niche_label_83', 'niche_label_82', 'niche_label_81', 'niche_label_80', 'niche_label_79', 'niche_label_78', 'niche_label_77', 'niche_label_76', 'niche_label_75', 'niche_label_74', 'niche_label_73', 'niche_label_72', 'niche_label_71', 'niche_label_70', 'niche_label_69', 'niche_label_68', 'niche_label_67', 'niche_label_66', 'niche_label_65', 'niche_label_64', 'niche_label_63', 'niche_label_62', 'niche_label_61', 'niche_label_60', 'niche_label_59', 'niche_label_58', 'niche_label_57', 'niche_label_56', 'niche_label_55', 'niche_label_54', 'niche_label_53', 'niche_label_52', 'niche_label_51', 'niche_label_50', 'niche_label_49', 'niche_label_48', 'niche_label_47', 'niche_label_46', 'niche_label_45', 'niche_label_44', 'niche_label_43', 'niche_label_42', 'niche_label_41', 'niche_label_40', 'niche_label_39', 'niche_label_38', 'niche_label_37', 'niche_label_36', 'niche_label_35', 'niche_label_34', 'niche_label_33', 'niche_label_32', 'niche_label_31', 'niche_label_30', 'niche_label_29', 'niche_label_28', 'niche_label_27', 'niche_label_26', 'niche_label_25', 'niche_label_24', 'niche_label_23', 'niche_label_22', 'niche_label_21', 'niche_label_20', 'niche_label_19', 'niche_label_18', 'niche_label_17', 'niche_label_16', 'niche_label_15', 'niche_label_14', 'niche_label_13', 'niche_label_12', 'niche_label_11', 'niche_label_10', 'niche_label_9', 'niche_label_8', 'niche_label_7', 'niche_label_6', 'niche_label_5', 'niche_label_4', 'niche_label_3', 'niche_label_2', 'niche_label'
var: 'n_cells'
uns: 'ct2idx', 'idx2ct', 'niche_cell_count', 'niche_dist', 'niche_label_summary', 'score_dict', 'log1p', 'rank_genes_groups'
obsm: 'micro_dist', 'onehot', 'spatial', 'spatial_ccf'
[39]:
n_top_genes = 5
result = adata_hf.uns['rank_genes_groups']
groups = result['names'].dtype.names
selected_genes = []
for group in groups:
genes = pd.Index(result["names"][group], name="gene")
pts_series = result["pts"][group]
pts_aligned = pts_series.reindex(genes)
df = pd.DataFrame({
'gene': result['names'][group],
'pvals_adj': result['pvals_adj'][group],
'logfoldchanges': result['logfoldchanges'][group],
'pts': pts_aligned.values,
})
df_filtered = df[(df['pvals_adj'] < 0.05) & (df['logfoldchanges'] > 1) & (df['pts'] > 0.3)]
topn = df_filtered.sort_values('logfoldchanges', ascending=False).head(n_top_genes)
selected_genes.extend(topn['gene'].values.tolist())
Plot the results
[40]:
result = sc.pl.heatmap(adata_hf, var_names=selected_genes, groupby='niche_label', cmap='viridis', use_raw=False, show=False, standard_scale='var')
fig = result['heatmap_ax'].get_figure()
fig.set_size_inches(6, 9)
ax = result['heatmap_ax']
ax.tick_params(axis='x', labelsize=16)
ax = result['groupby_ax']
ax.tick_params(axis='y', labelsize=16)
plt.tight_layout()
plt.show()
[41]:
result = sc.pl.dotplot(adata_hf, var_names=selected_genes, groupby='niche_label', use_raw=False, show=False, standard_scale="var", dot_max=0.6)
fig = result['mainplot_ax'].get_figure()
fig.set_size_inches(8, 3)
ax = result['mainplot_ax']
ax.tick_params(axis='x', labelsize=16)
ax.tick_params(axis='y', labelsize=16)
ax = result['size_legend_ax']
ax.tick_params(axis='x', labelsize=12)
ax = result['color_legend_ax']
ax.tick_params(axis='x', labelsize=12)
plt.tight_layout()
plt.show()
[42]:
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
print(template)
x = adata.obsm['spatial'][:, 0]
y = adata.obsm['spatial'][:, 1]
x_min, x_max = max(adata.obsm['spatial'][:, 0])-8150, max(adata.obsm['spatial'][:, 0])-5700
y_min, y_max = max(adata.obsm['spatial'][:, 1])-2300, max(adata.obsm['spatial'][:, 1])-900
mask = (x >= x_min) & (x <= x_max) & (y >= y_min) & (y <= y_max)
adata_crop = adata[mask].copy()
sc.pp.normalize_total(adata_crop, target_sum=1e3)
sc.pp.log1p(adata_crop)
fig, axes = plt.subplots(4, 5, figsize=(24, 14))
for i, marker in enumerate(selected_genes):
row = i // 5
col = i % 5
sc.pl.embedding(adata_crop, basis='spatial', color=marker,
ax=axes[row, col], s=30, show=False, frameon=False,
cmap='viridis', title=marker)
axes[row, col].set_title(marker, fontsize=20)
for row in axes:
for ax in row:
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.008
Main olfactory bulb
Show relative niches and corresponding annotations
[43]:
template = 'C57BL6J-3.008'
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
mob_niches = ['58', '2', '50', '14',]
print(template)
fig, axes = plt.subplots(1, 5, figsize=(24, 6))
sc.pl.embedding(adata, basis='spatial', color='niche_label', palette=niche_color_dict,
ax=axes[0], s=30, show=False, frameon=False, legend_loc=None)
axes[0].set_title('Main Olfactory Bulb', fontsize=20)
for i, niche in enumerate(mob_niches):
adata.obs['highlight'] = [label if label == niche else 'others' for label in adata.obs['niche_label']]
highlight_dict = {niche: niche_colors[int(niche)], 'others': 'lightgray'}
sc.pl.embedding(adata, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[i+1], s=30, show=False, frameon=False, legend_loc=None)
axes[i+1].set_title(region_dict[niche], fontsize=20)
for ax in axes:
ax.set_xlim(min(adata.obsm['spatial'][:, 0]), min(adata.obsm['spatial'][:, 0])+3200)
ax.set_ylim(min(adata.obsm['spatial'][:, 1])+800, min(adata.obsm['spatial'][:, 1])+4200)
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.008
Extract data of corresponding niches
[44]:
adata_mob = adata_new[adata_new.obs['niche_label'].isin(mob_niches), :].copy()
adata_mob
[44]:
AnnData object with n_obs × n_vars = 203023 × 1122
obs: 'donor_id', 'cell_type_ontology_term_id', 'tissue_ontology_term_id', 'subclass_transfer', 'major_brain_region', 'ccf_region_name', 'brain_section_label', 'cell_type', 'tissue', 'development_stage', 'slice_name', 'celltype_idx', 'n_neighbors', 'niche_label_jsd', 'niche_label_jsd_v2', 'niche_label_fmi', 'niche_label_ari', 'niche_label_nmi', 'niche_label_asw', 'niche_label_js_asw', 'niche_label_fisher', 'niche_label_chi', 'niche_label_dbi', 'niche_label_dass_min', 'niche_label_dass_mean', 'niche_label_dafisher', 'niche_label_dachi', 'niche_label_0.09', 'niche_label_0.11', 'niche_label_100', 'niche_label_99', 'niche_label_98', 'niche_label_97', 'niche_label_96', 'niche_label_95', 'niche_label_94', 'niche_label_93', 'niche_label_92', 'niche_label_91', 'niche_label_90', 'niche_label_89', 'niche_label_88', 'niche_label_87', 'niche_label_86', 'niche_label_85', 'niche_label_84', 'niche_label_83', 'niche_label_82', 'niche_label_81', 'niche_label_80', 'niche_label_79', 'niche_label_78', 'niche_label_77', 'niche_label_76', 'niche_label_75', 'niche_label_74', 'niche_label_73', 'niche_label_72', 'niche_label_71', 'niche_label_70', 'niche_label_69', 'niche_label_68', 'niche_label_67', 'niche_label_66', 'niche_label_65', 'niche_label_64', 'niche_label_63', 'niche_label_62', 'niche_label_61', 'niche_label_60', 'niche_label_59', 'niche_label_58', 'niche_label_57', 'niche_label_56', 'niche_label_55', 'niche_label_54', 'niche_label_53', 'niche_label_52', 'niche_label_51', 'niche_label_50', 'niche_label_49', 'niche_label_48', 'niche_label_47', 'niche_label_46', 'niche_label_45', 'niche_label_44', 'niche_label_43', 'niche_label_42', 'niche_label_41', 'niche_label_40', 'niche_label_39', 'niche_label_38', 'niche_label_37', 'niche_label_36', 'niche_label_35', 'niche_label_34', 'niche_label_33', 'niche_label_32', 'niche_label_31', 'niche_label_30', 'niche_label_29', 'niche_label_28', 'niche_label_27', 'niche_label_26', 'niche_label_25', 'niche_label_24', 'niche_label_23', 'niche_label_22', 'niche_label_21', 'niche_label_20', 'niche_label_19', 'niche_label_18', 'niche_label_17', 'niche_label_16', 'niche_label_15', 'niche_label_14', 'niche_label_13', 'niche_label_12', 'niche_label_11', 'niche_label_10', 'niche_label_9', 'niche_label_8', 'niche_label_7', 'niche_label_6', 'niche_label_5', 'niche_label_4', 'niche_label_3', 'niche_label_2', 'niche_label'
uns: 'ct2idx', 'idx2ct', 'niche_cell_count', 'niche_dist', 'niche_label_summary', 'score_dict'
obsm: 'micro_dist', 'onehot', 'spatial', 'spatial_ccf'
Filter out cell types that comprise less than 500 cell in these niches
[45]:
min_cells = 500
ct_counts = adata_mob.obs['subclass_transfer'].value_counts()
valid_ct = ct_counts[ct_counts >= min_cells].index
adata_mob = adata_mob[adata_mob.obs['subclass_transfer'].isin(valid_ct)].copy()
adata_mob.shape, len(list(set(adata_mob.obs['subclass_transfer'])))
[45]:
((198585, 1122), 20)
Compute the number of cells in each niche
[46]:
ct_onehot = label2onehot(adata_mob.obs['subclass_transfer'])
indices = [mob_niches.index(str(label)) for label in adata_mob.obs['niche_label']]
cn_onehot = index2onehot(indices)
cell_count_niche = cn_onehot.T.sum(axis=1)
mob_dist = (cn_onehot.T @ ct_onehot) / cell_count_niche
cell_count_niche = np.array(cell_count_niche).flatten()
mob_ct = sorted(set(adata_mob.obs['subclass_transfer']))
idx2ct = {str(i): mob_ct[i] for i in range(len(mob_ct))}
cell_count_niche
[46]:
array([ 11561., 46897., 17995., 122132.])
Run cell type enrichment analysis
[47]:
mob_df = ct_enrichment_test(mob_dist,
cell_count_niche,
idx2ct,
mob_niches,
method='fisher',
alpha=0.05,
fdr_method='fdr_by',
log2fc_threshold=1,
prop_threshold=0.01,
verbose=True,
)
mob_df.head()
4 niches and 20 cell types in total.
[47]:
| niche_idx | niche | celltype_idx | celltype | oddsratio | p-value | q-value | log2fc | prop | enrichment | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 58 | 0 | Astro-OLF NN | 0.783555 | 2.155552e-06 | 1.223240e-05 | -0.338580 | 0.033561 | False |
| 1 | 0 | 58 | 1 | Astro-TE NN | 1.036740 | 6.919693e-01 | 1.000000e+00 | 0.064140 | 0.009861 | False |
| 2 | 0 | 58 | 2 | Endo NN | 1.840272 | 4.981661e-52 | 3.733786e-51 | 0.830596 | 0.073610 | False |
| 3 | 0 | 58 | 3 | Microglia NN | 2.573269 | 6.137616e-52 | 4.514993e-51 | 1.334159 | 0.033042 | True |
| 4 | 0 | 58 | 4 | OB Dopa-Gaba | 0.434197 | 5.949043e-42 | 4.219979e-41 | -1.162588 | 0.018597 | False |
Plot the results
[48]:
matrix_df = pd.DataFrame(
data=mob_dist.toarray(),
index=mob_niches,
columns=mob_ct,
)
cn_dist_count = mob_dist.toarray() * cell_count_niche[:, np.newaxis]
cn_dist_norm = cn_dist_count / np.sum(cn_dist_count, axis=0)
matrix_df_norm = pd.DataFrame(
data=cn_dist_norm,
index=mob_niches,
columns=mob_ct,
)
mob_df['stars'] = mob_df['q-value'].apply(p2stars)
stars_df = pd.DataFrame(
'',
index=matrix_df.index,
columns=matrix_df.columns
)
for _, row in mob_df[mob_df['enrichment']].iterrows():
niche = row['niche']
ct = row['celltype']
if (niche in stars_df.index) and (ct in stars_df.columns):
stars_df.loc[niche, ct] = row['stars']
fig, axes = plt.subplots(1, 1, figsize=(20, 3))
sns_heatmap_0 = sns.heatmap(
matrix_df,
cmap='Blues',
# cbar_kws={'label': 'Cell type proportion'},
linewidths=0.5,
linecolor='gray',
# square=True,
ax=axes
)
for i, niche in enumerate(matrix_df.index):
for j, ct in enumerate(matrix_df.columns):
star = stars_df.iloc[i, j]
if star:
if matrix_df.iloc[i, j] > np.max(matrix_df.values) * 0.7:
color='white'
else:
color='black'
axes.text(j + 0.5, i + 0.6, star, ha='center', va='center', color=color, fontsize=20, fontweight='bold')
n_rows, n_cols = matrix_df.shape
axes.plot([0, n_cols], [n_rows, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.plot([n_cols, n_cols], [0, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.set_xticklabels(axes.get_xticklabels(), rotation=90, ha='center', fontsize=20)
axes.set_yticklabels(axes.get_yticklabels(), rotation=0, ha='right', fontsize=20)
axes.set_ylabel('Niche', fontsize=20)
axes.set_xlabel('Cell Type', fontsize=20)
axes.set_title('Cell Type Proportions', fontsize=20)
axes.collections[0].colorbar.ax.yaxis.label.set_size(20)
axes.collections[0].colorbar.ax.tick_params(labelsize=16)
axes.grid(False)
plt.tight_layout()
plt.show()
fig, axes = plt.subplots(1, 1, figsize=(20, 3))
sns_heatmap_1 = sns.heatmap(
matrix_df_norm,
cmap='Blues',
# cbar_kws={'label': 'Cell type proportion'},
linewidths=0.5,
linecolor='gray',
# square=True,
ax=axes
)
for i, niche in enumerate(matrix_df.index):
for j, ct in enumerate(matrix_df.columns):
star = stars_df.iloc[i, j]
if star:
if matrix_df_norm.iloc[i, j] > np.max(matrix_df_norm.values) * 0.7:
color='white'
else:
color='black'
axes.text(j + 0.5, i + 0.6, star, ha='center', va='center', color=color, fontsize=20, fontweight='bold')
n_rows, n_cols = matrix_df.shape
axes.plot([0, n_cols], [n_rows, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.plot([n_cols, n_cols], [0, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.set_xticklabels(axes.get_xticklabels(), rotation=90, ha='center', fontsize=20)
axes.set_yticklabels(axes.get_yticklabels(), rotation=0, ha='right', fontsize=20)
axes.set_ylabel('Niche', fontsize=20)
axes.set_xlabel('Cell Type', fontsize=20)
axes.set_title('Column Normalized Cell Type Proportions', fontsize=20)
axes.collections[0].colorbar.ax.yaxis.label.set_size(20)
axes.collections[0].colorbar.ax.tick_params(labelsize=16)
axes.grid(False)
plt.tight_layout()
plt.show()
Display enriched cell type for each niche
[49]:
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
enriched_ct = ['OEC NN', 'OB Dopa-Gaba', 'OB-out Frmd7 Gaba', 'OB-mi Frmd7 Gaba', 'OB-in Frmd7 Gaba']
print(template)
fig, axes = plt.subplots(1, 5, figsize=(24, 6))
for i, ct in enumerate(enriched_ct):
adata.obs['highlight'] = [label if label == ct else 'others' for label in adata.obs['subclass_transfer']]
highlight_dict = {ct: 'brown', 'others': 'lightgray'}
adata_tmp = adata[(adata.obs['highlight'].isin(['others']))].copy()
sc.pl.embedding(adata_tmp, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[i], s=30, show=False, frameon=False, legend_loc=None)
adata_tmp = adata[~(adata.obs['highlight'].isin(['others']))].copy()
sc.pl.embedding(adata_tmp, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[i], s=30, show=False, frameon=False, legend_loc=None)
axes[i].set_title(ct, fontsize=20)
for ax in axes:
ax.set_xlim(min(adata.obsm['spatial'][:, 0]), min(adata.obsm['spatial'][:, 0])+3200)
ax.set_ylim(min(adata.obsm['spatial'][:, 1])+800, min(adata.obsm['spatial'][:, 1])+4200)
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.008
Preprocess data and find differentially expressed genes for each niche
[50]:
# adata_mob = adata_new[adata_new.obs['niche_label'].isin(mob_niches), :].copy()
sc.pp.filter_genes(adata_mob, min_cells=100)
sc.pp.normalize_total(adata_mob, target_sum=1e3)
sc.pp.log1p(adata_mob)
sc.tl.rank_genes_groups(adata_mob, groupby='niche_label', method='wilcoxon', pts=True)
adata_mob
[50]:
AnnData object with n_obs × n_vars = 198585 × 1122
obs: 'donor_id', 'cell_type_ontology_term_id', 'tissue_ontology_term_id', 'subclass_transfer', 'major_brain_region', 'ccf_region_name', 'brain_section_label', 'cell_type', 'tissue', 'development_stage', 'slice_name', 'celltype_idx', 'n_neighbors', 'niche_label_jsd', 'niche_label_jsd_v2', 'niche_label_fmi', 'niche_label_ari', 'niche_label_nmi', 'niche_label_asw', 'niche_label_js_asw', 'niche_label_fisher', 'niche_label_chi', 'niche_label_dbi', 'niche_label_dass_min', 'niche_label_dass_mean', 'niche_label_dafisher', 'niche_label_dachi', 'niche_label_0.09', 'niche_label_0.11', 'niche_label_100', 'niche_label_99', 'niche_label_98', 'niche_label_97', 'niche_label_96', 'niche_label_95', 'niche_label_94', 'niche_label_93', 'niche_label_92', 'niche_label_91', 'niche_label_90', 'niche_label_89', 'niche_label_88', 'niche_label_87', 'niche_label_86', 'niche_label_85', 'niche_label_84', 'niche_label_83', 'niche_label_82', 'niche_label_81', 'niche_label_80', 'niche_label_79', 'niche_label_78', 'niche_label_77', 'niche_label_76', 'niche_label_75', 'niche_label_74', 'niche_label_73', 'niche_label_72', 'niche_label_71', 'niche_label_70', 'niche_label_69', 'niche_label_68', 'niche_label_67', 'niche_label_66', 'niche_label_65', 'niche_label_64', 'niche_label_63', 'niche_label_62', 'niche_label_61', 'niche_label_60', 'niche_label_59', 'niche_label_58', 'niche_label_57', 'niche_label_56', 'niche_label_55', 'niche_label_54', 'niche_label_53', 'niche_label_52', 'niche_label_51', 'niche_label_50', 'niche_label_49', 'niche_label_48', 'niche_label_47', 'niche_label_46', 'niche_label_45', 'niche_label_44', 'niche_label_43', 'niche_label_42', 'niche_label_41', 'niche_label_40', 'niche_label_39', 'niche_label_38', 'niche_label_37', 'niche_label_36', 'niche_label_35', 'niche_label_34', 'niche_label_33', 'niche_label_32', 'niche_label_31', 'niche_label_30', 'niche_label_29', 'niche_label_28', 'niche_label_27', 'niche_label_26', 'niche_label_25', 'niche_label_24', 'niche_label_23', 'niche_label_22', 'niche_label_21', 'niche_label_20', 'niche_label_19', 'niche_label_18', 'niche_label_17', 'niche_label_16', 'niche_label_15', 'niche_label_14', 'niche_label_13', 'niche_label_12', 'niche_label_11', 'niche_label_10', 'niche_label_9', 'niche_label_8', 'niche_label_7', 'niche_label_6', 'niche_label_5', 'niche_label_4', 'niche_label_3', 'niche_label_2', 'niche_label'
var: 'n_cells'
uns: 'ct2idx', 'idx2ct', 'niche_cell_count', 'niche_dist', 'niche_label_summary', 'score_dict', 'log1p', 'rank_genes_groups'
obsm: 'micro_dist', 'onehot', 'spatial', 'spatial_ccf'
[51]:
n_top_genes = 5
result = adata_mob.uns['rank_genes_groups']
groups = result['names'].dtype.names
selected_genes = []
for group in groups:
genes = pd.Index(result["names"][group], name="gene")
pts_series = result["pts"][group]
pts_aligned = pts_series.reindex(genes)
df = pd.DataFrame({
'gene': result['names'][group],
'pvals_adj': result['pvals_adj'][group],
'logfoldchanges': result['logfoldchanges'][group],
'pts': pts_aligned.values,
})
df_filtered = df[(df['pvals_adj'] < 0.05) & (df['logfoldchanges'] > 1) & (df['pts'] > 0.3)]
topn = df_filtered.sort_values('logfoldchanges', ascending=False).head(n_top_genes)
selected_genes.extend(topn['gene'].values.tolist())
Plot the results
[52]:
result = sc.pl.heatmap(adata_mob, var_names=selected_genes, groupby='niche_label', cmap='viridis', use_raw=False, show=False, standard_scale='var')
fig = result['heatmap_ax'].get_figure()
fig.set_size_inches(6, 9)
ax = result['heatmap_ax']
ax.tick_params(axis='x', labelsize=16)
ax = result['groupby_ax']
ax.tick_params(axis='y', labelsize=16)
plt.tight_layout()
plt.show()
[53]:
result = sc.pl.dotplot(adata_mob, var_names=selected_genes, groupby='niche_label', use_raw=False, show=False, standard_scale="var", dot_max=0.6)
fig = result['mainplot_ax'].get_figure()
fig.set_size_inches(8, 3)
ax = result['mainplot_ax']
ax.tick_params(axis='x', labelsize=16)
ax.tick_params(axis='y', labelsize=16)
ax = result['size_legend_ax']
ax.tick_params(axis='x', labelsize=12)
ax = result['color_legend_ax']
ax.tick_params(axis='x', labelsize=12)
plt.tight_layout()
plt.show()
[54]:
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
print(template)
x = adata.obsm['spatial'][:, 0]
y = adata.obsm['spatial'][:, 1]
x_min, x_max = min(adata.obsm['spatial'][:, 0]), min(adata.obsm['spatial'][:, 0])+3200
y_min, y_max = min(adata.obsm['spatial'][:, 1])+800, min(adata.obsm['spatial'][:, 1])+4200
mask = (x >= x_min) & (x <= x_max) & (y >= y_min) & (y <= y_max)
adata_crop = adata[mask].copy()
sc.pp.normalize_total(adata_crop, target_sum=1e3)
sc.pp.log1p(adata_crop)
fig, axes = plt.subplots(4, 5, figsize=(24, 18))
for i, marker in enumerate(selected_genes):
row = i // 5
col = i % 5
sc.pl.embedding(adata_crop, basis='spatial', color=marker,
ax=axes[row, col], s=20, show=False, frameon=False,
cmap='viridis', title=marker)
axes[row, col].set_title(marker, fontsize=20)
for row in axes:
for ax in row:
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.008
Ventricular systems
Show relative niches and corresponding annotations
[55]:
template = 'C57BL6J-3.008'
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
lv_niches = ['12', '23', '55']
print(template)
fig, axes = plt.subplots(1, 4, figsize=(22, 6))
sc.pl.embedding(adata, basis='spatial', color='niche_label', palette=niche_color_dict,
ax=axes[0], s=100, show=False, frameon=False, legend_loc=None)
axes[0].set_title('Ventricular System', fontsize=20)
for i, niche in enumerate(lv_niches):
adata.obs['highlight'] = [label if label == niche else 'others' for label in adata.obs['niche_label']]
highlight_dict = {niche: niche_colors[int(niche)], 'others': 'lightgray'}
sc.pl.embedding(adata, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[i+1], s=100, show=False, frameon=False, legend_loc=None)
axes[i+1].set_title(region_dict[niche], fontsize=20)
for ax in axes:
ax.set_xlim(max(adata.obsm['spatial'][:, 0])-9200, max(adata.obsm['spatial'][:, 0])-8000)
ax.set_ylim(max(adata.obsm['spatial'][:, 1])-3000, max(adata.obsm['spatial'][:, 1])-1400)
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.008
Extract data of corresponding niches
[56]:
adata_lv = adata_new[adata_new.obs['niche_label'].isin(lv_niches), :].copy()
adata_lv
[56]:
AnnData object with n_obs × n_vars = 41086 × 1122
obs: 'donor_id', 'cell_type_ontology_term_id', 'tissue_ontology_term_id', 'subclass_transfer', 'major_brain_region', 'ccf_region_name', 'brain_section_label', 'cell_type', 'tissue', 'development_stage', 'slice_name', 'celltype_idx', 'n_neighbors', 'niche_label_jsd', 'niche_label_jsd_v2', 'niche_label_fmi', 'niche_label_ari', 'niche_label_nmi', 'niche_label_asw', 'niche_label_js_asw', 'niche_label_fisher', 'niche_label_chi', 'niche_label_dbi', 'niche_label_dass_min', 'niche_label_dass_mean', 'niche_label_dafisher', 'niche_label_dachi', 'niche_label_0.09', 'niche_label_0.11', 'niche_label_100', 'niche_label_99', 'niche_label_98', 'niche_label_97', 'niche_label_96', 'niche_label_95', 'niche_label_94', 'niche_label_93', 'niche_label_92', 'niche_label_91', 'niche_label_90', 'niche_label_89', 'niche_label_88', 'niche_label_87', 'niche_label_86', 'niche_label_85', 'niche_label_84', 'niche_label_83', 'niche_label_82', 'niche_label_81', 'niche_label_80', 'niche_label_79', 'niche_label_78', 'niche_label_77', 'niche_label_76', 'niche_label_75', 'niche_label_74', 'niche_label_73', 'niche_label_72', 'niche_label_71', 'niche_label_70', 'niche_label_69', 'niche_label_68', 'niche_label_67', 'niche_label_66', 'niche_label_65', 'niche_label_64', 'niche_label_63', 'niche_label_62', 'niche_label_61', 'niche_label_60', 'niche_label_59', 'niche_label_58', 'niche_label_57', 'niche_label_56', 'niche_label_55', 'niche_label_54', 'niche_label_53', 'niche_label_52', 'niche_label_51', 'niche_label_50', 'niche_label_49', 'niche_label_48', 'niche_label_47', 'niche_label_46', 'niche_label_45', 'niche_label_44', 'niche_label_43', 'niche_label_42', 'niche_label_41', 'niche_label_40', 'niche_label_39', 'niche_label_38', 'niche_label_37', 'niche_label_36', 'niche_label_35', 'niche_label_34', 'niche_label_33', 'niche_label_32', 'niche_label_31', 'niche_label_30', 'niche_label_29', 'niche_label_28', 'niche_label_27', 'niche_label_26', 'niche_label_25', 'niche_label_24', 'niche_label_23', 'niche_label_22', 'niche_label_21', 'niche_label_20', 'niche_label_19', 'niche_label_18', 'niche_label_17', 'niche_label_16', 'niche_label_15', 'niche_label_14', 'niche_label_13', 'niche_label_12', 'niche_label_11', 'niche_label_10', 'niche_label_9', 'niche_label_8', 'niche_label_7', 'niche_label_6', 'niche_label_5', 'niche_label_4', 'niche_label_3', 'niche_label_2', 'niche_label'
uns: 'ct2idx', 'idx2ct', 'niche_cell_count', 'niche_dist', 'niche_label_summary', 'score_dict'
obsm: 'micro_dist', 'onehot', 'spatial', 'spatial_ccf'
Filter out cell types that comprise less than 100 cell in these niches
[57]:
min_cells = 100
ct_counts = adata_lv.obs['subclass_transfer'].value_counts()
valid_ct = ct_counts[ct_counts >= min_cells].index
adata_lv = adata_lv[adata_lv.obs['subclass_transfer'].isin(valid_ct)].copy()
adata_lv.shape, len(list(set(adata_lv.obs['subclass_transfer'])))
[57]:
((39333, 1122), 24)
Compute the number of cells in each niche
[58]:
ct_onehot = label2onehot(adata_lv.obs['subclass_transfer'])
indices = [lv_niches.index(str(label)) for label in adata_lv.obs['niche_label']]
cn_onehot = index2onehot(indices)
cell_count_niche = cn_onehot.T.sum(axis=1)
lv_dist = (cn_onehot.T @ ct_onehot) / cell_count_niche
cell_count_niche = np.array(cell_count_niche).flatten()
lv_ct = sorted(set(adata_lv.obs['subclass_transfer']))
idx2ct = {str(i): lv_ct[i] for i in range(len(lv_ct))}
cell_count_niche
[58]:
array([10878., 12856., 15599.])
Run cell type enrichment analysis
[59]:
lv_df = ct_enrichment_test(lv_dist,
cell_count_niche,
idx2ct,
lv_niches,
method='fisher',
alpha=0.05,
fdr_method='fdr_by',
log2fc_threshold=1,
prop_threshold=0.01,
verbose=True,
)
lv_df.head()
3 niches and 24 cell types in total.
[59]:
| niche_idx | niche | celltype_idx | celltype | oddsratio | p-value | q-value | log2fc | prop | enrichment | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 12 | 0 | Astro-NT NN | 14.870298 | 7.913046e-279 | 2.130303e-277 | 3.798464 | 0.068946 | True |
| 1 | 0 | 12 | 1 | Astro-OLF NN | 0.013757 | 7.214507e-98 | 1.262461e-96 | -6.148006 | 0.000368 | False |
| 2 | 0 | 12 | 2 | Astro-TE NN | 1.034157 | 3.157359e-01 | 1.000000e+00 | 0.042081 | 0.133480 | False |
| 3 | 0 | 12 | 3 | Astroependymal NN | 2.517070 | 1.279644e-09 | 7.997278e-09 | 1.324533 | 0.008274 | False |
| 4 | 0 | 12 | 4 | BAM NN | 0.320814 | 2.613817e-10 | 1.663235e-09 | -1.632631 | 0.002482 | False |
Plot the results
[60]:
matrix_df = pd.DataFrame(
data=lv_dist.toarray(),
index=lv_niches,
columns=lv_ct,
)
cn_dist_count = lv_dist.toarray() * cell_count_niche[:, np.newaxis]
cn_dist_norm = cn_dist_count / np.sum(cn_dist_count, axis=0)
matrix_df_norm = pd.DataFrame(
data=cn_dist_norm,
index=lv_niches,
columns=lv_ct,
)
lv_df['stars'] = lv_df['q-value'].apply(p2stars)
stars_df = pd.DataFrame(
'',
index=matrix_df.index,
columns=matrix_df.columns
)
for _, row in lv_df[lv_df['enrichment']].iterrows():
niche = row['niche']
ct = row['celltype']
if (niche in stars_df.index) and (ct in stars_df.columns):
stars_df.loc[niche, ct] = row['stars']
fig, axes = plt.subplots(1, 1, figsize=(20, 2))
sns_heatmap_0 = sns.heatmap(
matrix_df,
cmap='Blues',
# cbar_kws={'label': 'Cell type proportion'},
linewidths=0.5,
linecolor='gray',
# square=True,
ax=axes
)
for i, niche in enumerate(matrix_df.index):
for j, ct in enumerate(matrix_df.columns):
star = stars_df.iloc[i, j]
if star:
if matrix_df.iloc[i, j] > np.max(matrix_df.values) * 0.7:
color='white'
else:
color='black'
axes.text(j + 0.5, i + 0.6, star, ha='center', va='center', color=color, fontsize=20, fontweight='bold')
n_rows, n_cols = matrix_df.shape
axes.plot([0, n_cols], [n_rows, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.plot([n_cols, n_cols], [0, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.set_xticklabels(axes.get_xticklabels(), rotation=90, ha='center', fontsize=20)
axes.set_yticklabels(axes.get_yticklabels(), rotation=0, ha='right', fontsize=20)
axes.set_ylabel('Niche', fontsize=20)
axes.set_xlabel('Cell Type', fontsize=20)
axes.set_title('Cell Type Proportions', fontsize=20)
axes.collections[0].colorbar.ax.yaxis.label.set_size(20)
axes.collections[0].colorbar.ax.tick_params(labelsize=16)
axes.grid(False)
plt.tight_layout()
plt.show()
fig, axes = plt.subplots(1, 1, figsize=(20, 2))
sns_heatmap_1 = sns.heatmap(
matrix_df_norm,
cmap='Blues',
# cbar_kws={'label': 'Cell type proportion'},
linewidths=0.5,
linecolor='gray',
# square=True,
ax=axes
)
for i, niche in enumerate(matrix_df.index):
for j, ct in enumerate(matrix_df.columns):
star = stars_df.iloc[i, j]
if star:
if matrix_df_norm.iloc[i, j] > np.max(matrix_df_norm.values) * 0.7:
color='white'
else:
color='black'
axes.text(j + 0.5, i + 0.6, star, ha='center', va='center', color=color, fontsize=20, fontweight='bold')
n_rows, n_cols = matrix_df.shape
axes.plot([0, n_cols], [n_rows, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.plot([n_cols, n_cols], [0, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.set_xticklabels(axes.get_xticklabels(), rotation=90, ha='center', fontsize=20)
axes.set_yticklabels(axes.get_yticklabels(), rotation=0, ha='right', fontsize=20)
axes.set_ylabel('Niche', fontsize=20)
axes.set_xlabel('Cell Type', fontsize=20)
axes.set_title('Column Normalized Cell Type Proportions', fontsize=20)
axes.collections[0].colorbar.ax.yaxis.label.set_size(20)
axes.collections[0].colorbar.ax.tick_params(labelsize=16)
axes.grid(False)
plt.tight_layout()
plt.show()
Display enriched cell type for each niche
[61]:
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
enriched_ct = ['Ependymal NN', 'CHOR NN', 'OB-STR-CTX Inh IMN']
print(template)
fig, axes = plt.subplots(1, 4, figsize=(22, 6))
for i, ct in enumerate(enriched_ct):
adata.obs['highlight'] = [label if label == ct else 'others' for label in adata.obs['subclass_transfer']]
highlight_dict = {ct: 'brown', 'others': 'lightgray'}
adata_tmp = adata[(adata.obs['highlight'].isin(['others']))].copy()
sc.pl.embedding(adata_tmp, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[i], s=100, show=False, frameon=False, legend_loc=None)
adata_tmp = adata[~(adata.obs['highlight'].isin(['others']))].copy()
sc.pl.embedding(adata_tmp, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[i], s=100, show=False, frameon=False, legend_loc=None)
axes[i].set_title(ct, fontsize=20)
for ax in axes:
ax.set_xlim(max(adata.obsm['spatial'][:, 0])-9200, max(adata.obsm['spatial'][:, 0])-8000)
ax.set_ylim(max(adata.obsm['spatial'][:, 1])-3000, max(adata.obsm['spatial'][:, 1])-1400)
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.008
Preprocess data and find differentially expressed genes for each niche
[62]:
# adata_lv = adata_new[adata_new.obs['niche_label'].isin(lv_niches), :].copy()
sc.pp.filter_genes(adata_lv, min_cells=100)
sc.pp.normalize_total(adata_lv, target_sum=1e3)
sc.pp.log1p(adata_lv)
sc.tl.rank_genes_groups(adata_lv, groupby='niche_label', method='wilcoxon', pts=True)
adata_lv
[62]:
AnnData object with n_obs × n_vars = 39333 × 1122
obs: 'donor_id', 'cell_type_ontology_term_id', 'tissue_ontology_term_id', 'subclass_transfer', 'major_brain_region', 'ccf_region_name', 'brain_section_label', 'cell_type', 'tissue', 'development_stage', 'slice_name', 'celltype_idx', 'n_neighbors', 'niche_label_jsd', 'niche_label_jsd_v2', 'niche_label_fmi', 'niche_label_ari', 'niche_label_nmi', 'niche_label_asw', 'niche_label_js_asw', 'niche_label_fisher', 'niche_label_chi', 'niche_label_dbi', 'niche_label_dass_min', 'niche_label_dass_mean', 'niche_label_dafisher', 'niche_label_dachi', 'niche_label_0.09', 'niche_label_0.11', 'niche_label_100', 'niche_label_99', 'niche_label_98', 'niche_label_97', 'niche_label_96', 'niche_label_95', 'niche_label_94', 'niche_label_93', 'niche_label_92', 'niche_label_91', 'niche_label_90', 'niche_label_89', 'niche_label_88', 'niche_label_87', 'niche_label_86', 'niche_label_85', 'niche_label_84', 'niche_label_83', 'niche_label_82', 'niche_label_81', 'niche_label_80', 'niche_label_79', 'niche_label_78', 'niche_label_77', 'niche_label_76', 'niche_label_75', 'niche_label_74', 'niche_label_73', 'niche_label_72', 'niche_label_71', 'niche_label_70', 'niche_label_69', 'niche_label_68', 'niche_label_67', 'niche_label_66', 'niche_label_65', 'niche_label_64', 'niche_label_63', 'niche_label_62', 'niche_label_61', 'niche_label_60', 'niche_label_59', 'niche_label_58', 'niche_label_57', 'niche_label_56', 'niche_label_55', 'niche_label_54', 'niche_label_53', 'niche_label_52', 'niche_label_51', 'niche_label_50', 'niche_label_49', 'niche_label_48', 'niche_label_47', 'niche_label_46', 'niche_label_45', 'niche_label_44', 'niche_label_43', 'niche_label_42', 'niche_label_41', 'niche_label_40', 'niche_label_39', 'niche_label_38', 'niche_label_37', 'niche_label_36', 'niche_label_35', 'niche_label_34', 'niche_label_33', 'niche_label_32', 'niche_label_31', 'niche_label_30', 'niche_label_29', 'niche_label_28', 'niche_label_27', 'niche_label_26', 'niche_label_25', 'niche_label_24', 'niche_label_23', 'niche_label_22', 'niche_label_21', 'niche_label_20', 'niche_label_19', 'niche_label_18', 'niche_label_17', 'niche_label_16', 'niche_label_15', 'niche_label_14', 'niche_label_13', 'niche_label_12', 'niche_label_11', 'niche_label_10', 'niche_label_9', 'niche_label_8', 'niche_label_7', 'niche_label_6', 'niche_label_5', 'niche_label_4', 'niche_label_3', 'niche_label_2', 'niche_label'
var: 'n_cells'
uns: 'ct2idx', 'idx2ct', 'niche_cell_count', 'niche_dist', 'niche_label_summary', 'score_dict', 'log1p', 'rank_genes_groups'
obsm: 'micro_dist', 'onehot', 'spatial', 'spatial_ccf'
[63]:
n_top_genes = 5
result = adata_lv.uns['rank_genes_groups']
groups = result['names'].dtype.names
selected_genes = []
for group in groups:
genes = pd.Index(result["names"][group], name="gene")
pts_series = result["pts"][group]
pts_aligned = pts_series.reindex(genes)
df = pd.DataFrame({
'gene': result['names'][group],
'pvals_adj': result['pvals_adj'][group],
'logfoldchanges': result['logfoldchanges'][group],
'pts': pts_aligned.values,
})
df_filtered = df[(df['pvals_adj'] < 0.05) & (df['logfoldchanges'] > 1) & (df['pts'] > 0.3)]
topn = df_filtered.sort_values('logfoldchanges', ascending=False).head(n_top_genes)
selected_genes.extend(topn['gene'].values.tolist())
Plot the results
[64]:
result = sc.pl.heatmap(adata_lv, var_names=selected_genes, groupby='niche_label', cmap='viridis', use_raw=False, show=False, standard_scale='var')
fig = result['heatmap_ax'].get_figure()
fig.set_size_inches(6, 9)
ax = result['heatmap_ax']
ax.tick_params(axis='x', labelsize=16)
ax = result['groupby_ax']
ax.tick_params(axis='y', labelsize=16)
plt.tight_layout()
plt.show()
[65]:
result = sc.pl.dotplot(adata_lv, var_names=selected_genes, groupby='niche_label', use_raw=False, show=False, standard_scale="var", dot_max=0.6)
fig = result['mainplot_ax'].get_figure()
fig.set_size_inches(8, 3)
ax = result['mainplot_ax']
ax.tick_params(axis='x', labelsize=16)
ax.tick_params(axis='y', labelsize=16)
ax = result['size_legend_ax']
ax.tick_params(axis='x', labelsize=12)
ax = result['color_legend_ax']
ax.tick_params(axis='x', labelsize=12)
plt.tight_layout()
plt.show()
[66]:
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
print(template)
x = adata.obsm['spatial'][:, 0]
y = adata.obsm['spatial'][:, 1]
x_min, x_max = max(adata.obsm['spatial'][:, 0])-9200, max(adata.obsm['spatial'][:, 0])-8000
y_min, y_max = max(adata.obsm['spatial'][:, 1])-3000, max(adata.obsm['spatial'][:, 1])-1400
mask = (x >= x_min) & (x <= x_max) & (y >= y_min) & (y <= y_max)
adata_crop = adata[mask].copy()
sc.pp.normalize_total(adata_crop, target_sum=1e3)
sc.pp.log1p(adata_crop)
fig, axes = plt.subplots(3, 5, figsize=(24, 14))
for i, marker in enumerate(selected_genes):
row = i // 5
col = i % 5
sc.pl.embedding(adata_crop, basis='spatial', color=marker,
ax=axes[row, col], s=50, show=False, frameon=False,
cmap='viridis', title=marker)
axes[row, col].set_title(marker, fontsize=20)
for row in axes:
for ax in row:
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.008
Midbrain
Show relative niches and corresponding annotations
[67]:
template = 'C57BL6J-3.003'
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
mid_niches = ['4', '15', '26', '35', '48', '51']
print(template)
fig, axes = plt.subplots(2, 4, figsize=(20, 9))
axes = axes.flatten()
sc.pl.embedding(adata, basis='spatial', color='niche_label', palette=niche_color_dict,
ax=axes[0], s=40, show=False, frameon=False, legend_loc=None)
axes[0].set_title('Midbrain', fontsize=20)
for i, niche in enumerate(mid_niches):
adata.obs['highlight'] = [label if label == niche else 'others' for label in adata.obs['niche_label']]
highlight_dict = {niche: niche_colors[int(niche)], 'others': 'lightgray'}
sc.pl.embedding(adata, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[i+1], s=40, show=False, frameon=False, legend_loc=None)
axes[i+1].set_title(region_dict[niche], fontsize=20)
for ax in axes:
ax.set_xlim(max(adata.obsm['spatial'][:, 0])-7000, max(adata.obsm['spatial'][:, 0])-3300)
ax.set_ylim(max(adata.obsm['spatial'][:, 1])-3700, max(adata.obsm['spatial'][:, 1])-400)
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.003
Extract data of corresponding niches
[68]:
adata_mid = adata_new[adata_new.obs['niche_label'].isin(mid_niches), :].copy()
adata_mid
[68]:
AnnData object with n_obs × n_vars = 86678 × 1122
obs: 'donor_id', 'cell_type_ontology_term_id', 'tissue_ontology_term_id', 'subclass_transfer', 'major_brain_region', 'ccf_region_name', 'brain_section_label', 'cell_type', 'tissue', 'development_stage', 'slice_name', 'celltype_idx', 'n_neighbors', 'niche_label_jsd', 'niche_label_jsd_v2', 'niche_label_fmi', 'niche_label_ari', 'niche_label_nmi', 'niche_label_asw', 'niche_label_js_asw', 'niche_label_fisher', 'niche_label_chi', 'niche_label_dbi', 'niche_label_dass_min', 'niche_label_dass_mean', 'niche_label_dafisher', 'niche_label_dachi', 'niche_label_0.09', 'niche_label_0.11', 'niche_label_100', 'niche_label_99', 'niche_label_98', 'niche_label_97', 'niche_label_96', 'niche_label_95', 'niche_label_94', 'niche_label_93', 'niche_label_92', 'niche_label_91', 'niche_label_90', 'niche_label_89', 'niche_label_88', 'niche_label_87', 'niche_label_86', 'niche_label_85', 'niche_label_84', 'niche_label_83', 'niche_label_82', 'niche_label_81', 'niche_label_80', 'niche_label_79', 'niche_label_78', 'niche_label_77', 'niche_label_76', 'niche_label_75', 'niche_label_74', 'niche_label_73', 'niche_label_72', 'niche_label_71', 'niche_label_70', 'niche_label_69', 'niche_label_68', 'niche_label_67', 'niche_label_66', 'niche_label_65', 'niche_label_64', 'niche_label_63', 'niche_label_62', 'niche_label_61', 'niche_label_60', 'niche_label_59', 'niche_label_58', 'niche_label_57', 'niche_label_56', 'niche_label_55', 'niche_label_54', 'niche_label_53', 'niche_label_52', 'niche_label_51', 'niche_label_50', 'niche_label_49', 'niche_label_48', 'niche_label_47', 'niche_label_46', 'niche_label_45', 'niche_label_44', 'niche_label_43', 'niche_label_42', 'niche_label_41', 'niche_label_40', 'niche_label_39', 'niche_label_38', 'niche_label_37', 'niche_label_36', 'niche_label_35', 'niche_label_34', 'niche_label_33', 'niche_label_32', 'niche_label_31', 'niche_label_30', 'niche_label_29', 'niche_label_28', 'niche_label_27', 'niche_label_26', 'niche_label_25', 'niche_label_24', 'niche_label_23', 'niche_label_22', 'niche_label_21', 'niche_label_20', 'niche_label_19', 'niche_label_18', 'niche_label_17', 'niche_label_16', 'niche_label_15', 'niche_label_14', 'niche_label_13', 'niche_label_12', 'niche_label_11', 'niche_label_10', 'niche_label_9', 'niche_label_8', 'niche_label_7', 'niche_label_6', 'niche_label_5', 'niche_label_4', 'niche_label_3', 'niche_label_2', 'niche_label'
uns: 'ct2idx', 'idx2ct', 'niche_cell_count', 'niche_dist', 'niche_label_summary', 'score_dict'
obsm: 'micro_dist', 'onehot', 'spatial', 'spatial_ccf'
Filter out cell types that comprise less than 500 cell in these niches
[69]:
min_cells = 500
ct_counts = adata_mid.obs['subclass_transfer'].value_counts()
valid_ct = ct_counts[ct_counts >= min_cells].index
adata_mid = adata_mid[adata_mid.obs['subclass_transfer'].isin(valid_ct)].copy()
adata_mid.shape, len(list(set(adata_mid.obs['subclass_transfer'])))
[69]:
((79578, 1122), 30)
Compute the number of cells in each niche
[70]:
ct_onehot = label2onehot(adata_mid.obs['subclass_transfer'])
indices = [mid_niches.index(str(label)) for label in adata_mid.obs['niche_label']]
cn_onehot = index2onehot(indices)
cell_count_niche = cn_onehot.T.sum(axis=1)
mid_dist = (cn_onehot.T @ ct_onehot) / cell_count_niche
cell_count_niche = np.array(cell_count_niche).flatten()
mid_ct = sorted(set(adata_mid.obs['subclass_transfer']))
idx2ct = {str(i): mid_ct[i] for i in range(len(mid_ct))}
cell_count_niche
[70]:
array([ 9318., 14076., 13817., 27257., 4697., 10413.])
Run cell type enrichment analysis
[71]:
mid_df = ct_enrichment_test(mid_dist,
cell_count_niche,
idx2ct,
mid_niches,
method='fisher',
alpha=0.05,
fdr_method='fdr_by',
log2fc_threshold=1,
prop_threshold=0.01,
verbose=True,
)
mid_df.head()
6 niches and 30 cell types in total.
[71]:
| niche_idx | niche | celltype_idx | celltype | oddsratio | p-value | q-value | log2fc | prop | enrichment | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 4 | 0 | APN C1ql2 Glut | 0.000000 | 3.307996e-53 | 4.841464e-52 | -27.036237 | 0.000000 | False |
| 1 | 0 | 4 | 1 | APN C1ql4 Glut | 0.037798 | 9.886438e-46 | 1.317090e-44 | -4.705975 | 0.000537 | False |
| 2 | 0 | 4 | 2 | Astro-NT NN | 1.397420 | 2.220405e-29 | 2.480958e-28 | 0.406630 | 0.182443 | False |
| 3 | 0 | 4 | 3 | Endo NN | 0.759335 | 1.163221e-13 | 1.051077e-12 | -0.356740 | 0.089719 | False |
| 4 | 0 | 4 | 4 | IC Six3 En2 Gaba | 0.019660 | 5.589897e-75 | 9.681054e-74 | -5.638057 | 0.000429 | False |
Plot the results
[72]:
matrix_df = pd.DataFrame(
data=mid_dist.toarray(),
index=mid_niches,
columns=mid_ct,
)
cn_dist_count = mid_dist.toarray() * cell_count_niche[:, np.newaxis]
cn_dist_norm = cn_dist_count / np.sum(cn_dist_count, axis=0)
matrix_df_norm = pd.DataFrame(
data=cn_dist_norm,
index=mid_niches,
columns=mid_ct,
)
mid_df['stars'] = mid_df['q-value'].apply(p2stars)
stars_df = pd.DataFrame(
'',
index=matrix_df.index,
columns=matrix_df.columns
)
for _, row in mid_df[mid_df['enrichment']].iterrows():
niche = row['niche']
ct = row['celltype']
if (niche in stars_df.index) and (ct in stars_df.columns):
stars_df.loc[niche, ct] = row['stars']
fig, axes = plt.subplots(1, 1, figsize=(20, 3))
sns_heatmap_0 = sns.heatmap(
matrix_df,
cmap='Blues',
# cbar_kws={'label': 'Cell type proportion'},
linewidths=0.5,
linecolor='gray',
# square=True,
ax=axes
)
for i, niche in enumerate(matrix_df.index):
for j, ct in enumerate(matrix_df.columns):
star = stars_df.iloc[i, j]
if star:
if matrix_df.iloc[i, j] > np.max(matrix_df.values) * 0.7:
color='white'
else:
color='black'
axes.text(j + 0.5, i + 0.6, star, ha='center', va='center', color=color, fontsize=20, fontweight='bold')
n_rows, n_cols = matrix_df.shape
axes.plot([0, n_cols], [n_rows, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.plot([n_cols, n_cols], [0, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.set_xticklabels(axes.get_xticklabels(), rotation=90, ha='center', fontsize=20)
axes.set_yticklabels(axes.get_yticklabels(), rotation=0, ha='right', fontsize=20)
axes.set_ylabel('Niche', fontsize=20)
axes.set_xlabel('Cell Type', fontsize=20)
axes.set_title('Cell Type Proportions', fontsize=20)
axes.collections[0].colorbar.ax.yaxis.label.set_size(20)
axes.collections[0].colorbar.ax.tick_params(labelsize=16)
axes.grid(False)
plt.tight_layout()
plt.show()
fig, axes = plt.subplots(1, 1, figsize=(20, 3))
sns_heatmap_1 = sns.heatmap(
matrix_df_norm,
cmap='Blues',
# cbar_kws={'label': 'Cell type proportion'},
linewidths=0.5,
linecolor='gray',
# square=True,
ax=axes
)
for i, niche in enumerate(matrix_df.index):
for j, ct in enumerate(matrix_df.columns):
star = stars_df.iloc[i, j]
if star:
if matrix_df_norm.iloc[i, j] > np.max(matrix_df_norm.values) * 0.7:
color='white'
else:
color='black'
axes.text(j + 0.5, i + 0.6, star, ha='center', va='center', color=color, fontsize=20, fontweight='bold')
n_rows, n_cols = matrix_df.shape
axes.plot([0, n_cols], [n_rows, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.plot([n_cols, n_cols], [0, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.set_xticklabels(axes.get_xticklabels(), rotation=90, ha='center', fontsize=20)
axes.set_yticklabels(axes.get_yticklabels(), rotation=0, ha='right', fontsize=20)
axes.set_ylabel('Niche', fontsize=20)
axes.set_xlabel('Cell Type', fontsize=20)
axes.set_title('Column Normalized Cell Type Proportions', fontsize=20)
axes.collections[0].colorbar.ax.yaxis.label.set_size(20)
axes.collections[0].colorbar.ax.tick_params(labelsize=16)
axes.grid(False)
plt.tight_layout()
plt.show()
Display enriched cell type for each niche
[73]:
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
enriched_ct = ['SCsg Pde5a Glut', 'SCsg Gabrr2 Gaba', 'SCig Foxb1 Glut', 'PRT Tcf7l2 Gaba',
'IC Tfap2d Maf Glut', 'SCig-an-PPT Foxb1 Glut', 'PAG Pou4f2 Glut']
print(template)
fig, axes = plt.subplots(2, 4, figsize=(20, 9))
axes = axes.flatten()
for i, ct in enumerate(enriched_ct):
adata.obs['highlight'] = [label if label == ct else 'others' for label in adata.obs['subclass_transfer']]
highlight_dict = {ct: 'brown', 'others': 'lightgray'}
adata_tmp = adata[(adata.obs['highlight'].isin(['others']))].copy()
sc.pl.embedding(adata_tmp, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[i], s=30, show=False, frameon=False, legend_loc=None)
adata_tmp = adata[~(adata.obs['highlight'].isin(['others']))].copy()
sc.pl.embedding(adata_tmp, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[i], s=30, show=False, frameon=False, legend_loc=None)
axes[i].set_title(ct, fontsize=20)
for ax in axes:
ax.set_xlim(max(adata.obsm['spatial'][:, 0])-7000, max(adata.obsm['spatial'][:, 0])-3300)
ax.set_ylim(max(adata.obsm['spatial'][:, 1])-3700, max(adata.obsm['spatial'][:, 1])-400)
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.003
Hypothalamus
Show relative niches and corresponding annotations
[74]:
template = 'C57BL6J-3.001'
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
hth_niches = ['8', '27', '34', '40', '47', '49', '59', '61']
print(template)
fig, axes = plt.subplots(3, 3, figsize=(24, 15))
axes = axes.flatten()
sc.pl.embedding(adata, basis='spatial', color='niche_label', palette=niche_color_dict,
ax=axes[0], s=80, show=False, frameon=False, legend_loc=None)
axes[0].set_title('Hypothalamus', fontsize=20)
for i, niche in enumerate(hth_niches):
adata.obs['highlight'] = [label if label == niche else 'others' for label in adata.obs['niche_label']]
highlight_dict = {niche: niche_colors[int(niche)], 'others': 'lightgray'}
sc.pl.embedding(adata, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[i+1], s=80, show=False, frameon=False, legend_loc=None)
axes[i+1].set_title(region_dict[niche], fontsize=20)
for ax in axes:
ax.set_xlim(max(adata.obsm['spatial'][:, 0])-8500, max(adata.obsm['spatial'][:, 0])-4400)
ax.set_ylim(max(adata.obsm['spatial'][:, 1])-6000, max(adata.obsm['spatial'][:, 1])-3800)
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.001
Extract data of corresponding niches
[75]:
adata_hth = adata_new[adata_new.obs['niche_label'].isin(hth_niches), :].copy()
adata_hth
[75]:
AnnData object with n_obs × n_vars = 99219 × 1122
obs: 'donor_id', 'cell_type_ontology_term_id', 'tissue_ontology_term_id', 'subclass_transfer', 'major_brain_region', 'ccf_region_name', 'brain_section_label', 'cell_type', 'tissue', 'development_stage', 'slice_name', 'celltype_idx', 'n_neighbors', 'niche_label_jsd', 'niche_label_jsd_v2', 'niche_label_fmi', 'niche_label_ari', 'niche_label_nmi', 'niche_label_asw', 'niche_label_js_asw', 'niche_label_fisher', 'niche_label_chi', 'niche_label_dbi', 'niche_label_dass_min', 'niche_label_dass_mean', 'niche_label_dafisher', 'niche_label_dachi', 'niche_label_0.09', 'niche_label_0.11', 'niche_label_100', 'niche_label_99', 'niche_label_98', 'niche_label_97', 'niche_label_96', 'niche_label_95', 'niche_label_94', 'niche_label_93', 'niche_label_92', 'niche_label_91', 'niche_label_90', 'niche_label_89', 'niche_label_88', 'niche_label_87', 'niche_label_86', 'niche_label_85', 'niche_label_84', 'niche_label_83', 'niche_label_82', 'niche_label_81', 'niche_label_80', 'niche_label_79', 'niche_label_78', 'niche_label_77', 'niche_label_76', 'niche_label_75', 'niche_label_74', 'niche_label_73', 'niche_label_72', 'niche_label_71', 'niche_label_70', 'niche_label_69', 'niche_label_68', 'niche_label_67', 'niche_label_66', 'niche_label_65', 'niche_label_64', 'niche_label_63', 'niche_label_62', 'niche_label_61', 'niche_label_60', 'niche_label_59', 'niche_label_58', 'niche_label_57', 'niche_label_56', 'niche_label_55', 'niche_label_54', 'niche_label_53', 'niche_label_52', 'niche_label_51', 'niche_label_50', 'niche_label_49', 'niche_label_48', 'niche_label_47', 'niche_label_46', 'niche_label_45', 'niche_label_44', 'niche_label_43', 'niche_label_42', 'niche_label_41', 'niche_label_40', 'niche_label_39', 'niche_label_38', 'niche_label_37', 'niche_label_36', 'niche_label_35', 'niche_label_34', 'niche_label_33', 'niche_label_32', 'niche_label_31', 'niche_label_30', 'niche_label_29', 'niche_label_28', 'niche_label_27', 'niche_label_26', 'niche_label_25', 'niche_label_24', 'niche_label_23', 'niche_label_22', 'niche_label_21', 'niche_label_20', 'niche_label_19', 'niche_label_18', 'niche_label_17', 'niche_label_16', 'niche_label_15', 'niche_label_14', 'niche_label_13', 'niche_label_12', 'niche_label_11', 'niche_label_10', 'niche_label_9', 'niche_label_8', 'niche_label_7', 'niche_label_6', 'niche_label_5', 'niche_label_4', 'niche_label_3', 'niche_label_2', 'niche_label'
uns: 'ct2idx', 'idx2ct', 'niche_cell_count', 'niche_dist', 'niche_label_summary', 'score_dict'
obsm: 'micro_dist', 'onehot', 'spatial', 'spatial_ccf'
Filter out cell types that comprise less than 500 cell in these niches
[76]:
min_cells = 500
ct_counts = adata_hth.obs['subclass_transfer'].value_counts()
valid_ct = ct_counts[ct_counts >= min_cells].index
adata_hth = adata_hth[adata_hth.obs['subclass_transfer'].isin(valid_ct)].copy()
adata_hth.shape, len(list(set(adata_hth.obs['subclass_transfer'])))
[76]:
((82698, 1122), 35)
Compute the number of cells in each niche
[77]:
ct_onehot = label2onehot(adata_hth.obs['subclass_transfer'])
indices = [hth_niches.index(str(label)) for label in adata_hth.obs['niche_label']]
cn_onehot = index2onehot(indices)
cell_count_niche = cn_onehot.T.sum(axis=1)
hth_dist = (cn_onehot.T @ ct_onehot) / cell_count_niche
cell_count_niche = np.array(cell_count_niche).flatten()
hth_ct = sorted(set(adata_hth.obs['subclass_transfer']))
idx2ct = {str(i): hth_ct[i] for i in range(len(hth_ct))}
cell_count_niche
[77]:
array([ 5043., 3898., 7237., 3272., 12874., 6625., 41467., 2282.])
Run cell type enrichment analysis
[78]:
hth_df = ct_enrichment_test(hth_dist,
cell_count_niche,
idx2ct,
hth_niches,
method='fisher',
alpha=0.05,
fdr_method='fdr_by',
log2fc_threshold=1,
prop_threshold=0.01,
verbose=True,
)
hth_df.head()
8 niches and 35 cell types in total.
[78]:
| niche_idx | niche | celltype_idx | celltype | oddsratio | p-value | q-value | log2fc | prop | enrichment | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 8 | 0 | AHN Onecut3 Gaba | 0.406090 | 2.187859e-05 | 1.730259e-04 | -1.292618 | 0.003569 | False |
| 1 | 0 | 8 | 1 | AHN-RCH-LHA Otp Fezf1 Glut | 0.194680 | 9.176093e-16 | 1.191427e-14 | -2.245100 | 0.002974 | False |
| 2 | 0 | 8 | 2 | AHN-SBPV-PVHd Pdrm12 Gaba | 0.299865 | 4.671227e-11 | 5.048004e-10 | -1.723654 | 0.004164 | False |
| 3 | 0 | 8 | 3 | Astro-NT NN | 0.776393 | 3.591985e-12 | 4.138778e-11 | -0.289933 | 0.185802 | False |
| 4 | 0 | 8 | 4 | BST-MPN Six3 Nrgn Gaba | 0.549145 | 5.213317e-03 | 3.672246e-02 | -0.860052 | 0.003966 | False |
Plot the results
[79]:
matrix_df = pd.DataFrame(
data=hth_dist.toarray(),
index=hth_niches,
columns=hth_ct,
)
cn_dist_count = hth_dist.toarray() * cell_count_niche[:, np.newaxis]
cn_dist_norm = cn_dist_count / np.sum(cn_dist_count, axis=0)
matrix_df_norm = pd.DataFrame(
data=cn_dist_norm,
index=hth_niches,
columns=hth_ct,
)
hth_df['stars'] = hth_df['q-value'].apply(p2stars)
stars_df = pd.DataFrame(
'',
index=matrix_df.index,
columns=matrix_df.columns
)
for _, row in hth_df[hth_df['enrichment']].iterrows():
niche = row['niche']
ct = row['celltype']
if (niche in stars_df.index) and (ct in stars_df.columns):
stars_df.loc[niche, ct] = row['stars']
fig, axes = plt.subplots(1, 1, figsize=(24, 9))
sns_heatmap_0 = sns.heatmap(
matrix_df,
cmap='Blues',
# cbar_kws={'label': 'Cell type proportion'},
linewidths=0.5,
linecolor='gray',
# square=True,
ax=axes
)
for i, niche in enumerate(matrix_df.index):
for j, ct in enumerate(matrix_df.columns):
star = stars_df.iloc[i, j]
if star:
if matrix_df.iloc[i, j] > np.max(matrix_df.values) * 0.7:
color='white'
else:
color='black'
axes.text(j + 0.5, i + 0.6, star, ha='center', va='center', color=color, fontsize=20, fontweight='bold')
n_rows, n_cols = matrix_df.shape
axes.plot([0, n_cols], [n_rows, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.plot([n_cols, n_cols], [0, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.set_xticklabels(axes.get_xticklabels(), rotation=90, ha='center', fontsize=20)
axes.set_yticklabels(axes.get_yticklabels(), rotation=0, ha='right', fontsize=20)
axes.set_ylabel('Niche', fontsize=20)
axes.set_xlabel('Cell Type', fontsize=20)
axes.set_title('Cell Type Proportions', fontsize=20)
axes.collections[0].colorbar.ax.yaxis.label.set_size(20)
axes.collections[0].colorbar.ax.tick_params(labelsize=16)
axes.grid(False)
plt.tight_layout()
plt.show()
fig, axes = plt.subplots(1, 1, figsize=(24, 9))
sns_heatmap_1 = sns.heatmap(
matrix_df_norm,
cmap='Blues',
# cbar_kws={'label': 'Cell type proportion'},
linewidths=0.5,
linecolor='gray',
# square=True,
ax=axes
)
for i, niche in enumerate(matrix_df.index):
for j, ct in enumerate(matrix_df.columns):
star = stars_df.iloc[i, j]
if star:
if matrix_df_norm.iloc[i, j] > np.max(matrix_df_norm.values) * 0.7:
color='white'
else:
color='black'
axes.text(j + 0.5, i + 0.6, star, ha='center', va='center', color=color, fontsize=20, fontweight='bold')
n_rows, n_cols = matrix_df.shape
axes.plot([0, n_cols], [n_rows, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.plot([n_cols, n_cols], [0, n_rows], color='gray', linewidth=0.5, clip_on=False)
axes.set_xticklabels(axes.get_xticklabels(), rotation=90, ha='center', fontsize=20)
axes.set_yticklabels(axes.get_yticklabels(), rotation=0, ha='right', fontsize=20)
axes.set_ylabel('Niche', fontsize=20)
axes.set_xlabel('Cell Type', fontsize=20)
axes.set_title('Column Normalized Cell Type Proportions', fontsize=20)
axes.collections[0].colorbar.ax.yaxis.label.set_size(20)
axes.collections[0].colorbar.ax.tick_params(labelsize=16)
axes.grid(False)
plt.tight_layout()
plt.show()
Display enriched cell type for each niche
[80]:
adata = adata_new[adata_new.obs['slice_name'] == template, :].copy()
enriched_ct = ['DMH Hmx2 Gaba',
'MPO-ADP Lhx8 Gaba', 'PVpo-VMPO-MPN Hmx2 Gaba',
'PH Pitx2 Glut',
'VMH Nr5a1 Glut',
'ZI Pax6 Gaba',
'AHN Onecut3 Gaba', 'AHN-RCH-LHA Otp Fezf1 Glut', 'AHN-SBPV-PVHd Pdrm12 Gaba',
'MM Foxb1 Glut',
'BST-MPN Six3 Nrgn Gaba', 'MEA-BST Lhx6 Nfib Gaba',]
print(template)
fig, axes = plt.subplots(4, 3, figsize=(24, 20))
axes = axes.flatten()
for i, ct in enumerate(enriched_ct):
adata.obs['highlight'] = [label if label == ct else 'others' for label in adata.obs['subclass_transfer']]
highlight_dict = {ct: 'brown', 'others': 'lightgray'}
adata_tmp = adata[(adata.obs['highlight'].isin(['others']))].copy()
sc.pl.embedding(adata_tmp, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[i], s=80, show=False, frameon=False, legend_loc=None)
adata_tmp = adata[~(adata.obs['highlight'].isin(['others']))].copy()
sc.pl.embedding(adata_tmp, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[i], s=80, show=False, frameon=False, legend_loc=None)
axes[i].set_title(ct, fontsize=20)
for ax in axes:
ax.set_xlim(max(adata.obsm['spatial'][:, 0])-8500, max(adata.obsm['spatial'][:, 0])-4400)
ax.set_ylim(max(adata.obsm['spatial'][:, 1])-6000, max(adata.obsm['spatial'][:, 1])-3800)
ax.axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.001
Non-neuronal cells in major brain regions
[82]:
n_ct = len(sorted(set(adata_new.obs['subclass_transfer'])))
mbr = sorted(set(adata_new.obs['major_brain_region']))
n_niches = len(mbr)
niche_dist, niche_cell_count = calculate_distribution(adata_new.obs['major_brain_region'].tolist(),
adata_new.obs['celltype_idx'].tolist(),
label_summary=mbr,
n_niches=n_niches,
n_celltypes=n_ct,
change2str=True,
sparse=True,
)
adata_new.uns['niche_dist_mbr'] = niche_dist
adata_new.uns['niche_cell_count_mbr'] = niche_cell_count
adata_new.uns['niche_label_summary_mbr'] = mbr
[83]:
ct_df = ct_enrichment_test(adata_new.uns['niche_dist_mbr'],
adata_new.uns['niche_cell_count_mbr'],
adata_new.uns['idx2ct'],
adata_new.uns['niche_label_summary_mbr'],
method='fisher',
alpha=0.05,
fdr_method='fdr_by',
log2fc_threshold=1,
prop_threshold=0.01,
verbose=True,
)
ct_df.head()
15 niches and 338 cell types in total.
[83]:
| niche_idx | niche | celltype_idx | celltype | oddsratio | p-value | q-value | log2fc | prop | enrichment | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | Cerebellum | 0 | ABC NN | 0.868828 | 4.040883e-05 | 5.461122e-04 | -0.202176 | 0.003132 | False |
| 1 | 0 | Cerebellum | 1 | ACB-BST-FS D1 Gaba | 0.112662 | 2.853784e-37 | 8.060349e-36 | -3.149330 | 0.000052 | False |
| 2 | 0 | Cerebellum | 2 | AD Serpinb7 Glut | 0.000000 | 9.522971e-22 | 2.090149e-20 | -20.726676 | 0.000000 | False |
| 3 | 0 | Cerebellum | 3 | ADP-MPO Trp73 Glut | 0.026801 | 4.937920e-27 | 1.191386e-25 | -5.221203 | 0.000007 | False |
| 4 | 0 | Cerebellum | 4 | AHN Onecut3 Gaba | 0.000000 | 6.243249e-58 | 2.219484e-56 | -22.160329 | 0.000000 | False |
[84]:
ct_labels = sorted(adata_new.obs['subclass_transfer'].unique())
# cn_dist_count = adata_new.uns['niche_dist_mbr'].toarray() * adata_new.uns['niche_cell_count_mbr'][:, np.newaxis]
# cn_dist_norm = cn_dist_count / np.sum(cn_dist_count, axis=0)
matrix_df_norm = pd.DataFrame(
data=adata_new.uns['niche_dist_mbr'].toarray().T,
index=ct_labels,
columns=mbr,
)
ct_df['stars'] = ct_df['q-value'].apply(p2stars)
stars_df = pd.DataFrame(
'',
index=matrix_df_norm.index,
columns=matrix_df_norm.columns
)
for _, row in ct_df[ct_df['enrichment']].iterrows():
niche = row['niche']
ct = row['celltype']
if (ct in stars_df.index) and (niche in stars_df.columns):
stars_df.loc[ct, niche] = row['stars']
fig, ax = plt.subplots(1, 1, figsize=(10, 80))
sns_heatmap = sns.heatmap(
matrix_df_norm,
cmap='Blues',
# cbar_kws={'label': 'Cell type proportion'},
linewidths=0.5,
linecolor='gray',
# square=True,
ax=ax,
cbar_kws={
'label': 'Cell type proportion',
'shrink': 0.4,
'aspect': 50
}
)
for i, ct in enumerate(matrix_df_norm.index):
for j, niche in enumerate(matrix_df_norm.columns):
star = stars_df.iloc[i, j]
if star:
if matrix_df_norm.iloc[i, j] > np.max(matrix_df_norm.values) * 0.7:
color='white'
else:
color='black'
ax.text(j + 0.5, i + 0.6, star, ha='center', va='center', color=color, fontsize=16, fontweight='bold')
n_rows, n_cols = matrix_df_norm.shape
ax.plot([0, n_cols], [n_rows, n_rows], color='gray', linewidth=0.5, clip_on=False)
ax.plot([n_cols, n_cols], [0, n_rows], color='gray', linewidth=0.5, clip_on=False)
ax.set_xticks(np.arange(len(matrix_df_norm.columns)) + 0.5)
ax.set_xticklabels(matrix_df_norm.columns, rotation=90, ha='center', fontsize=16)
ax.set_yticks(np.arange(len(matrix_df_norm.index)) + 0.5)
ax.set_yticklabels(matrix_df_norm.index, rotation=0, ha='right', fontsize=16)
ax.set_ylabel('Cell Type', fontsize=16)
ax.set_xlabel('Major Brain Region', fontsize=16)
ax.set_title('Column Normalized Cell Type Proportions', fontsize=16)
ax.collections[0].colorbar.ax.yaxis.label.set_size(16)
ax.collections[0].colorbar.ax.tick_params(labelsize=14)
ax.grid(False)
plt.tight_layout()
plt.show()
[85]:
dist_mbr = adata_new.uns['niche_dist_mbr'].toarray()
ct2idx = {ct: i for i, ct in enumerate(ct_labels)}
ct_oligo = 'Oligo NN'
ct_astro = 'Astro-NT NN'
ct_endo = 'Endo NN'
special = {ct_oligo, ct_astro, ct_endo}
idx_oligo = int(ct2idx[ct_oligo])
idx_astro = int(ct2idx[ct_astro])
idx_endo = int(ct2idx[ct_endo])
nn_other_cts = [ct for ct in ct_labels if ct.endswith(' NN') and ct not in special]
idx_nn_other = [int(ct2idx[ct]) for ct in nn_other_cts]
neuronal_cts = [ct for ct in ct_labels if not ct.endswith(' NN')]
idx_neuronal = [int(ct2idx[ct]) for ct in neuronal_cts]
parts = {
'Oligo NN': [idx_oligo],
'Astro-NT NN': [idx_astro],
'Endo NN': [idx_endo],
'Other NN': idx_nn_other,
'Neuronal': idx_neuronal,
}
prop = pd.DataFrame(index=list(mbr))
for name, idxs in parts.items():
prop[name] = dist_mbr[:, idxs].sum(axis=1) if len(idxs) > 0 else 0.0
prop.head()
[85]:
| Oligo NN | Astro-NT NN | Endo NN | Other NN | Neuronal | |
|---|---|---|---|---|---|
| Cerebellum | 0.036677 | 0.014696 | 0.050525 | 0.121680 | 0.776423 |
| Cortical_subplate | 0.053915 | 0.007892 | 0.085002 | 0.238079 | 0.615112 |
| Fiber_tracts | 0.360179 | 0.041909 | 0.062665 | 0.190795 | 0.344452 |
| Hippocampus | 0.052771 | 0.006340 | 0.097792 | 0.272723 | 0.570373 |
| Hypothalamus | 0.159600 | 0.142769 | 0.091421 | 0.097755 | 0.508455 |
[86]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.cluster.hierarchy as sch
import scipy.spatial.distance as ssd
from matplotlib.gridspec import GridSpec
ctoi_list = ['Oligo NN', 'Astro-NT NN', 'Endo NN', 'Other NN', 'Neuronal']
color_dict = {
'Oligo NN': "#FD6C6F",
'Astro-NT NN': "#19BBDF",
'Endo NN': "#F7F70A",
'Other NN': "#6EE75B",
'Neuronal': "#D3D3D3",
}
X = prop[ctoi_list].values
dist = ssd.pdist(X, metric='cosine')
link = sch.linkage(dist, method='average')
link_opt = sch.optimal_leaf_ordering(link, ssd.squareform(dist))
d = sch.dendrogram(link_opt, no_plot=True)
region_order = prop.index[d['leaves']]
prop_ord = prop.loc[region_order, ctoi_list]
fig = plt.figure(figsize=(12, 6))
gs = GridSpec(2, 1, height_ratios=[1.5, 3], hspace=0.05)
ax_d = fig.add_subplot(gs[0, 0])
sch.dendrogram(
link_opt,
orientation='top',
no_labels=True,
ax=ax_d,
color_threshold=0,
above_threshold_color='black'
)
for spine in ax_d.spines.values():
spine.set_visible(False)
ax_d.set_xticks([])
ax_d.set_yticks([])
ax_d.set_ylim(0, ax_d.get_ylim()[1] * 1.6)
ax = fig.add_subplot(gs[1, 0])
n = prop_ord.shape[0]
x = np.arange(n) * 10 + 5
bottom = np.zeros(n)
for col in ctoi_list:
vals = prop_ord[col].values
ax.bar(
x, vals, bottom=bottom,
width=8,
color=color_dict[col],
edgecolor='none',
label=col
)
bottom += vals
ax.set_xticks(x)
ax.set_xticklabels(prop_ord.index.tolist(), rotation=90, ha='center', fontsize=14)
ax.set_ylabel('Proportion', fontsize=14)
ax.set_xlabel('Major Brain Region', fontsize=14)
ax.set_ylim(0, 1.0)
for spine in ax.spines.values():
spine.set_visible(False)
ax.tick_params(axis='both', length=0)
ax.grid(False)
ax.set_xlim(ax_d.get_xlim())
ax.legend(frameon=False, fontsize=14, loc='upper left', bbox_to_anchor=(1.02, 1.0))
plt.tight_layout()
plt.show()
[87]:
celltype_col = "subclass_transfer"
s = adata_new.obs[celltype_col].astype(str)
mask_nn = s.str.endswith(" NN")
nn = s[mask_nn]
cnt = nn.value_counts()
nn_prop = cnt / cnt.sum()
df_nn = pd.DataFrame({
"celltype_NN": cnt.index,
"count": cnt.values,
"prop_of_all_NN": nn_prop.values
})
df_nn.head(5)
[87]:
| celltype_NN | count | prop_of_all_NN | |
|---|---|---|---|
| 0 | Oligo NN | 249026 | 0.281224 |
| 1 | Endo NN | 178248 | 0.201295 |
| 2 | Astro-TE NN | 110937 | 0.125281 |
| 3 | Astro-NT NN | 82078 | 0.092690 |
| 4 | Microglia NN | 57906 | 0.065393 |
[88]:
ctoi_list = ['Oligo NN', 'Astro-NT NN', 'Endo NN']
highlight_dict = {
'Oligo NN': "#FD6C6F",
'Astro-NT NN': "#19BBDF",
'Endo NN': "#F7F70A",
'others': "#D3D3D3",
}
for i in range(len(slice_name_list)):
print(slice_name_list[i])
adata = adata_new[adata_new.obs['slice_name'] == slice_name_list[i], :].copy()
fig, axes = plt.subplots(1, 2, figsize=(24, 6))
sc.pl.embedding(adata, basis='spatial', color='niche_label', palette=niche_color_dict,
ax=axes[0], s=10, show=False, frameon=False, title='Cell Niche', legend_fontsize=16)
axes[0].set_title('Cell Niche', fontsize=20)
adata.obs['highlight'] = [label if label in ctoi_list else 'others' for label in adata.obs['subclass_transfer']]
adata_tmp = adata[(adata.obs['highlight'].isin(['others']))].copy()
sc.pl.embedding(adata_tmp, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[1], s=10, show=False, frameon=False,)
adata_tmp = adata[~(adata.obs['highlight'].isin(['others']))].copy()
sc.pl.embedding(adata_tmp, basis='spatial', color='highlight', palette=highlight_dict,
ax=axes[1], s=10, show=False, frameon=False,)
axes[1].set_title('Top 3 enriched NN in pons and medulla', fontsize=20)
x_mean = (max(adata.obsm['spatial'][:, 0]) + min(adata.obsm['spatial'][:, 0])) / 2
y_mean = (max(adata.obsm['spatial'][:, 1]) + min(adata.obsm['spatial'][:, 1])) / 2
if x_mean-7800 > min(adata.obsm['spatial'][:, 0]) or x_mean+7800 < max(adata.obsm['spatial'][:, 0]):
raise ValueError('Please adjust the x_lim !')
if y_mean-3600 > min(adata.obsm['spatial'][:, 1]) or y_mean+3600 < max(adata.obsm['spatial'][:, 1]):
raise ValueError('Please adjust the y_lim !')
for ax in axes:
ax.set_xlim(x_mean-7800, x_mean+7800)
ax.set_ylim(y_mean-3600, y_mean+3600)
niche_legend_elements = [
Line2D([0], [0], marker='o', color='w', label=label,
markerfacecolor=color, markersize=10)
for label, color in niche_color_dict.items()
]
axes[0].legend(handles=niche_legend_elements, loc=(1.05, -0.1), frameon=False, ncol=3)
axes[0].axis('off')
plt.tight_layout()
plt.show()
C57BL6J-3.001
C57BL6J-3.002
C57BL6J-3.003
C57BL6J-3.004
C57BL6J-3.005
C57BL6J-3.006
C57BL6J-3.007
C57BL6J-3.008
C57BL6J-3.009
C57BL6J-3.010
C57BL6J-3.011
C57BL6J-3.012
C57BL6J-3.013
C57BL6J-3.015
C57BL6J-3.016
C57BL6J-3.017
C57BL6J-3.019
C57BL6J-3.020
C57BL6J-3.021
C57BL6J-3.022
C57BL6J-3.023
C57BL6J-3.024
C57BL6J-3.025
