Exploration of hormonal receptor genes in organoid atlas dataset - HNOCA extended¶

Reference paper

1. Environment Set Up¶

1.1 Library upload¶

In [1]:
import numpy as np
import pandas as pd
import scanpy as sc
import seaborn as sns
import igraph as ig
import matplotlib.pyplot as plt 
from scipy.sparse import csr_matrix, isspmatrix
from datetime import datetime
import sys
sys.path.append('../')
import functions as fn

print(np.__version__)
print(pd.__version__)
print(sc.__version__)

1.23.5
2.0.0
1.9.3
In [2]:
sc.settings.verbosity = 3
sc.settings.set_figure_params(dpi=100)

1.2 Starting computations: timestamp¶

In [3]:
print(datetime.now())

2026-03-11 12:54:18.288343

2. Read input files¶

2.1 adata loading¶

Adata subset in SubsetHNOCA.ipynb by protol and organoid age

In [4]:
adata = sc.read("../../../../DataDir/ExternalData/SingleCellData/HNOCA_protocol_age.h5ad")
In [5]:
adata
Out[5]:
AnnData object with n_obs × n_vars = 336027 × 35725
    obs: 'assay_differentiation', 'assay_type_differentiation', 'bio_sample', 'cell_line', 'cell_type_original', 'gm', 'id', 'individual', 'state_exact', 'suspension_type', 'tech_sample', 'treatment', 'organoid_age_days', 'publication', 'doi', 'batch', 'annot_level_1', 'annot_level_2', 'annot_level_3_rev2', 'annot_level_4_rev2', 'annot_region_rev2', 'annot_ntt_rev2', 'Hallmark_Glycolysis', 'hnoca_core', 'annot_level_2_extended', 'tissue_type', 'sex_ontology_term_id', 'donor_id', 'assay_ontology_term_id', 'self_reported_ethnicity_ontology_term_id', 'tissue_ontology_term_id', 'disease_ontology_term_id', 'development_stage_ontology_term_id', 'cell_type_ontology_term_id', 'is_primary_data', 'cell_type', 'assay', 'disease', 'sex', 'tissue', 'self_reported_ethnicity', 'development_stage', 'observation_joinid'
    var: 'gene_length', 'highly_variable', 'highly_variable_rank', 'highly_variable_nbatches', 'feature_is_filtered', 'feature_name', 'feature_reference', 'feature_biotype', 'feature_length', 'feature_type'
    uns: 'batch_condition', 'citation', 'default_embedding', 'organism', 'organism_ontology_term_id', 'schema_reference', 'schema_version', 'title'
    obsm: 'X_scpoli', 'X_umap_scpoli'
    obsp: 'knn_scpoli_connectivities', 'knn_scpoli_distances'
In [6]:
adata.var
Out[6]:
gene_length highly_variable highly_variable_rank highly_variable_nbatches feature_is_filtered feature_name feature_reference feature_biotype feature_length feature_type
ENSG00000000003 3796 False 2243.0 118 False TSPAN6 NCBITaxon:9606 gene 2396 protein_coding
ENSG00000000005 1205 True 738.5 114 False TNMD NCBITaxon:9606 gene 873 protein_coding
ENSG00000000419 3004 False 2032.0 9 False DPM1 NCBITaxon:9606 gene 1262 protein_coding
ENSG00000000457 6308 False 2532.0 15 False SCYL3 NCBITaxon:9606 gene 2916 protein_coding
ENSG00000000460 4355 False 2444.0 37 False FIRRM NCBITaxon:9606 gene 2661 protein_coding
... ... ... ... ... ... ... ... ... ... ...
ENSG00000288721 6172 False NaN 0 False ENSG00000288721 NCBITaxon:9606 gene 3684 protein_coding
ENSG00000288722 1707 False 2566.5 6 False F8A1 NCBITaxon:9606 gene 1707 protein_coding
ENSG00000288723 1015 False NaN 0 False ENSG00000288723 NCBITaxon:9606 gene 542 lncRNA
ENSG00000288724 625 False NaN 0 False ENSG00000288724 NCBITaxon:9606 gene 494 lncRNA
ENSG00000288725 3810 False NaN 0 False ENSG00000288725 NCBITaxon:9606 gene 3810 protein_coding

35725 rows × 10 columns

In [7]:
adata.obs.columns
Out[7]:
Index(['assay_differentiation', 'assay_type_differentiation', 'bio_sample',
       'cell_line', 'cell_type_original', 'gm', 'id', 'individual',
       'state_exact', 'suspension_type', 'tech_sample', 'treatment',
       'organoid_age_days', 'publication', 'doi', 'batch', 'annot_level_1',
       'annot_level_2', 'annot_level_3_rev2', 'annot_level_4_rev2',
       'annot_region_rev2', 'annot_ntt_rev2', 'Hallmark_Glycolysis',
       'hnoca_core', 'annot_level_2_extended', 'tissue_type',
       'sex_ontology_term_id', 'donor_id', 'assay_ontology_term_id',
       'self_reported_ethnicity_ontology_term_id', 'tissue_ontology_term_id',
       'disease_ontology_term_id', 'development_stage_ontology_term_id',
       'cell_type_ontology_term_id', 'is_primary_data', 'cell_type', 'assay',
       'disease', 'sex', 'tissue', 'self_reported_ethnicity',
       'development_stage', 'observation_joinid'],
      dtype='object')

2.2 Receptors signature loading¶

Loading of hormonal receptor gene signature.

In [8]:
signatures = '../../../../DataDir/ExternalData/Receptors/ReceptorsComplete.txt'
In [9]:
sig = pd.read_csv(signatures, sep="\t", keep_default_na=False)  #keep_default_na=False: remove Na values
print(sig.shape)
sig

(39, 2)
Out[9]:
GeneName Signature
0 THRB Thyroid
1 THRA Thyroid
2 THRAP3 Thyroid
3 DIO1 Thyroid
4 DIO2 Thyroid
5 DIO3 Thyroid
6 SLC16A10 Thyroid
7 SLC16A2 Thyroid
8 SLC7A5 Thyroid
9 KLF9 Thyroid
10 THRSP Thyroid
11 ESRRG Estrogen
12 ESRRA Estrogen
13 GPER1 Estrogen
14 ESR1 Estrogen
15 ESR2 Estrogen
16 ESRRB Estrogen
17 CYP19A1 Estrogen
18 AR Androgen
19 RBP4 Retinoic Acid
20 RARA Retinoic Acid
21 RARB Retinoic Acid
22 RARG Retinoic Acid
23 RXRA Retinoic Acid
24 RXRB Retinoic Acid
25 RXRG Retinoic Acid
26 AHR AhHyd
27 NR3C1 GC
28 NR1H2 LivX
29 NR1H3 LivX
30 PTGER1 PGE2
31 PTGER2 PGE2
32 PTGER3 PGE2
33 PTGER4 PGE2
34 PPARA PPAR
35 PPARD PPAR
36 PPARG PPAR
37 PGR Progesterone
38 VDR Vitamine D
In [10]:
genes = sig["GeneName"].values.tolist()

3. Visualizations¶

3.1 Counts from adata¶

In [11]:
adata.obs
Out[11]:
assay_differentiation assay_type_differentiation bio_sample cell_line cell_type_original gm id individual state_exact suspension_type ... cell_type_ontology_term_id is_primary_data cell_type assay disease sex tissue self_reported_ethnicity development_stage observation_joinid
new_index
homosapiens_None_2020_10x3v2_bhaduriaparna_001_d10_1038_s41586_020_1962_0_1_H28126SWeek3_TCGGTAAAGGCATGTG Velasco, 2019 (doi: 10.1038/s41586-019-1289-x) guided homosapiens_None_2020_10x3v2_bhaduriaparna_001... custom_H28126 Radial Glia_Hindbrain RG unknown homosapiens_None_2020_10x3v2_bhaduriaparna_001... unknown unknown cell ... CL:0000681 True radial glial cell 10x 3' v2 normal male telencephalon unknown unknown <t`y&lfVy(
homosapiens_None_2020_10x3v2_bhaduriaparna_001_d10_1038_s41586_020_1962_0_33_H28126SWeek3_ATTCTACGTCATATGC Velasco, 2019 (doi: 10.1038/s41586-019-1289-x) guided homosapiens_None_2020_10x3v2_bhaduriaparna_001... custom_H28126 Radial Glia_Hindbrain RG unknown homosapiens_None_2020_10x3v2_bhaduriaparna_001... unknown unknown cell ... CL:0000681 True radial glial cell 10x 3' v2 normal male telencephalon unknown unknown LdHroRYE4j
homosapiens_None_2020_10x3v2_bhaduriaparna_001_d10_1038_s41586_020_1962_0_54_H28126SWeek3_AAAGATGGTTACGACT Velasco, 2019 (doi: 10.1038/s41586-019-1289-x) guided homosapiens_None_2020_10x3v2_bhaduriaparna_001... custom_H28126 Mixed Neuron/Radial Glia_Mixed unknown homosapiens_None_2020_10x3v2_bhaduriaparna_001... unknown unknown cell ... CL:0002319 True neural cell 10x 3' v2 normal male telencephalon unknown unknown yf?|Ue-li!
homosapiens_None_2020_10x3v2_bhaduriaparna_001_d10_1038_s41586_020_1962_0_97_H28126SWeek3_GCGCAGTTCATGTAGC Velasco, 2019 (doi: 10.1038/s41586-019-1289-x) guided homosapiens_None_2020_10x3v2_bhaduriaparna_001... custom_H28126 Radial Glia_Pan-radial glia unknown homosapiens_None_2020_10x3v2_bhaduriaparna_001... unknown unknown cell ... CL:0000681 True radial glial cell 10x 3' v2 normal male telencephalon unknown unknown `hGx_2?OL{
homosapiens_None_2020_10x3v2_bhaduriaparna_001_d10_1038_s41586_020_1962_0_100_H28126SWeek3_CCAGCGATCATCTGTT Velasco, 2019 (doi: 10.1038/s41586-019-1289-x) guided homosapiens_None_2020_10x3v2_bhaduriaparna_001... custom_H28126 Radial Glia_Pan-radial glia unknown homosapiens_None_2020_10x3v2_bhaduriaparna_001... unknown unknown cell ... CL:0000681 True radial glial cell 10x 3' v2 normal male telencephalon unknown unknown n;^Bhpg>KM
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
homosapiens_telencephalon_2022_10x3v3_uzquianoana_003_d10_1016_j_cell_2022_09_010_38570_3_TTTGTTGGTGATGAAT-1_2_1.5m Velasco, 2019 (doi: 10.1038/s41586-019-1289-x) guided homosapiens_telencephalon_2022_10x3v3_uzquiano... GM23338 aRG unknown homosapiens_telencephalon_2022_10x3v3_uzquiano... homosapiens_telencephalon_2022_10x3v3_uzquiano... unknown cell ... CL:0000681 True radial glial cell 10x 3' v3 normal male telencephalon unknown 55-year-old stage v$^8QG1sL8
homosapiens_telencephalon_2022_10x3v3_uzquianoana_003_d10_1016_j_cell_2022_09_010_38571_3_TTTGTTGTCAAGCCAT-1_2_1.5m Velasco, 2019 (doi: 10.1038/s41586-019-1289-x) guided homosapiens_telencephalon_2022_10x3v3_uzquiano... GM23338 Immature DL PN unknown homosapiens_telencephalon_2022_10x3v3_uzquiano... homosapiens_telencephalon_2022_10x3v3_uzquiano... unknown cell ... CL:0000598 True pyramidal neuron 10x 3' v3 normal male telencephalon unknown 55-year-old stage g@%ZeUo^AT
homosapiens_telencephalon_2022_10x3v3_uzquianoana_003_d10_1016_j_cell_2022_09_010_38572_3_TTTGTTGTCAGAGCGA-1_2_1.5m Velasco, 2019 (doi: 10.1038/s41586-019-1289-x) guided homosapiens_telencephalon_2022_10x3v3_uzquiano... GM23338 Immature DL PN unknown homosapiens_telencephalon_2022_10x3v3_uzquiano... homosapiens_telencephalon_2022_10x3v3_uzquiano... unknown cell ... CL:0000598 True pyramidal neuron 10x 3' v3 normal male telencephalon unknown 55-year-old stage lt!(%sk)!3
homosapiens_telencephalon_2022_10x3v3_uzquianoana_003_d10_1016_j_cell_2022_09_010_38573_3_TTTGTTGTCTAGACCA-1_2_1.5m Velasco, 2019 (doi: 10.1038/s41586-019-1289-x) guided homosapiens_telencephalon_2022_10x3v3_uzquiano... GM23338 aRG unknown homosapiens_telencephalon_2022_10x3v3_uzquiano... homosapiens_telencephalon_2022_10x3v3_uzquiano... unknown cell ... CL:0000681 True radial glial cell 10x 3' v3 normal male telencephalon unknown 55-year-old stage >q;KQf*b8g
homosapiens_telencephalon_2022_10x3v3_uzquianoana_003_d10_1016_j_cell_2022_09_010_38574_3_TTTGTTGTCTGACAGT-1_2_1.5m Velasco, 2019 (doi: 10.1038/s41586-019-1289-x) guided homosapiens_telencephalon_2022_10x3v3_uzquiano... GM23338 CFuPN unknown homosapiens_telencephalon_2022_10x3v3_uzquiano... homosapiens_telencephalon_2022_10x3v3_uzquiano... unknown cell ... CL:4023009 True extratelencephalic-projecting glutamatergic co... 10x 3' v3 normal male telencephalon unknown 55-year-old stage Z+NPwXB(&`

336027 rows × 43 columns

In [12]:
adata.obsm
Out[12]:
AxisArrays with keys: X_scpoli, X_umap_scpoli

3.2 Clusters annotation¶

In [13]:
sc.pl.embedding(adata,  basis="X_umap_scpoli", color=['cell_type', 'organoid_age_days', 'assay_differentiation', 'tissue'], ncols=1)

/usr/local/lib/python3.8/dist-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
/usr/local/lib/python3.8/dist-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
/usr/local/lib/python3.8/dist-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
In [ ]:
#sc.tl.draw_graph(adata, random_state=1, neighbors_key="harmony")
In [ ]:
#adata.obsm["X_draw_graph_fa_harmony"] = adata.obsm["X_draw_graph_fa"].copy()
#del adata.obsm["X_draw_graph_fa"]
In [14]:
print(adata.X.data[:10])

[6.0723834 6.0723834 6.0723834 6.0723834 6.764377  6.0723834 6.0723834
 6.0723834 6.0723834 6.0723834]
In [15]:
adata.raw.X.data[:20]   # se sparse
Out[15]:
array([1., 1., 1., 1., 2., 1., 1., 1., 1., 1., 2., 1., 1., 7., 1., 1., 1.,
       1., 1., 2.], dtype=float32)
In [16]:
adata.X = adata.raw.X.copy()
In [17]:
adata.layers['counts'] = adata.X.copy()
In [19]:
adata.obs['assay_differentiation'].value_counts()
Out[19]:
assay_differentiation
Velasco, 2019 (doi: 10.1038/s41586-019-1289-x)                   277411
Bhaduri, 2020 (doi: 10.1038/s41586-020-1962-0); most directed     28522
Bhaduri, 2020 (doi: 10.1038/s41586-020-1962-0); directed          13272
Pellegrini, 2020 (doi: 10.1126/science.aaz5626); hChPO             9637
Trujillo, 2019 (doi: 10.1016/j.stem.2019.08.002)                   4427
Pellegrini, 2020 (doi: 10.1126/science.aaz5626); hCO               2758
Name: count, dtype: int64
In [21]:
# Velasco
adata_velasco = adata[adata.obs['assay_differentiation'].str.contains("Velasco")].copy()
adata_velasco
Out[21]:
AnnData object with n_obs × n_vars = 277411 × 35725
    obs: 'assay_differentiation', 'assay_type_differentiation', 'bio_sample', 'cell_line', 'cell_type_original', 'gm', 'id', 'individual', 'state_exact', 'suspension_type', 'tech_sample', 'treatment', 'organoid_age_days', 'publication', 'doi', 'batch', 'annot_level_1', 'annot_level_2', 'annot_level_3_rev2', 'annot_level_4_rev2', 'annot_region_rev2', 'annot_ntt_rev2', 'Hallmark_Glycolysis', 'hnoca_core', 'annot_level_2_extended', 'tissue_type', 'sex_ontology_term_id', 'donor_id', 'assay_ontology_term_id', 'self_reported_ethnicity_ontology_term_id', 'tissue_ontology_term_id', 'disease_ontology_term_id', 'development_stage_ontology_term_id', 'cell_type_ontology_term_id', 'is_primary_data', 'cell_type', 'assay', 'disease', 'sex', 'tissue', 'self_reported_ethnicity', 'development_stage', 'observation_joinid'
    var: 'gene_length', 'highly_variable', 'highly_variable_rank', 'highly_variable_nbatches', 'feature_is_filtered', 'feature_name', 'feature_reference', 'feature_biotype', 'feature_length', 'feature_type'
    uns: 'batch_condition', 'citation', 'default_embedding', 'organism', 'organism_ontology_term_id', 'schema_reference', 'schema_version', 'title', 'cell_type_colors', 'assay_differentiation_colors', 'tissue_colors'
    obsm: 'X_scpoli', 'X_umap_scpoli'
    layers: 'counts'
    obsp: 'knn_scpoli_connectivities', 'knn_scpoli_distances'
In [22]:
# Bhaduri (both directed and most directed)
adata_bhaduri = adata[adata.obs['assay_differentiation'].str.contains("Bhaduri")].copy()
adata_bhaduri
Out[22]:
AnnData object with n_obs × n_vars = 41794 × 35725
    obs: 'assay_differentiation', 'assay_type_differentiation', 'bio_sample', 'cell_line', 'cell_type_original', 'gm', 'id', 'individual', 'state_exact', 'suspension_type', 'tech_sample', 'treatment', 'organoid_age_days', 'publication', 'doi', 'batch', 'annot_level_1', 'annot_level_2', 'annot_level_3_rev2', 'annot_level_4_rev2', 'annot_region_rev2', 'annot_ntt_rev2', 'Hallmark_Glycolysis', 'hnoca_core', 'annot_level_2_extended', 'tissue_type', 'sex_ontology_term_id', 'donor_id', 'assay_ontology_term_id', 'self_reported_ethnicity_ontology_term_id', 'tissue_ontology_term_id', 'disease_ontology_term_id', 'development_stage_ontology_term_id', 'cell_type_ontology_term_id', 'is_primary_data', 'cell_type', 'assay', 'disease', 'sex', 'tissue', 'self_reported_ethnicity', 'development_stage', 'observation_joinid'
    var: 'gene_length', 'highly_variable', 'highly_variable_rank', 'highly_variable_nbatches', 'feature_is_filtered', 'feature_name', 'feature_reference', 'feature_biotype', 'feature_length', 'feature_type'
    uns: 'batch_condition', 'citation', 'default_embedding', 'organism', 'organism_ontology_term_id', 'schema_reference', 'schema_version', 'title', 'cell_type_colors', 'assay_differentiation_colors', 'tissue_colors'
    obsm: 'X_scpoli', 'X_umap_scpoli'
    layers: 'counts'
    obsp: 'knn_scpoli_connectivities', 'knn_scpoli_distances'
In [23]:
# Trujillo
adata_trujillo = adata[adata.obs['assay_differentiation'].str.contains("Trujillo")].copy()
adata_trujillo
Out[23]:
AnnData object with n_obs × n_vars = 4427 × 35725
    obs: 'assay_differentiation', 'assay_type_differentiation', 'bio_sample', 'cell_line', 'cell_type_original', 'gm', 'id', 'individual', 'state_exact', 'suspension_type', 'tech_sample', 'treatment', 'organoid_age_days', 'publication', 'doi', 'batch', 'annot_level_1', 'annot_level_2', 'annot_level_3_rev2', 'annot_level_4_rev2', 'annot_region_rev2', 'annot_ntt_rev2', 'Hallmark_Glycolysis', 'hnoca_core', 'annot_level_2_extended', 'tissue_type', 'sex_ontology_term_id', 'donor_id', 'assay_ontology_term_id', 'self_reported_ethnicity_ontology_term_id', 'tissue_ontology_term_id', 'disease_ontology_term_id', 'development_stage_ontology_term_id', 'cell_type_ontology_term_id', 'is_primary_data', 'cell_type', 'assay', 'disease', 'sex', 'tissue', 'self_reported_ethnicity', 'development_stage', 'observation_joinid'
    var: 'gene_length', 'highly_variable', 'highly_variable_rank', 'highly_variable_nbatches', 'feature_is_filtered', 'feature_name', 'feature_reference', 'feature_biotype', 'feature_length', 'feature_type'
    uns: 'batch_condition', 'citation', 'default_embedding', 'organism', 'organism_ontology_term_id', 'schema_reference', 'schema_version', 'title', 'cell_type_colors', 'assay_differentiation_colors', 'tissue_colors'
    obsm: 'X_scpoli', 'X_umap_scpoli'
    layers: 'counts'
    obsp: 'knn_scpoli_connectivities', 'knn_scpoli_distances'
In [24]:
# Pellegrini (hChPO and hCO)
adata_pellegrini = adata[adata.obs['assay_differentiation'].str.contains("Pellegrini")].copy()
adata_pellegrini
Out[24]:
AnnData object with n_obs × n_vars = 12395 × 35725
    obs: 'assay_differentiation', 'assay_type_differentiation', 'bio_sample', 'cell_line', 'cell_type_original', 'gm', 'id', 'individual', 'state_exact', 'suspension_type', 'tech_sample', 'treatment', 'organoid_age_days', 'publication', 'doi', 'batch', 'annot_level_1', 'annot_level_2', 'annot_level_3_rev2', 'annot_level_4_rev2', 'annot_region_rev2', 'annot_ntt_rev2', 'Hallmark_Glycolysis', 'hnoca_core', 'annot_level_2_extended', 'tissue_type', 'sex_ontology_term_id', 'donor_id', 'assay_ontology_term_id', 'self_reported_ethnicity_ontology_term_id', 'tissue_ontology_term_id', 'disease_ontology_term_id', 'development_stage_ontology_term_id', 'cell_type_ontology_term_id', 'is_primary_data', 'cell_type', 'assay', 'disease', 'sex', 'tissue', 'self_reported_ethnicity', 'development_stage', 'observation_joinid'
    var: 'gene_length', 'highly_variable', 'highly_variable_rank', 'highly_variable_nbatches', 'feature_is_filtered', 'feature_name', 'feature_reference', 'feature_biotype', 'feature_length', 'feature_type'
    uns: 'batch_condition', 'citation', 'default_embedding', 'organism', 'organism_ontology_term_id', 'schema_reference', 'schema_version', 'title', 'cell_type_colors', 'assay_differentiation_colors', 'tissue_colors'
    obsm: 'X_scpoli', 'X_umap_scpoli'
    layers: 'counts'
    obsp: 'knn_scpoli_connectivities', 'knn_scpoli_distances'
In [26]:
sc.pp.normalize_total(adata, target_sum=1e4, exclude_highly_expressed=True)
sc.pp.log1p(adata)
adata.layers['lognorm'] = adata.X.copy()
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5, batch_key='batch')

normalizing counts per cell The following highly-expressed genes are not considered during normalization factor computation:
['ENSG00000026025', 'ENSG00000054523', 'ENSG00000080824', 'ENSG00000081051', 'ENSG00000087086', 'ENSG00000091651', 'ENSG00000099797', 'ENSG00000101353', 'ENSG00000102003', 'ENSG00000105372', 'ENSG00000106261', 'ENSG00000111341', 'ENSG00000114735', 'ENSG00000118137', 'ENSG00000118271', 'ENSG00000119946', 'ENSG00000132002', 'ENSG00000133636', 'ENSG00000137135', 'ENSG00000141933', 'ENSG00000142002', 'ENSG00000147437', 'ENSG00000157005', 'ENSG00000157978', 'ENSG00000158874', 'ENSG00000160180', 'ENSG00000163220', 'ENSG00000163832', 'ENSG00000166426', 'ENSG00000167996', 'ENSG00000179965', 'ENSG00000183395', 'ENSG00000198712', 'ENSG00000198804', 'ENSG00000198899', 'ENSG00000198938', 'ENSG00000205542', 'ENSG00000251562', 'ENSG00000254505', 'ENSG00000278156']
    finished (0:00:02)
extracting highly variable genes
    finished (0:00:27)
--> added
    'highly_variable', boolean vector (adata.var)
    'means', float vector (adata.var)
    'dispersions', float vector (adata.var)
    'dispersions_norm', float vector (adata.var)
In [27]:
sc.tl.pca(adata, use_highly_variable=True)

computing PCA
    on highly variable genes
    with n_comps=50
    finished (0:00:26)
In [28]:
sc.pl.pca_variance_ratio(adata, log=True)
In [29]:
N_NB = int(0.5 * len(adata) ** 0.5)
if N_NB > 100:
    N_NB = 100
print(N_NB)

100
In [30]:
sc.pp.neighbors(adata, n_neighbors=N_NB, n_pcs=12, key_added="pca")

computing neighbors
    using 'X_pca' with n_pcs = 12

2026-03-11 13:10:06.891220: I tensorflow/core/util/port.cc:110] oneDNN custom operations are on. You may see slightly different numerical results due to floating-point round-off errors from different computation orders. To turn them off, set the environment variable `TF_ENABLE_ONEDNN_OPTS=0`.
2026-03-11 13:10:09.186013: I tensorflow/core/platform/cpu_feature_guard.cc:182] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.
To enable the following instructions: AVX2 AVX512F AVX512_VNNI FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.
2026-03-11 13:10:12.281137: W tensorflow/compiler/tf2tensorrt/utils/py_utils.cc:38] TF-TRT Warning: Could not find TensorRT

    finished: added to `.uns['pca']`
    `.obsp['pca_distances']`, distances for each pair of neighbors
    `.obsp['pca_connectivities']`, weighted adjacency matrix (0:06:06)
In [31]:
sc.tl.umap(adata, random_state=1, neighbors_key="pca")
# store coordinates in a named slot so to avoid confusion with batch-corrected
adata.obsm["X_umap_nocorr"] = adata.obsm["X_umap"].copy()
del adata.obsm["X_umap"]

computing UMAP
    finished: added
    'X_umap', UMAP coordinates (adata.obsm) (0:10:19)
In [ ]:
#sc.tl.draw_graph(adata, random_state=1, neighbors_key="pca")
#adata.obsm["X_draw_graph_fa_nocorr"] = adata.obsm["X_draw_graph_fa"].copy()
#del adata.obsm["X_draw_graph_fa"]
In [32]:
sc.pl.embedding(adata,  basis="X_umap_nocorr", color=['cell_type', 'organoid_age_days', 'assay_differentiation', 'tissue'], ncols=1)

/usr/local/lib/python3.8/dist-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
/usr/local/lib/python3.8/dist-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
/usr/local/lib/python3.8/dist-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
In [33]:
sc.external.pp.harmony_integrate(adata, "batch", random_state=5)

2026-03-11 13:26:33,666 - harmonypy - INFO - Computing initial centroids with sklearn.KMeans...
2026-03-11 13:27:45,801 - harmonypy - INFO - sklearn.KMeans initialization complete.
2026-03-11 13:27:47,263 - harmonypy - INFO - Iteration 1 of 10
2026-03-11 13:30:46,306 - harmonypy - INFO - Iteration 2 of 10
2026-03-11 13:33:47,004 - harmonypy - INFO - Iteration 3 of 10
2026-03-11 13:37:04,431 - harmonypy - INFO - Iteration 4 of 10
2026-03-11 13:40:12,136 - harmonypy - INFO - Iteration 5 of 10
2026-03-11 13:42:51,247 - harmonypy - INFO - Iteration 6 of 10
2026-03-11 13:44:27,083 - harmonypy - INFO - Iteration 7 of 10
2026-03-11 13:45:56,508 - harmonypy - INFO - Iteration 8 of 10
2026-03-11 13:47:20,135 - harmonypy - INFO - Converged after 8 iterations
In [34]:
sc.pp.neighbors(adata, n_neighbors=N_NB, n_pcs=12, use_rep='X_pca_harmony', key_added='harmony')

computing neighbors
    finished: added to `.uns['harmony']`
    `.obsp['harmony_distances']`, distances for each pair of neighbors
    `.obsp['harmony_connectivities']`, weighted adjacency matrix (0:05:44)
In [35]:
sc.tl.umap(adata, random_state=1, neighbors_key="harmony")
adata.obsm["X_umap_harmony"] = adata.obsm["X_umap"].copy()
del adata.obsm["X_umap"]

computing UMAP
    finished: added
    'X_umap', UMAP coordinates (adata.obsm) (0:11:48)
In [36]:
sc.pl.embedding(adata,  basis="X_umap_harmony", color=['cell_type', 'organoid_age_days', 'assay_differentiation', 'tissue'], ncols=1)

/usr/local/lib/python3.8/dist-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
/usr/local/lib/python3.8/dist-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
/usr/local/lib/python3.8/dist-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
In [37]:
sc.tl.draw_graph(adata, random_state=1, neighbors_key="harmony")
adata.obsm["X_draw_graph_fa_harmony"] = adata.obsm["X_draw_graph_fa"].copy()
del adata.obsm["X_draw_graph_fa"]

drawing single-cell graph using layout 'fa'
    finished: added
    'X_draw_graph_fa', graph_drawing coordinates (adata.obsm) (1:27:53)
In [38]:
sc.pl.embedding(adata,  basis="X_draw_graph_fa_harmony", color=['cell_type', 'organoid_age_days', 'assay_differentiation', 'tissue'], ncols=1)

/usr/local/lib/python3.8/dist-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
/usr/local/lib/python3.8/dist-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(
/usr/local/lib/python3.8/dist-packages/scanpy/plotting/_tools/scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  cax = scatter(

3.3 Visualization of receptors on UMAP¶

In [41]:
fn.CustomUmap(adata, genes, embedding="X_draw_graph_fa_harmony", var_col = "feature_name", gene_symbols="feature_name")

WARNING: saving figure to file figures/X_draw_graph_fa_harmonydraw_graph.png
In [42]:
fn.CustomUmap(adata, genes, embedding="X_umap_harmony", var_col = "feature_name", gene_symbols="feature_name")
In [43]:
available_genes = [gene for gene in genes if gene in adata.var['feature_name'].values]

if available_genes:
    sc.pl.dotplot(adata, available_genes, groupby='cell_type', gene_symbols='feature_name')
else:
    print("None of the specified genes are found in adata.var_names.")

/usr/local/lib/python3.8/dist-packages/scanpy/plotting/_dotplot.py:749: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored
  dot_ax.scatter(x, y, **kwds)

4. Savings¶

4.1 Timestamp finished computations¶

In [44]:
print(datetime.now())

2026-03-11 16:11:01.416770

4.2 Save adata¶

In [40]:
base_path = "../../../../DataDir/ExternalData/SingleCellData/"

adata.write(base_path +"HNOCA_adata_processed.h5ad")

adata_velasco.write(base_path + "HNOCA_velasco.h5ad")
adata_bhaduri.write(base_path + "HNOCA_bhaduri.h5ad")
adata_trujillo.write(base_path + "HNOCA_trujillo.h5ad")
adata_pellegrini.write(base_path + "HNOCA_pellegrini.h5ad")

4.3 Save python and html version of notebook¶

In [46]:
%%bash

# save also html and python versions for git
jupyter nbconvert ExplorationHNOCA.ipynb --to="python" --output="ExplorationHNOCA"
jupyter nbconvert ExplorationHNOCA.ipynb --to="html" --output="ExplorationHNOCA"

[NbConvertApp] Converting notebook ExplorationHNOCA.ipynb to python
[NbConvertApp] Writing 6087 bytes to ExplorationHNOCA.py
[NbConvertApp] Converting notebook ExplorationHNOCA.ipynb to html
[NbConvertApp] Writing 33634764 bytes to ExplorationHNOCA.html
In [ ]: