Exploration of hormonal receptor genes in Castaldi et al organoid dataset¶

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-04-27 13:51:09.155479

2. Read input files¶

2.1 adata loading¶

Let's start from the raw data. Not the processed adata with a selection of genes as in Castaldi_FA.ipynb

In [4]:
path = '../../../../Castaldi_multiplexingCBO/'
input_file_raw = path + 'adataPagaRaw.h5ad'
In [5]:
adata = sc.read(input_file_raw)
In [6]:
adata
Out[6]:
AnnData object with n_obs × n_vars = 14913 × 33538
    obs: 'dataset', 'cellID', 'cellID_newName', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'stage', 'type', 'id_stage', 'cellID_newName_type', 'S_score', 'G2M_score', 'phase', 'leidenAnnotated', 'leiden_1.2', 'endpoint_GlutamatergicNeurons_late', 'endpoint_GlutamatergicNeurons_early', 'endpoint_MigratingNeurons', 'endpoint_OuterRadialGliaAstrocytes', 'endpoint_Interneurons', 'endpoint_Interneurons_GAD2', 'endpoint_CajalR_like', 'Exc_Lineage', 'endpoint_GlutamatergicNeurons_both'
    var: 'highly_variable'
    uns: 'cellID_colors', 'cellID_newName_colors', 'cellID_newName_type_colors', 'dataset_colors', 'stage_colors', 'type_colors'
In [7]:
print('Loaded AnnData object: number of cells', adata.n_obs)
print('Loaded AnnData object: number of genes', adata.n_vars)

# To see the columns of the metadata (information available for each cell)  
print('Available metadata for each cell: ', adata.obs.columns)

Loaded AnnData object: number of cells 14913
Loaded AnnData object: number of genes 33538
Available metadata for each cell:  Index(['dataset', 'cellID', 'cellID_newName', 'n_genes_by_counts',
       'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts',
       'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt',
       'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo',
       'stage', 'type', 'id_stage', 'cellID_newName_type', 'S_score',
       'G2M_score', 'phase', 'leidenAnnotated', 'leiden_1.2',
       'endpoint_GlutamatergicNeurons_late',
       'endpoint_GlutamatergicNeurons_early', 'endpoint_MigratingNeurons',
       'endpoint_OuterRadialGliaAstrocytes', 'endpoint_Interneurons',
       'endpoint_Interneurons_GAD2', 'endpoint_CajalR_like', 'Exc_Lineage',
       'endpoint_GlutamatergicNeurons_both'],
      dtype='object')
In [8]:
np.unique(adata.obs.values[:,14])
Out[8]:
array(['downstream', 'upstream'], dtype=object)

2.2 Receptors signature loading¶

Loading of hormonal receptor gene signature.

In [9]:
signatures = '../../../../DataDir/ExternalData/Receptors/EndocrineKeyGenes.txt'
In [10]:
sig = pd.read_csv(signatures, sep="\t", keep_default_na=False)  
print(sig.shape)
sig

(30, 2)
Out[10]:
GeneName Signature
0 THRB THY
1 THRA THY
2 THRAP3 THY
3 DIO1 THY
4 DIO2 THY
5 DIO3 THY
6 SLC16A10 THY
7 SLC16A2 THY
8 SLC7A5 THY
9 KLF9 THY
10 THRSP THY
11 ESRRG EST
12 ESRRA EST
13 GPER1 EST
14 ESR1 EST
15 ESR2 EST
16 ESRRB EST
17 CYP19A1 EST
18 AR AND
19 RBP4 RA
20 RARA RA
21 RARB RA
22 RARG RA
23 RXRA RA
24 RXRB RA
25 RXRG RA
26 AHR AH
27 NR3C1 GC
28 NR1H2 LX
29 NR1H3 LX
In [11]:
genes = sig["GeneName"].values.tolist()

3. Visualizations¶

3.1 Counts from adata¶

In [12]:
adata.obsm
Out[12]:
AxisArrays with keys:
In [13]:
adata.X.data[:10]   # se sparse
Out[13]:
array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.], dtype=float32)
In [14]:
adata.layers['counts'] = adata.X.copy()
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)

normalizing counts per cell The following highly-expressed genes are not considered during normalization factor computation:
['MTRNR2L12', 'HBB', 'FTH1', 'MALAT1', 'HSP90AA1', 'MT-CO1', 'MT-CO2', 'MT-ATP6', 'MT-CO3', 'MT-ND4', 'MT-CYB']
    finished (0:00:00)
extracting highly variable genes
    finished (0:00:00)
--> added
    'highly_variable', boolean vector (adata.var)
    'means', float vector (adata.var)
    'dispersions', float vector (adata.var)
    'dispersions_norm', float vector (adata.var)
In [15]:
sc.tl.pca(adata, use_highly_variable=True)
sc.pl.pca_variance_ratio(adata, log=True)

computing PCA
    on highly variable genes
    with n_comps=50
    finished (0:00:03)
In [16]:
N_NB = int(0.5 * len(adata) ** 0.5)
if N_NB > 100:
    N_NB = 100
print(N_NB) 
sc.pp.neighbors(adata, n_neighbors=N_NB, n_pcs=12, key_added="pca")

61
computing neighbors
    using 'X_pca' with n_pcs = 12

2026-04-27 13:51:21.944483: 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-04-27 13:51:21.984035: 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-04-27 13:51:23.814156: 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:00:33)
In [17]:
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:00:12)

3.2 Clusters annotation¶

In [18]:
sc.pl.embedding(adata,  basis="X_umap_nocorr", color=['leidenAnnotated', 'dataset'], 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(
In [19]:
sc.external.pp.harmony_integrate(adata, "dataset", random_state=5)

2026-04-27 13:52:03,598 - harmonypy - INFO - Computing initial centroids with sklearn.KMeans...
2026-04-27 13:52:07,314 - harmonypy - INFO - sklearn.KMeans initialization complete.
2026-04-27 13:52:07,365 - harmonypy - INFO - Iteration 1 of 10
2026-04-27 13:52:09,963 - harmonypy - INFO - Iteration 2 of 10
2026-04-27 13:52:12,567 - harmonypy - INFO - Converged after 2 iterations
In [20]:
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:00:07)
In [21]:
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:00:12)
In [22]:
sc.pl.embedding(adata,  basis="X_umap_harmony", color=['leidenAnnotated', 'dataset'], 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(
In [23]:
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) (0:01:45)
In [31]:
sc.settings.figdir = "../../../../FigPaper/"
 
sc.set_figure_params(dpi_save=600)
 
sc.pl.embedding(
    adata,
    basis="X_draw_graph_fa_harmony",
    color=['leidenAnnotated'],
    ncols=1,
    save="CastaldiAllV2_highres.png"
)

WARNING: saving figure to file ../../../../FigPaper/X_draw_graph_fa_harmonyCastaldiAllV2_highres.png

/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 [32]:
fn.CustomUmap(adata, genes, embedding="X_umap_harmony")
In [33]:
fn.CustomUmap(adata, genes, embedding="X_draw_graph_fa_harmony")
In [34]:
available_genes = [gene for gene in genes if gene in adata.var_names]

if available_genes:
    sc.pl.dotplot(adata, available_genes, groupby='leidenAnnotated')
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. Save Notebooks¶

Adata is not saved, since no new computations have been performed. I just save the notebooks.

4.1 Timestamp finished computations¶

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

2026-04-27 14:40:03.639114

4.2 Save adata¶

In [36]:
adata.write('../../../../DataDir/ExternalData/SingleCellData/Castaldi_adata.h5ad')
In [ ]: