Exploration of hormonal receptor and metabolism genes¶

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 ipynbname
import os
import sys
sys.path.append('../0_HormonalGenes/exploration_scRNASeq/')
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-14 14:08:50.580066

2. Read input files¶

2.1 adata loading¶

In [4]:
adata = sc.read('../../DataDir/scRNASeq/cellxgene/e1bc953c-a4e9-4250-95a6-28b1ebe610d8.h5ad')
In [5]:
adata
Out[5]:
AnnData object with n_obs × n_vars = 163714 × 35474
    obs: 'cmo', 'condition', 'replicate', 'run', 'sample_id', 'days_of_differentiation', 'Singlet', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'total_counts_hb', 'log1p_total_counts_hb', 'pct_counts_hb', 'log_pct_counts_hb', 'log_total_counts', 'controls_substudy', 'phase', 'S_score', 'G2M_score', 'Consensus_nIntersection', 'Consensus_call', 'Leiden_1', 'donor_id', 'institute', 'library_id', 'manner_of_death', 'sample_source', 'sex_ontology_term_id', 'cell_type_ontology_term_id', 'sample_collection_method', 'tissue_type', 'sampled_site_condition', 'tissue_ontology_term_id', 'sample_preservation_method', 'suspension_type', 'cell_enrichment', 'assay_ontology_term_id', 'library_preparation_batch', 'library_sequencing_run', 'sequenced_fragment', 'is_primary_data', 'reference_genome', 'gene_annotation_version', 'alignment_software', 'disease_ontology_term_id', 'self_reported_ethnicity_ontology_term_id', 'development_stage_ontology_term_id', 'radial_tissue_term', 'cell_type', 'assay', 'disease', 'sex', 'tissue', 'self_reported_ethnicity', 'development_stage', 'observation_joinid'
    var: 'gene_symbols', 'feature_is_filtered', 'feature_name', 'feature_reference', 'feature_biotype', 'feature_length', 'feature_type'
    uns: 'ambient_count_correction', 'batch_condition', 'citation', 'doublet_detection', 'organism', 'organism_ontology_term_id', 'schema_reference', 'schema_version', 'study_pi', 'title'
    obsm: 'X_pca', 'X_pca_harmony', 'X_umap'
In [6]:
print('Loaded Normalizes AnnData object: number of cells', adata.n_obs)
print('Loaded Normalizes 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 Normalizes AnnData object: number of cells 163714
Loaded Normalizes AnnData object: number of genes 35474
Available metadata for each cell:  Index(['cmo', 'condition', 'replicate', 'run', 'sample_id',
       'days_of_differentiation', 'Singlet', 'n_genes_by_counts',
       'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts',
       'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes',
       'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes',
       'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt',
       'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo',
       'total_counts_hb', 'log1p_total_counts_hb', 'pct_counts_hb',
       'log_pct_counts_hb', 'log_total_counts', 'controls_substudy', 'phase',
       'S_score', 'G2M_score', 'Consensus_nIntersection', 'Consensus_call',
       'Leiden_1', 'donor_id', 'institute', 'library_id', 'manner_of_death',
       'sample_source', 'sex_ontology_term_id', 'cell_type_ontology_term_id',
       'sample_collection_method', 'tissue_type', 'sampled_site_condition',
       'tissue_ontology_term_id', 'sample_preservation_method',
       'suspension_type', 'cell_enrichment', 'assay_ontology_term_id',
       'library_preparation_batch', 'library_sequencing_run',
       'sequenced_fragment', 'is_primary_data', 'reference_genome',
       'gene_annotation_version', 'alignment_software',
       'disease_ontology_term_id', 'self_reported_ethnicity_ontology_term_id',
       'development_stage_ontology_term_id', 'radial_tissue_term', 'cell_type',
       'assay', 'disease', 'sex', 'tissue', 'self_reported_ethnicity',
       'development_stage', 'observation_joinid'],
      dtype='object')
In [7]:
adata_male = adata[adata.obs['sex'] == 'male'].copy()
adata_female = adata[adata.obs['sex'] == 'female'].copy()
In [8]:
adata.obs['sex'] = adata.obs['sex'].astype('category')

2.2 Receptors and metabolism signatures loading¶

Loading of hormonal receptor and metabolism gene signatures.

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_horm = sig["GeneName"].values.tolist()
In [12]:
signatures = '../../DataDir/bulkRNASeq/4.GeneSignatures/GeneSets/MetabolismGenes.txt'
In [13]:
sig = pd.read_csv(signatures, sep="\t", keep_default_na=False)  
print(sig.shape)
sig

(42, 3)
Out[13]:
GeneName Signature BiosyntheticPathway
0 CYP11A1 Cholesterol side-chain cleavage Steroid biogenesis
1 CYP17A1 17-alpha hydroxylase Steroid biogenesis
2 CYP11B1 11-beta-hydroxylase Steroid biogenesis
3 CYP11B2 18-hydroxylase Steroid biogenesis
4 CYP21A1 21-alpha hydroxylase Steroid biogenesis
5 HSD3B2 3-beta-hydroxysteroid dehydrogenase Steroid biogenesis
6 AKR1C4 3-alpha-HSD Steroid biogenesis
7 HSD11B1 11-beta-HSD Steroid biogenesis
8 HSD11B2 11-beta-HSD Steroid biogenesis
9 HSD17B3 17-beta HSD Steroid biogenesis
10 HSD3B1 3-beta-hydroxysteroid dehydrogenase Steroid biogenesis
11 AKR1C1 20-alpha-HSD Steroid biogenesis
12 AKR1C2 20-alpha-HSD Steroid biogenesis
13 AKR1C3 20-alpha-HSD Steroid biogenesis
14 AKR1C4 20-beta-HSD Steroid biogenesis
15 SRD5A1 5-alpha- reductase Steroid biogenesis
16 SRD5A2 5-alpha- reductase Steroid biogenesis
17 STAR Steroidogenic acute regulatory protein Steroid biogenesis
18 SRD5A3 5-alpha- reductase Steroid biogenesis
19 CYP19A1 aromatase Steroid biogenesis
20 ACAT1 Thiolase Cholesterol synthesis
21 ACAT2 Thiolase Cholesterol synthesis
22 HMGCS1 HMG-CoA synthase (soluble) Cholesterol synthesis
23 HMGCS2 HMG-CoA synthase (mitochondrial) Cholesterol synthesis
24 HMGCR HMG CoA reductase Cholesterol synthesis
25 MVK Mevalonate kinase Cholesterol synthesis
26 PMVK phosphomevalonate kinase Cholesterol synthesis
27 MVD Mevalonate 5-pyrophosphate decarboxylase Cholesterol synthesis
28 IDI1 isopentenyl-PP-isomerase Cholesterol synthesis
29 IDI2 isopentenyl-PP-isomerase Cholesterol synthesis
30 FDPS Farnesyl-PP-synthase Cholesterol synthesis
31 FDFT1 squalene synthase Cholesterol synthesis
32 SQLE squalene monoxygenase Cholesterol synthesis
33 SQLE squalene epoxydase Cholesterol synthesis
34 GGPS1 geranylgeranyl pyrophosphate synthase Cholesterol synthesis
35 LSS Lanosterol synthase Cholesterol synthesis
36 DHCR24 24-dehydrocholesterol reductase Cholesterol synthesis
37 DHCR7 7-dehydrocholesterol reductase Cholesterol synthesis
38 HSD17B7 hydroxysteroid dehydrogenase 7 Cholesterol synthesis
39 MSMO1 methylsterol monooxygenase 1 Cholesterol synthesis
40 NSDHL NAD-dependent steroid dehydrogenase-like Cholesterol synthesis
41 SC5D sterol-C5-desaturase Cholesterol synthesis
In [14]:
genes_met = sig["GeneName"].values.tolist()

3. Visualizations¶

3.1 Counts from adata¶

In [15]:
adata.obsm
Out[15]:
AxisArrays with keys: X_pca, X_pca_harmony, X_umap
In [16]:
sc.pl.embedding(adata, basis="X_umap", color=['n_genes_by_counts',"total_counts", 'pct_counts_mt', 'pct_counts_ribo'])

3.2 Clusters annotation¶

In [17]:
adata.obs['condition'].value_counts()
Out[17]:
condition
ARYL_HYD_ANTAGONIST     17206
ANDROGEN_ANTAGONIST     14460
LIVER-X_ANTAGONIST      12290
ESTROGEN_AGONIST        12242
RET_AGONIST             12205
THYROID_ANTAGONIST      11168
LIVER-X_AGONIST         10934
ARYL_HYD_AGONIST        10763
ESTROGEN_ANTAGONIST      9934
GLUCOCORT_ANTAGONIST     9472
DMSO                     9334
ANDROGEN_AGONIST         9309
GLUCOCORT_AGONIST        8703
THYROID_AGONIST          8373
RET_ANTAGONIST           7321
Name: count, dtype: int64
In [18]:
# Map original → your desired labels
condition_map = {
    "DMSO": "DMSO",

    "ANDROGEN_AGONIST": "AND_AGONIST",
    "ANDROGEN_ANTAGONIST": "AND_INHIBITOR",

    "ESTROGEN_AGONIST": "EST_AGONIST",
    "ESTROGEN_ANTAGONIST": "EST_INHIBITOR",

    "GLUCOCORT_AGONIST": "GC_AGONIST",
    "GLUCOCORT_ANTAGONIST": "GC_INHIBITOR",

    "THYROID_AGONIST": "THY_AGONIST",
    "THYROID_ANTAGONIST": "THY_INHIBITOR",

    "RET_AGONIST": "RA_AGONIST",
    "RET_ANTAGONIST": "RA_INHIBITOR",

    "LIVER-X_AGONIST": "LX_AGONIST",
    "LIVER-X_ANTAGONIST": "LX_INHIBITOR",

    "ARYL_HYD_AGONIST": "AH_AGONIST",
    "ARYL_HYD_ANTAGONIST": "AH_INHIBITOR",
}

adata.obs["exposure"] = adata.obs["condition"].map(condition_map)

# Order (optional but recommended)
cat_order = [
    "DMSO",
    "AH_AGONIST", "AH_INHIBITOR",
    "AND_AGONIST", "AND_INHIBITOR",
    "EST_AGONIST", "EST_INHIBITOR",
    "GC_AGONIST", "GC_INHIBITOR",
    "LX_AGONIST", "LX_INHIBITOR",
    "RA_AGONIST", "RA_INHIBITOR",
    "THY_AGONIST", "THY_INHIBITOR",
]

adata.obs["exposure"] = pd.Categorical(
    adata.obs["exposure"],
    categories=cat_order,
    ordered=True
)

# SAME colors, just reassigned to your labels
condition_colors = {
    'DMSO': '#4D4D4D',

    'AH_AGONIST': '#F8766D',
    'AH_INHIBITOR': '#F8766D50',

    'AND_AGONIST': '#fccb17',
    'AND_INHIBITOR': '#C49A0050',

    'EST_AGONIST': '#53B400',
    'EST_INHIBITOR': '#53B40050',

    'GC_AGONIST': '#00C094',
    'GC_INHIBITOR': '#00C09450',

    'LX_AGONIST': '#00B6EB',
    'LX_INHIBITOR': '#00B6EB50',

    'RA_AGONIST': '#A58AFF',
    'RA_INHIBITOR': '#A58AFF50',

    'THY_AGONIST': '#FB61D7',
    'THY_INHIBITOR': '#FB61D750'
}
In [19]:
outdir = "../../../ENDpoiNTs/FigPaper/"

# create folder if it doesn't exist
os.makedirs(outdir, exist_ok=True)

# tell Scanpy where to save
sc.settings.figdir = outdir

# Plot
sc.pl.embedding(
    adata,
    basis="X_umap",
    color="exposure",
    palette=condition_colors,
    ncols=1,
    size=5,
    save="umap_exposure.png"
)

WARNING: saving figure to file ../../../ENDpoiNTs/FigPaper/X_umapumap_exposure.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(
In [20]:
# Plot
sc.pl.embedding(
    adata,
    basis="X_umap",
    color="cell_type",
    #palette=condition_colors,
    ncols=1,
    size=5,
    save="_celltype.png"
)

WARNING: saving figure to file ../../../ENDpoiNTs/FigPaper/X_umap_celltype.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(
In [21]:
sc.pl.embedding(adata,  basis="X_umap", color=['condition', 'donor_id', 'cell_type'], 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 genes¶

In [22]:
adata.var.head()
Out[22]:
gene_symbols feature_is_filtered feature_name feature_reference feature_biotype feature_length feature_type
ENSG00000243485 MIR1302-2HG False MIR1302-2HG NCBITaxon:9606 gene 517 lncRNA
ENSG00000237613 FAM138A False FAM138A NCBITaxon:9606 gene 1015 lncRNA
ENSG00000186092 OR4F5 False OR4F5 NCBITaxon:9606 gene 2618 protein_coding
ENSG00000239945 AL627309.3 False ENSG00000239945 NCBITaxon:9606 gene 1319 lncRNA
ENSG00000239906 AL627309.2 False ENSG00000239906 NCBITaxon:9606 gene 323 lncRNA
In [23]:
adata.var['gene_symbols'].duplicated().sum()
Out[23]:
0
In [24]:
adata.var['gene_symbols'].isna().sum()
Out[24]:
0
In [25]:
fn.CustomUmap(adata, genes_horm, embedding="X_umap", var_col = "gene_symbols", gene_symbols="gene_symbols")
In [26]:
available_genes = [gene for gene in genes_horm if gene in adata.var['gene_symbols'].values]

if available_genes:
    sc.pl.dotplot(adata, available_genes, groupby='cell_type', gene_symbols='gene_symbols')
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)
In [27]:
available_genes = [gene for gene in genes_horm if gene in adata.var['gene_symbols'].values]

if available_genes:
    sc.pl.dotplot(adata, available_genes, groupby=['sex', 'cell_type'], gene_symbols='gene_symbols')
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)
In [28]:
sc.pl.matrixplot(
    adata,
    var_names=available_genes,
    groupby=['cell_type', 'sex'],
    gene_symbols='gene_symbols',
    cmap='viridis'
)
In [29]:
for ct in adata.obs['cell_type'].cat.categories:
    adata_ct = adata[adata.obs['cell_type'] == ct]
    sc.pl.dotplot(
        adata_ct,
        available_genes,
        groupby='sex',
        gene_symbols='gene_symbols',
        title=ct
    )

/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)

/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)

MALE LINE: CTL08¶

In [30]:
available_genes = [gene for gene in genes_horm if gene in adata_male.var['gene_symbols'].values]

if available_genes:
    sc.pl.dotplot(adata_male, available_genes, groupby='cell_type', gene_symbols='gene_symbols')
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)

FEMALE LINE: CTL04¶

In [31]:
available_genes = [gene for gene in genes_horm if gene in adata_female.var['gene_symbols'].values]

if available_genes:
    sc.pl.dotplot(adata_female, available_genes, groupby='cell_type', gene_symbols='gene_symbols')
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)

3.4 Visualization of metabolism genes¶

In [42]:
fn.CustomUmap(adata, genes_met, embedding="X_umap", var_col = "gene_symbols", gene_symbols="gene_symbols")

Missing: {'CYP21A1'}
In [43]:
available_genes = [gene for gene in genes_met if gene in adata.var['gene_symbols'].values]

if available_genes:
    sc.pl.dotplot(adata, available_genes, groupby='cell_type', gene_symbols='gene_symbols')
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)
In [44]:
available_genes = [gene for gene in genes_met if gene in adata.var['gene_symbols'].values]

if available_genes:
    sc.pl.dotplot(adata, available_genes, groupby=['sex', 'cell_type'], gene_symbols='gene_symbols')
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)
In [45]:
sc.pl.matrixplot(
    adata,
    var_names=available_genes,
    groupby=['cell_type', 'sex'],
    gene_symbols='gene_symbols',
    cmap='viridis'
)
In [46]:
for ct in adata.obs['cell_type'].cat.categories:
    adata_ct = adata[adata.obs['cell_type'] == ct]
    sc.pl.dotplot(
        adata_ct,
        available_genes,
        groupby='sex',
        gene_symbols='gene_symbols',
        title=ct
    )

/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)

/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)

MALE LINE: CTL08¶

In [47]:
available_genes = [gene for gene in genes_met if gene in adata_male.var['gene_symbols'].values]

if available_genes:
    sc.pl.dotplot(adata_male, available_genes, groupby='cell_type', gene_symbols='gene_symbols')
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)

FEMALE LINE: CTL04¶

In [48]:
available_genes = [gene for gene in genes_met if gene in adata_female.var['gene_symbols'].values]

if available_genes:
    sc.pl.dotplot(adata_female, available_genes, groupby='cell_type', gene_symbols='gene_symbols')
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 [49]:
print(datetime.now())

2026-04-14 14:44:47.376742

4.2 Save python and html version of notebook¶

In [50]:
nb_fname = ipynbname.name()
nb_fname
Out[50]:
'7.GeneSignatures_exploration'
In [51]:
%%bash -s "$nb_fname"
jupyter nbconvert "$1".ipynb --to="python"
jupyter nbconvert "$1".ipynb --to="html"

[NbConvertApp] Converting notebook 7.GeneSignatures_exploration.ipynb to python
[NbConvertApp] Writing 8943 bytes to 7.GeneSignatures_exploration.py
[NbConvertApp] Converting notebook 7.GeneSignatures_exploration.ipynb to html
[NbConvertApp] Writing 17920953 bytes to 7.GeneSignatures_exploration.html
In [ ]: