Exploration of hormonal receptor genes in Castaldi et al organoid dataset - Metacell calculation by SeaCells (only day 50)¶

Reference paper

QC filter on cells and cell annotation already performed by authors of the paper.

Here in the first part:

  • subsetting adata with raw counts with only day 50 organoids
  • Normalization and log10 transformation by Scanpy functions
  • Feature selection (HVG) by Scanpy functions
  • Dimensionality reduction
  • Exploration of hormonal receptor gene expressions in single cells

and in the second part:

  • Metacell calculation by SEACells (following this tutorial)
  • Exploration of hormonal receptor gene expressions in metacells

1. Environment¶

1.1 Modules¶

In [1]:
import os
import sys

import numpy as np
import pandas as pd
import scanpy as sc

import pickle

#Plotting
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns

#utils
import ipynbname
from datetime import datetime

# SeaCell
import SEACells

#import custom functions
sys.path.append('../../')
import functions as fn

findfont: Font family ['Raleway'] not found. Falling back to DejaVu Sans.
findfont: Font family ['Lato'] not found. Falling back to DejaVu Sans.
In [2]:
# Some plotting aesthetics
%matplotlib inline

sns.set_style('ticks')
matplotlib.rcParams['figure.figsize'] = [3.5, 3.5]
matplotlib.rcParams['figure.dpi'] = 100
In [3]:
print("Scanpy version: ", sc.__version__)
print("Pandas version: ", pd.__version__)
print("SEACell version: ", SEACells.__version__)

Scanpy version:  1.9.3
Pandas version:  2.0.0
SEACell version:  0.3.3

1.2 Settings¶

In [4]:
sc.settings.verbosity = 3
sc.settings.set_figure_params(dpi=80)

1.3 Files and parameters¶

In [5]:
path = '../../../../../Castaldi_multiplexingCBO/'
input_file_raw = path + 'adataPagaRaw.h5ad'
input_file_processed = path + 'adataPaga.h5ad'

output_file = '../../../../../DataDir/ExternalData/SingleCellData/CastaldiAdata_Metacells_day50.h5ad'

1.4 Start computations¶

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

2026-01-08 13:33:20.218432

2. Data Load¶

2.1 Read adata file¶

  • SEACells uses as input filtered unnormalized counts from scRNASeq data, that are subjected to normalization and log-transformation before metacell assignment
  • The sum of expression levels from cells to metacells will occur at the count level
In [6]:
adata_raw = sc.read(input_file_raw)
In [7]:
adata_raw
Out[7]:
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 [8]:
print(adata_raw.X)

  (0, 21)	1.0
  (0, 32)	1.0
  (0, 51)	1.0
  (0, 55)	1.0
  (0, 66)	1.0
  (0, 74)	1.0
  (0, 77)	1.0
  (0, 78)	1.0
  (0, 86)	1.0
  (0, 91)	1.0
  (0, 132)	1.0
  (0, 153)	1.0
  (0, 154)	8.0
  (0, 161)	3.0
  (0, 178)	3.0
  (0, 190)	3.0
  (0, 198)	1.0
  (0, 201)	3.0
  (0, 220)	1.0
  (0, 226)	1.0
  (0, 228)	1.0
  (0, 229)	1.0
  (0, 240)	1.0
  (0, 244)	1.0
  (0, 257)	1.0
  :	:
  (14912, 33250)	1.0
  (14912, 33252)	1.0
  (14912, 33254)	3.0
  (14912, 33257)	1.0
  (14912, 33286)	1.0
  (14912, 33297)	2.0
  (14912, 33326)	6.0
  (14912, 33327)	1.0
  (14912, 33376)	2.0
  (14912, 33400)	1.0
  (14912, 33443)	1.0
  (14912, 33445)	2.0
  (14912, 33446)	3.0
  (14912, 33474)	1.0
  (14912, 33493)	1.0
  (14912, 33496)	9.0
  (14912, 33497)	17.0
  (14912, 33498)	35.0
  (14912, 33499)	26.0
  (14912, 33501)	22.0
  (14912, 33502)	23.0
  (14912, 33503)	7.0
  (14912, 33505)	19.0
  (14912, 33506)	5.0
  (14912, 33508)	8.0
In [9]:
adata_raw.obs
Out[9]:
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 ... 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
AAACCTGAGAGACTAT-1_DownD250 DownD250 MIFF1 CTL02A 2586 7.858254 6358.0 8.757627 175.0 5.170484 2.752438 ... 0 1 0 0 0 0 0 0 GlutamatergicNeurons_late 1
AAACCTGCATGGTTGT-1_DownD250 DownD250 KOLF CTL08A 1372 7.224753 2540.0 7.840313 73.0 4.304065 2.874016 ... 0 1 0 0 0 0 0 0 GlutamatergicNeurons_late 1
AAACCTGTCAGTTAGC-1_DownD250 DownD250 3391B CTL01 1216 7.104144 2859.0 7.958577 43.0 3.784190 1.504022 ... 6 0 0 0 0 1 0 0 Other 0
AAACGGGCAGCTATTG-1_DownD250 DownD250 MIFF1 CTL02A 1682 7.428333 3906.0 8.270525 117.0 4.770685 2.995392 ... 6 0 0 0 0 1 0 0 Other 0
AAACGGGGTGCACGAA-1_DownD250 DownD250 MIFF1 CTL02A 930 6.836259 1733.0 7.458186 166.0 5.117994 9.578765 ... 11 1 1 1 0 1 1 1 Shared 1
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
TTTGTCAGTGACGGTA-1_UpD50 UpD50 809 CTL04E 1644 7.405496 3618.0 8.193954 36.0 3.610918 0.995025 ... 11 1 1 1 0 1 1 1 Shared 1
TTTGTCAGTTGGACCC-1_UpD50 UpD50 809 CTL04E 1278 7.153834 2861.0 7.959276 30.0 3.433987 1.048584 ... 8 0 1 1 0 0 1 0 glutamatergic_neurons_early 1
TTTGTCATCGCCTGTT-1_UpD50 UpD50 3391B CTL01 1361 7.216709 2898.0 7.972121 20.0 3.044523 0.690131 ... 12 0 0 0 0 0 0 1 Other 0
TTTGTCATCGCGTAGC-1_UpD50 UpD50 3391B CTL01 1203 7.093405 2296.0 7.739359 11.0 2.484907 0.479094 ... 10 1 1 1 1 1 1 1 Shared 1
TTTGTCATCGTTTATC-1_UpD50 UpD50 809 CTL04E 2815 7.943073 8304.0 9.024613 171.0 5.147494 2.059248 ... 2 1 1 1 1 1 1 1 Shared 1

14913 rows × 31 columns

2.2 Subset and process adata¶

Subsetting adata to have only 'DownD50' and 'UpD50'

In [10]:
adata = adata_raw[adata_raw.obs['dataset'].isin(['DownD50', 'UpD50'])].copy()
In [11]:
adata
Out[11]:
AnnData object with n_obs × n_vars = 4987 × 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 [13]:
adata.layers['counts'] = adata.X.copy()
  1. Normalization and log10 transformation by Scanpy functions
  2. Feature selection (HVG) by Scanpy functions
  3. Dimensionality reduction (PCA and then UMAP)
In [14]:
sc.pp.normalize_total(adata)
sc.pp.log1p(adata)
adata.layers['lognorm'] = adata.X.copy()

sc.pp.highly_variable_genes(adata, inplace=True)
sc.tl.pca(adata, use_highly_variable=True)
sc.pp.neighbors(adata, use_rep='X_pca')
sc.tl.umap(adata)
sc.tl.draw_graph(adata, layout='fa')

normalizing counts per cell
    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)
computing PCA
    on highly variable genes
    with n_comps=50
    finished (0:00:00)
computing neighbors

2026-01-08 13:33:27.777781: 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-01-08 13:33:27.819303: 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-01-08 13:33:29.230686: W tensorflow/compiler/tf2tensorrt/utils/py_utils.cc:38] TF-TRT Warning: Could not find TensorRT

    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:10)
computing UMAP
    finished: added
    'X_umap', UMAP coordinates (adata.obsm) (0:00:06)
drawing single-cell graph using layout 'fa'
    finished: added
    'X_draw_graph_fa', graph_drawing coordinates (adata.obsm) (0:00:28)
In [15]:
sc.pl.embedding(adata,  basis="X_umap", color=['leidenAnnotated'], 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(
In [16]:
sc.pl.embedding(adata,  basis="X_draw_graph_fa", color=['leidenAnnotated'], 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(

2.3 Explore hormonal gene expressions¶

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

(30, 2)
In [18]:
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)
In [19]:
fn.CustomUmap(adata, genes, embedding="X_umap")

The following marker genes are missing:  set()
In [20]:
fn.CustomUmap(adata, genes, embedding="X_draw_graph_fa")

The following marker genes are missing:  set()

3. Run SEACells¶

We follow the workflow described in the SEACells tutorial

3.1 Define parameters¶

Parameters to be defined:

  • Number of metacells (graining level): choosing one metacell for every 60 single-cells.
  • X_pca as the key in .obsm for computing in metacells (to be changed to X_svd in case od ATAC data)
  • Number of eigenvalues for metacell initialization

We then employ the SEACells.core.SEACells function to initialize the model.

In [21]:
## Core parameters 
n_SEACells = round(adata.n_obs/60)
print(n_SEACells)

build_kernel_on = 'X_pca' # key in ad.obsm to use for computing metacells
                          # This would be replaced by 'X_svd' for ATAC data

## Additional parameters
n_waypoint_eigs = 10 # Number of eigenvalues to consider when initializing metacells

83
In [22]:
#help(SEACells.core.SEACells)

model = SEACells.core.SEACells(adata, 
                  build_kernel_on=build_kernel_on, 
                  n_SEACells=n_SEACells, 
                  n_waypoint_eigs=n_waypoint_eigs,
                  convergence_epsilon = 1e-5)

Welcome to SEACells!

3.2 Construct kernel matrix¶

construct_kernel_matrix function constructs the kernel matrix from data matrix using PCA/SVD and nearest neighbors. Key parameters are:

  • n_neighbors: (int) number of nearest neighbors to use for graph construction (default 15). Increasing this number results in a more homogeneous distribution of metacell size, but may impact on the preservation of rare cell types. Usually stable performance in the range 5-30.
  • graph_construction: method for graph construction.
In [23]:
model.construct_kernel_matrix(n_neighbors=20)
M = model.kernel_matrix

Computing kNN graph using scanpy NN ...
computing neighbors
    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:00)
Computing radius for adaptive bandwidth kernel...

  0%|          | 0/4987 [00:00<?, ?it/s]
Making graph symmetric...
Parameter graph_construction = union being used to build KNN graph...
Computing RBF kernel...

  0%|          | 0/4987 [00:00<?, ?it/s]
Building similarity LIL matrix...

  0%|          | 0/4987 [00:00<?, ?it/s]
Constructing CSR matrix...

3.3 Initialize archetypes¶

  • initialize_archetypes function initializes the B matrix which defines cells as SEACells.
  • The function uses waypoint analysis for initialization into to fully cover the phenotype space, and then greedily selects the remaining cells.
In [24]:
# Initialize archetypes
model.initialize_archetypes()

Building kernel on X_pca
Computing diffusion components from X_pca for waypoint initialization ... 
Determing nearest neighbor graph...
computing neighbors
    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:00)
Done.
Sampling waypoints ...
Done.
Selecting 78 cells from waypoint initialization.
Initializing residual matrix using greedy column selection
Initializing f and g...

100%|██████████| 15/15 [00:00<00:00, 109.11it/s]
Selecting 5 cells from greedy initialization.
























In [25]:



     


# Plot the initilization to ensure they are spread across phenotypic space
SEACells.plot.plot_initialization(adata, model)



     




















    
    



































3.4 Fit Model¶
We apply the model.fit function.



















In [26]:



     


model.fit(min_iter=10, max_iter=100)
model



     




















    
    






Randomly initialized A matrix.
Setting convergence threshold at 0.00137
Starting iteration 1.
Completed iteration 1.
Starting iteration 10.
Completed iteration 10.
Starting iteration 20.
Completed iteration 20.
Starting iteration 30.
Completed iteration 30.
Converged after 30 iterations.








    
    
Out[26]:








<SEACells.core.SEACells at 0x155457404700>




























4. Access SEACells results¶






















4.1 Check model convergence¶


















In [27]:



     


# Plot model convergence
model.plot_convergence()



     




















    
    



































4.2 Inspect SEACells soft assignment¶
    

  • SEACells analysis returns soft assignments of cells to SEACells (full assignment matrix can be accessed at model.A_)
    • 

  • The majority of single-cells are assigned to no more than 4 archetypes with non-trivial weight
    • 

  • The algorithm returns the top 5 metacell assignments as well as the corresponding assignment weights in the function model.get_soft_assignments()
    • 




















In [28]:



     


#model.A_

plt.figure(figsize=(4,3))
sns.distplot((model.A_.T > 0.1).sum(axis=1), kde=False)
plt.title(f'Non-trivial (> 0.1) assignments per cell')
plt.xlabel('# Non-trivial SEACell Assignments')
plt.ylabel('# Cells')
plt.show()

plt.figure(figsize=(4,3))
b = np.partition(model.A_.T, -5)    
sns.heatmap(np.sort(b[:,-5:])[:, ::-1], cmap='viridis', vmin=0)
plt.title('Strength of top 5 strongest assignments')
plt.xlabel('$n^{th}$ strongest assignment')
plt.show()



     




















    
    






/tmp/ipykernel_2622269/2880491417.py:4: UserWarning: 

`distplot` is a deprecated function and will be removed in seaborn v0.14.0.

Please adapt your code to use either `displot` (a figure-level function with
similar flexibility) or `histplot` (an axes-level function for histograms).

For a guide to updating your code to use the new functions, please see
https://gist.github.com/mwaskom/de44147ed2974457ad6372750bbe5751

  sns.distplot((model.A_.T > 0.1).sum(axis=1), kde=False)








    
    















    
    































In [29]:



     


labels,weights = model.get_soft_assignments()
#labels.head()



     























4.3 Inspect and retrieve SEACells hard assignment¶
We will use hard assignment of each cell to a single metacell to define our metacells for downstream steps. Hard assingment can be retreived in two ways:


    

  • in the modified anndata object in .obs['SEACell']
    • 

  • from the model using .get_hard_assignments()
    • 




















In [30]:



     


adata.obs[['SEACell']].head()

# Alternatively: 
# model.get_hard_assignments().head()



     




















    
    
Out[30]:











  
    

      
      SEACell
    

    

      index
      
    

  
  
    

      AAACCTGCACGAGAGT-1_DownD50
      SEACell-27
    

    

      AAACCTGCAGCGTCCA-1_DownD50
      SEACell-42
    

    

      AAACGGGAGCTAACAA-1_DownD50
      SEACell-39
    

    

      AAACGGGCAGCTGGCT-1_DownD50
      SEACell-6
    

    

      AAACGGGGTAACGCGA-1_DownD50
      SEACell-77
    

  































5. Summarizing data by hard assignment¶
    

  • The original dataset can be summarized by SEACell by aggregating cells within each SEACell, summing over all raw data for all cells belonging to a SEACell.
    • 

  • The output of this function is an anndata object of shape n_metacells x original_data_dimension.
    • 

  • Data is unnormalized and raw aggregated counts are stored in X.
    • 

  • Attributes associated with variables (.var) are copied over, but relevant per SEACell attributes must be manually copied, since certain attributes may need to be summed, or averaged etc, depending on the attribute.
    • 




















In [31]:



     


# By default, ad.raw is used for summarization. Other layers present in the anndata can be specified using the parameter summarize_layer parameter
SEA_adata = SEACells.core.summarize_by_SEACell(adata, SEACells_label='SEACell', summarize_layer='counts', celltype_label='leidenAnnotated')
SEA_adata



     




















    
    






100%|██████████| 83/83 [00:00<00:00, 187.54it/s]








    
    
Out[31]:








AnnData object with n_obs × n_vars = 83 × 33538
    obs: 'leidenAnnotated', 'leidenAnnotated_purity'
    layers: 'raw'
























In [32]:



     


SEA_adata.obs.head()



     




















    
    
Out[32]:











  
    

      
      leidenAnnotated
      leidenAnnotated_purity
    

  
  
    

      SEACell-27
      RadialGliaProgenitors
      0.953488
    

    

      SEACell-42
      GlutamatergicNeurons_early
      1.000000
    

    

      SEACell-39
      Neurons
      0.428571
    

    

      SEACell-6
      Neurons
      1.000000
    

    

      SEACell-77
      GlutamatergicNeurons_early
      0.796875
    

  



























In [33]:



     


print(SEA_adata.X[40:43, 40:43])



     




















    
    






  (0, 0)	1.0
  (1, 0)	5.0
  (1, 1)	2.0
  (2, 0)	4.0




























7. Visualize results¶
    

  • plot_2D():  visualize metacell assignments on UMAP or any 2-dimensional embedding froma adata.obsm. Plots can also be coloured by metacell assignment. The visualization allows to assess the representativeness of the calculated metacells: ability to represent the global structure of the original dataset. Better representation corresponds to a more uniform coverage of the dataset.
    

    • 

  • plot_SEACell_sizes(): distribution of number of cells assigned to each metacell. Given to the non-uniform distribution of cell type density, metacells are expected to vary in size: larger metacells are retrieved from denser regions, and smaller metacells from sparser ones. However outlier values may require further inspection.
    

    • 




















In [34]:



     


SEACells.plot.plot_2D(adata, key='X_umap', colour_metacells=False)



     




















    
    






/usr/local/lib/python3.8/dist-packages/seaborn/relational.py:573: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  points = ax.scatter(x=x, y=y, **kws)
/usr/local/lib/python3.8/dist-packages/seaborn/relational.py:573: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  points = ax.scatter(x=x, y=y, **kws)








    
    































In [35]:



     


SEACells.plot.plot_2D(adata, key='X_umap', colour_metacells=True)



     




















    
    






/usr/local/lib/python3.8/dist-packages/seaborn/relational.py:573: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  points = ax.scatter(x=x, y=y, **kws)
/usr/local/lib/python3.8/dist-packages/seaborn/relational.py:573: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
  points = ax.scatter(x=x, y=y, **kws)








    
    































In [36]:



     


SEACells.plot.plot_SEACell_sizes(adata, bins=10, figsize=(10, 5))



     




















    
    






/usr/local/lib/python3.8/dist-packages/SEACells/plot.py:121: UserWarning: 

`distplot` is a deprecated function and will be removed in seaborn v0.14.0.

Please adapt your code to use either `displot` (a figure-level function with
similar flexibility) or `histplot` (an axes-level function for histograms).

For a guide to updating your code to use the new functions, please see
https://gist.github.com/mwaskom/de44147ed2974457ad6372750bbe5751

  sns.distplot(label_df.groupby('SEACell').count().iloc[:, 0], bins=bins)








    
    















    
    
Out[36]:











  
    

      
      size
    

    

      SEACell
      
    

  
  
    

      SEACell-0
      17
    

    

      SEACell-1
      100
    

    

      SEACell-10
      21
    

    

      SEACell-11
      110
    

    

      SEACell-12
      112
    

    

      ...
      ...
    

    

      SEACell-8
      65
    

    

      SEACell-80
      69
    

    

      SEACell-81
      28
    

    

      SEACell-82
      96
    

    

      SEACell-9
      29
    

  



83 rows × 1 columns






























8. Quantifying results¶
    

  • Purity: compute_celltype_purity(ad, col_name) computes the purity of different celltype labels within a SEACell metacell.
    

    • 

  • Compactness: per-SEAcell variance in diffusion components (typically 'X_pca' for RNA). Lower values of compactness suggest more compact/lower variance metacells.
    

    • 

  • Separation: distance between a SEACell and its nth_nbr. If cluster is provided as a string, e.g. 'celltype', nearest neighbors are restricted to have the same celltype value. Higher values of separation suggest better distinction between metacells.
    

    • 
























8.1 Purity¶
Cell Population¶


















In [37]:



     


SEACell_purity = SEACells.evaluate.compute_celltype_purity(adata, 'leidenAnnotated')

plt.figure(figsize=(6,4))
sns.boxplot(data=SEACell_purity, y='leidenAnnotated_purity')
plt.title('leidenAnnotated Purity')
sns.despine()
plt.show()
plt.close()

SEACell_purity.head()



     




















    
    















    
    
Out[37]:











  
    

      
      leidenAnnotated
      leidenAnnotated_purity
    

    

      SEACell
      
      
    

  
  
    

      SEACell-0
      ProliferatingProgenitors
      1.000000
    

    

      SEACell-1
      MigratingNeurons
      0.660000
    

    

      SEACell-10
      ProliferatingProgenitors
      1.000000
    

    

      SEACell-11
      GlutamatergicNeurons_early
      0.618182
    

    

      SEACell-12
      Neurons
      1.000000
    

  































8.2 Compactness¶


















In [38]:



     


compactness = SEACells.evaluate.compactness(adata, 'X_pca')

plt.figure(figsize=(6,4))
sns.boxplot(data=compactness, y='compactness')
plt.title('Compactness')
sns.despine()
plt.show()
plt.close()

compactness.head()



     




















    
    






Determing nearest neighbor graph...
computing neighbors
    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:00)








    
    















    
    
Out[38]:











  
    

      
      compactness
    

    

      SEACell
      
    

  
  
    

      SEACell-0
      0.008562
    

    

      SEACell-1
      0.001151
    

    

      SEACell-10
      0.023260
    

    

      SEACell-11
      0.002424
    

    

      SEACell-12
      0.001706
    

  































8.3 Separation¶


















In [39]:



     


separation = SEACells.evaluate.separation(adata, 'X_pca')

plt.figure(figsize=(6,4))
sns.boxplot(data=separation, y='separation')
plt.title('Separation')
sns.despine()
plt.show()
plt.close()

separation.head()



     




















    
    






Determing nearest neighbor graph...
computing neighbors
    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:00)








    
    















    
    
Out[39]:











  
    

      
      separation
    

    

      SEACell
      
    

  
  
    

      SEACell-0
      0.226287
    

    

      SEACell-1
      0.098583
    

    

      SEACell-10
      0.195324
    

    

      SEACell-11
      0.090437
    

    

      SEACell-12
      0.084924
    

  































9. Hard assignment: downstream processing¶
    

  • normalization and log-transformation
    • 

  • dimensionality reduction: PCA and UMAP on SEACells
    • 




















In [40]:



     


SEA_adata.layers['counts'] = SEA_adata.X.copy()

sc.pp.normalize_total(SEA_adata)
sc.pp.log1p(SEA_adata)
SEA_adata.layers['lognorm'] = SEA_adata.X.copy()

sc.pp.highly_variable_genes(SEA_adata, inplace=True)
sc.tl.pca(SEA_adata, use_highly_variable=True)
sc.pp.neighbors(SEA_adata, use_rep='X_pca')
sc.tl.umap(SEA_adata)



     




















    
    






normalizing counts per cell
    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)
computing PCA
    on highly variable genes
    with n_comps=50
    finished (0:00:00)
computing neighbors
    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:00)
computing UMAP
    finished: added
    'X_umap', UMAP coordinates (adata.obsm) (0:00:00)
























In [41]:



     


sc.pl.umap(SEA_adata, color=['leidenAnnotated'], s=75)



     




















    
    






/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 [42]:



     


sc.pl.umap(SEA_adata, color=['leidenAnnotated_purity'], s=75)



     




















    
    



































10. Explore hormonal gene expressions¶


















In [17]:



     


outdir = "../../../../../FigPaper/"

os.makedirs(outdir, exist_ok=True)

sc.settings.figdir = outdir

gene_dict = sig.groupby("Signature")["GeneName"].apply(list).to_dict()


gene_dict_filtered = {
    sig_name: [g for g in genes if g in adata.var_names]
    for sig_name, genes in gene_dict.items()
}

gene_dict_filtered = {k: v for k, v in gene_dict_filtered.items() if v}

if gene_dict_filtered:
    sc.pl.dotplot(
        SEA_adata,
        gene_dict_filtered,
        groupby='leidenAnnotated',
        save="dotplot_Castaldi_day50.png"
    )
else:
    print("None of the specified genes are found.")



     




















    
    






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








    
    






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








    
    





































11. Saving¶






















11.1 Save SEACell Adata¶


















In [44]:



     


SEA_adata



     




















    
    
Out[44]:








AnnData object with n_obs × n_vars = 83 × 33538
    obs: 'leidenAnnotated', 'leidenAnnotated_purity'
    var: 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
    uns: 'log1p', 'hvg', 'pca', 'neighbors', 'umap', 'leidenAnnotated_colors'
    obsm: 'X_pca', 'X_umap'
    varm: 'PCs'
    layers: 'raw', 'counts', 'lognorm'
    obsp: 'distances', 'connectivities'
























In [45]:



     


SEA_adata.X



     




















    
    
Out[45]:








<83x33538 sparse matrix of type '<class 'numpy.float64'>'
	with 928181 stored elements in Compressed Sparse Row format>
























In [46]:



     


SEA_adata.write(output_file)



     























11.2 Save in other formats¶


















In [47]:



     


nb_fname = ipynbname.name()
nb_fname



     




















    
    
Out[47]:








'SEACellsCastaldi-day50'
























In [50]:



     


%%bash -s "$nb_fname"
jupyter nbconvert "$1".ipynb --to="python"
jupyter nbconvert "$1".ipynb --to="html"



     




















    
    






[NbConvertApp] Converting notebook SEACellsCastaldi-day50.ipynb to python
[NbConvertApp] Writing 13000 bytes to SEACellsCastaldi-day50.py
[NbConvertApp] Converting notebook SEACellsCastaldi-day50.ipynb to html
[NbConvertApp] Writing 21710081 bytes to SEACellsCastaldi-day50.html




























11.3 Finished computations: timestamp¶


















In [49]:



     


print(datetime.now())



     




















    
    






2026-01-08 13:39:01.597743