| id | gene_panel_selection |
| name | Gene Panel Selection Workflow |
| description | End-to-end workflow for gene panel design in scRNA-seq and spatial transcriptomics, that should be **STRICTLY** followed:
dataset understanding + smart downsampling + train/test splits,
algorithmic selection (HVG/DE/RF/scGeneFit/SpaPROS),
optimal sub-panel discovery (ARI vs size),
biological completion with a stability gate (Completion Rule), consensus scoring and completion (only if there is still room),
and benchmarking on test splits (ARI/NMI/Silhouette + UMAP similarity).
|
| tags | ["gene-panel","selection","scrna-seq","spatial","scanpy","scverse","benchmarking","spapros","scgenefit","random-forest"] |
Gene Panel Selection Workflow
This skill is used when you need to construct biologically meaningful and algorithmically robust gene panels. You will receive context from the leaderagent , use this context and STRICLTY follow this Gene Panel Selection Workflow
Workflow Enforcement (MANDATORY)
Determine which stage of the workflow (Steps 1–5) is required for the current task,
and STRICTLY follow the corresponding step(s).
Once a step is entered, all its mandatory sub-steps must be executed.
No partial execution or silent degradation is allowed.
Workdir
Always work in the workdir provided by the leader agent.
Calling other agents
You can call other agents by using the call_agent(agent_name, instruction) function.
-
Call the browser_use agent for information collection:
When you encounter software or biological knowledge you are not familiar with, call browser_use to search the web and collect the necessary information.
-
Call the system_manager agent for software environment installation:
When you need to install software packages, call system_manager to install them.
-
Call the biologist agent for results interpretation:
When you plot figures, compute a panel, or have intermediate results, call biologist to ask for interpretations and include them in your report.
-
Call the reporteragent when all the results are obtained, to make a well written pdf.
Visual understanding
Use the observe_images function in the file_manager toolset to examine images and figures.
If a figure is not publication-quality, replot it.
Reporting
At the end of the task, write a markdown report named:
report_analysis.md
The report must include:
- Summary
- Data (inputs, key parameters, outputs)
- Results (figures + tables)
- Key findings
- Next steps
Large datasets
If the dataset is large, perform smart downsampling while preserving all cell types.
Gene Panel Selection Hyperparameters
Constant (in scripts/gene_panel_helpers.py) | Default | Description |
|---|
SCGENEFIT_MAX_CONSTRAINTS | 1000 | Max LP constraints for scGeneFit |
SPAPROS_N_HVG | 3000 | HVG pre-filter size for SpaPROS |
RF_N_ESTIMATORS | 300 | Random Forest tree count |
SPAPROS_RUNTIME_WARNING_MINUTES | 5.0 | Leader prompts the user via notify_user if SpaPROS estimate exceeds this |
SPAPROS_RUNTIME_SKIP_MINUTES | 30.0 | Strong-warning threshold; user may choose to skip SpaPROS |
ARI_DROP_THRESHOLD | 0.05 | Max acceptable ARI degradation during panel completion |
DOWNSAMPLE_MAX_CELLS | 500000 | Above this cell count, downsampling is required before selection |
GENE_COUNT_THRESHOLD | 30000 | Above this gene count, gene subsetting is required before selection |
N_TRAINING_SPLITS | 1 | Number of training datasets to build during the train/test split |
N_TEST_SPLITS | 5 | Minimum number of test splits to build (more is fine) |
SPLIT_CELL_LIMIT | 50000 | Target cells per test split (soft cap, preserve diversity) |
Workflow (IMPORTANT : STRICLY FOLLOW NEEDED STEPS)
0. Dataset
If the user provided an AnnData object / dataset path → skip to Step 1.
If no dataset was provided, you must search and retrieve a relevant dataset
before proceeding. Follow the sub-steps below in order.
[!IMPORTANT]
Before starting, read the database access skill index:
.pantheon/skills/omics/database_access/SKILL.md
(or use glob with pattern="**/database_access/SKILL.md")
and the relevant skill files it references (especially cellxgene_census.md and gget.md).
0.1 Parse the user query
Extract search parameters from the leader-provided context:
- Organism: e.g., "Homo sapiens", "Mus musculus"
- Tissue / organ: e.g., "lung", "brain", "bone marrow", "tumor"
- Disease context: e.g., "COVID-19", "cancer", "normal"
- Cell types of interest: e.g., "immune cells", "T cells", "neurons"
- Assay preference: e.g., scRNA-seq, spatial transcriptomics
- Scope: Is this a focused task (single tissue/disease/system) or a broad task
(multi-tissue, pan-disease, cross-system)? This is critical for dataset selection.
[!CRITICAL]
Match dataset scope to task scope. The dataset(s) you retrieve must be
representative of the full biological diversity the panel is designed for.
- Focused task (e.g., "brain cortex panel", "kidney disease panel") →
fetch data from that specific tissue/disease/system
- Broad / cross-tissue task (e.g., "pan-cancer panel", "whole-body immune panel",
"multi-organ developmental panel") → you must include data from all relevant
tissues, diseases, or biological contexts so the panel captures both shared and
context-specific biology. Do NOT narrow down to a single tissue or disease.
- In general: the biological diversity in the retrieved dataset should reflect the
biological diversity that the final gene panel must resolve. If the panel needs
to distinguish 10 tissues, the dataset must contain cells from those 10 tissues.
0.2 Search CELLxGENE Census (PRIMARY source)
CELLxGENE Census is the largest curated single-cell collection (217M+ cells)
and returns AnnData objects directly — always try this first.
Read the full skill: .pantheon/skills/omics/database_access/cellxgene_census.md
Strategy:
A) First, look for existing atlases / large integrated datasets that already match
the task scope. CELLxGENE hosts many curated cross-tissue and disease-specific atlases
(e.g., Tabula Sapiens, Human Cell Atlas collections, organ-specific atlases,
disease-focused atlases). A single well-curated atlas is far better than stitching
together cells from separate studies (avoids batch effects, inconsistent annotations, etc.).
import cellxgene_census
with cellxgene_census.open_soma() as census:
datasets = census["census_info"]["datasets"].read().concat().to_pandas()
print(datasets[["dataset_id", "collection_name", "dataset_title"]].head(30))
obs_df = cellxgene_census.get_obs(
census, "<organism>",
value_filter="is_primary_data == True",
column_names=["dataset_id", "tissue_general", "disease", "cell_type"],
)
diversity = obs_df.groupby("dataset_id").agg(
n_cells=("cell_type", "size"),
n_tissues=("tissue_general", "nunique"),
n_diseases=("disease", "nunique"),
n_cell_types=("cell_type", "nunique"),
).sort_values("n_tissues", ascending=False)
print(diversity.head(20))
Pick the dataset that best matches the task scope. For broad tasks, prioritize datasets
with highest tissue/disease/cell-type diversity. For focused tasks, prioritize relevance
to the specific tissue/disease. Prefer datasets with >50k cells and existing cell type annotations.
B) If no single atlas suffices, build a composite query across multiple tissues/diseases:
-
Explore available data — query cell metadata to estimate dataset sizes:
with cellxgene_census.open_soma() as census:
obs_df = cellxgene_census.get_obs(
census, "<organism>",
value_filter="tissue_general == '<tissue>' and is_primary_data == True",
column_names=["cell_type", "tissue", "tissue_general", "disease", "assay", "dataset_id"],
)
print(f"Total cells: {len(obs_df)}")
print(obs_df["cell_type"].value_counts().head(20))
print(obs_df["disease"].value_counts().head(10))
print(obs_df["tissue_general"].value_counts().head(15))
-
Refine filters — but preserve the task scope:
- For broad tasks: keep multiple tissues/diseases/contexts in the filter.
Sample a balanced number of cells per category to avoid one dominating.
- For focused tasks: narrow to the specific tissue/disease/context.
- Always check the diversity of cell types, tissues, and diseases after filtering
to confirm the dataset matches the task scope.
-
Download the dataset as AnnData:
with cellxgene_census.open_soma() as census:
adata = cellxgene_census.get_anndata(
census,
organism="<organism>",
obs_value_filter="<refined_filter> and is_primary_data == True",
column_names={
"obs": ["cell_type", "tissue", "tissue_general", "disease", "sex",
"assay", "donor_id", "dataset_id", "development_stage"],
},
)
-
Always filter is_primary_data == True to avoid duplicate cells
-
If the dataset is very large (above DOWNSAMPLE_MAX_CELLS, default 500000), downsample per category rather than
dropping entire tissues/diseases. For example, sample up to N cells per
(tissue, disease) combination to keep diversity while controlling size.
Alternatively, use the streaming API (ExperimentAxisQuery) — see the skill file.
0.3 Alternative sources (if Census is insufficient)
If CELLxGENE Census does not have suitable data (e.g., rare tissue, specific organism,
spatial data needed), try these alternatives in order of preference:
- gget.cellxgene — query CZ CELLxGENE Discover for specific datasets:
Read:
.pantheon/skills/omics/database_access/gget.md
import gget
gget.setup("cellxgene")
adata = gget.cellxgene(species="homo_sapiens", tissue="<tissue>",
cell_type=["<cell_type1>", "<cell_type2>"])
- GEO / ArrayExpress — call
browser_use to search for accession numbers,
then download via gget or direct URL
- Human Cell Atlas (HCA) / Tabula Sapiens / Broad Single Cell Portal
— call
browser_use for specific dataset URLs
Prefer datasets that already provide processed count matrices
(h5ad, loom, mtx format) with cell type annotations and metadata.
0.4 Validate the retrieved dataset
Before proceeding to Step 1, verify:
[!NOTE]
Document in the notebook which database was queried, what filters were used,
and why this dataset was selected. This information goes into the final report (Step 6).
1) Dataset Understanding and Splitting
Start with exploratory inspection using an integrated notebook.
1.1 Basic structure
Inspect:
- file format (h5ad or other)
- number of cells / genes
- batches / conditions
.obs, .var, .obsm, .uns
- whether dataset has spatial or multimodal components
Checklist:
1.2 Downsampling (CRITICAL)
Thresholds come from the module-level constants in
scripts/gene_panel_helpers.py (see the Hyperparameters table at the top).
Rules:
- If
adata.n_obs > DOWNSAMPLE_MAX_CELLS (default 500000): downsample to below that limit, preserving all cell types.
- If
adata.n_vars > GENE_COUNT_THRESHOLD (default 30000): reduce to ≤ GENE_COUNT_THRESHOLD via QC/HVG for compute-heavy steps.
- Save downsampled
adata to a new file in workdir via file_manager.
[!IMPORTANT]
Prefer stratified downsampling by label_key if available; otherwise stratify by clustering labels.
Use the constants imported from gene_panel_helpers (or the numeric defaults documented above) — don't re-hardcode them elsewhere.
1.3 Splitting
If provided one dataset, split to preserve all cell type distribution across all datasets:
N_TRAINING_SPLITS training dataset(s), diversified (default: 1).
- at least
N_TEST_SPLITS test batches (default: 5).
- constraint: each split should target
SPLIT_CELL_LIMIT cells (default: 50000) to preserve diversity — treat this as a soft cap, go slightly under rather than well under.
- make splits as non-redundant as possible and represent all cell types.
1.4 Disk Space Management (MANDATORY)
[!CRITICAL]
Gene panel selection generates large intermediate files. You must minimize disk usage
to avoid running out of space mid-pipeline.
Rules:
1.5 Preprocessing status
Check:
- normalization
- PCA
- UMAP
- clustering
Recompute only if missing or invalid.
1.6 Preprocessing (if needed)
- QC (follow the QC skill if available)
- Normalize / log1p / scale
- PCA / neighbors / UMAP
- Batch correction (if needed)
- Leiden clustering
- DEG & marker detection
- Cell type annotation
- Marker plots (dotplots, heatmaps)
[!IMPORTANT]
If heavy steps are slow or unstable on notebook use python code
2) Algorithmic Gene Panel Selection
2.1 Pre-established methods
Algorithmic Methods = {HVG, DE, Random Forest, scGeneFit, SpaPROS}
- Use true cell type as
label_key whenever available.
- Implement HVG / DE via Scanpy directly (
sc.pp.highly_variable_genes,
sc.tl.rank_genes_groups).
- For Random Forest / scGeneFit / SpaPROS use the helper script shipped
with this skill:
scripts/gene_panel_helpers.py. It is a plain Python
module (no registered toolset) with four functions:
select_spapros, select_random_forest, select_scgenefit,
estimate_spapros_runtime.
Load the helper script
Open the helper from this skill's directory (.pantheon/skills/omics/gene_panel_selection/scripts/).
In a notebook cell:
import sys
from pathlib import Path
skill_scripts = Path.cwd()
for candidate in (
Path.home() / ".pantheon/skills/omics/gene_panel_selection/scripts",
Path.cwd() / ".pantheon/skills/omics/gene_panel_selection/scripts",
):
if candidate.exists():
skill_scripts = candidate
break
sys.path.insert(0, str(skill_scripts))
from gene_panel_helpers import (
estimate_spapros_runtime,
select_spapros,
select_random_forest,
select_scgenefit,
SCGENEFIT_MAX_CONSTRAINTS,
SPAPROS_N_HVG,
RF_N_ESTIMATORS,
SPAPROS_RUNTIME_WARNING_MINUTES,
SPAPROS_RUNTIME_SKIP_MINUTES,
)
If the above lookup fails (non-default install), read the file via
skill_view(name='omics/gene_panel_selection', file_path='scripts/gene_panel_helpers.py')
to discover its location, then adjust the sys.path insertion accordingly.
[!CAUTION]
SpaPROS runtime gate (MANDATORY). SpaPROS can run for tens of minutes
to hours on large datasets. Call estimate_spapros_runtime(...) first
and inspect the severity tier of the returned dict:
"fast" → run select_spapros(...) directly, no user confirmation needed.
"slow" or "very_slow" → stop, return the estimate dict verbatim
to the leader so it can call notify_user with a Run/Skip choice. The
leader then dispatches back with one of:
"SpaPROS APPROVED by user" → run select_spapros(...) as normal.
"SpaPROS SKIPPED by user" → do not call select_spapros;
report the skip in report_analysis.md and continue with the other
methods only.
estimate = estimate_spapros_runtime(
adata_path=adata_path,
num_markers=200,
n_hvg=SPAPROS_N_HVG,
warning_minutes=SPAPROS_RUNTIME_WARNING_MINUTES,
skip_minutes=SPAPROS_RUNTIME_SKIP_MINUTES,
)
print(estimate)
select_scgenefit(
adata_path=adata_path,
label_key="cell_type",
return_scores=True,
max_constraints=SCGENEFIT_MAX_CONSTRAINTS,
workdir=workdir,
)
select_spapros(
adata_path=adata_path,
label_key="cell_type",
num_markers=200,
n_hvg=SPAPROS_N_HVG,
return_scores=True,
workdir=workdir,
)
select_random_forest(
adata_path=adata_path,
label_key="cell_type",
n_estimators=RF_N_ESTIMATORS,
return_scores=True,
workdir=workdir,
)
- Always request gene scores (
return_scores=True).
- Outputs land in
workdir/gene_panels/{spapros,random_forest,scgenefit}/.
- Each method writes a scores CSV that is the single source of truth
for ranking — Step 3 consumes these files directly.
- To override any cap per call, pass the kwarg (
max_constraints=,
n_hvg=, n_estimators=, etc.); don't edit the defaults in the helper
unless the new value should apply project-wide.
3) Optimal SEED panel Discovery
For each method independently (HVG, DE, Scgenefit, RF, SpapROS):
Let N be the target final panel size requested by the leader.
[!CRITICAL]
The ARI vs K sweep is a pandas slicing operation, not a re-run of
the algorithm. Each method's scores CSV (from Step 2.1) already ranks
every gene. Call each algorithm once; slice top-K in memory.
- Load the method-specific gene score CSV and rank genes (descending score):
import pandas as pd
scores = pd.read_csv(scores_csv).sort_values("score", ascending=False)
- Build candidate sub-panels of sizes K ∈ {100, 200, …, N} by taking the top-K:
panel_K = scores.head(K)["gene"].tolist()
- For each method and each K:
- Subset the dataset to panel genes only:
adata_K = adata[:, panel_genes]
- Recompute neighbors + Leiden on
adata_K (same preprocessing policy across K)
- Compute ARI between Leiden clusters and true cell types (
label_key).
- Plot ARI vs K for each method.
- Pick the seed panel = (method, K*) with the best ARI.
Note: SEED STEP is performed using the training adata. It is IMPORTANT you investigate ARI vs panel size for all methods (HVG, DE, Scgenefit, RF, SpapROS) when possible, to make sure you take the best one!
4) Curation Logic
4.1 Curation pipeline (STRICT ORDER)
Final panel is built in two phases:
Phase 1 — Seed-panel (algorithmic)
- Use the optimal Seed-panel identified in Step 3 as seed subpanel
- Do not change genes in the seed
Phase 2 — Completion (biological lookup is the PRIMARY mechanism)
CRITICAL: Biological curation is the MAIN completion mechanism, NOT consensus fill.
The purpose of completion is to add biologically meaningful genes that algorithmic methods may have missed.
Consensus fill is ONLY a small last-resort gap filler. If you find yourself adding more consensus-fill genes
than biological genes, you have NOT done enough biological lookup.
0) Completion Rule
Before adding a batch of genes:
- test whether it makes ARI drop considerably or become less stable (training)
- If completing the panel up to size N degrades performance substantially (eg ARI drop >
ARI_DROP_THRESHOLD), propose:
- an optimal stable panel (< N)
- a supplemental gene list to reach N if required
- a modest ARI drop is acceptable if it adds important biological coverage
Check this on the training dataset.
1) Assess Seed Coverage First
Before biological lookup, inspect genes in the seed panel:
- Map seed gene IDs to symbols
- Identify which biological categories from the leader's context are already covered
- Note which categories are MISSING or under-represented
2) Exhaustive Biological Lookup (CRITICAL — MUST BE THOROUGH)
Derive the relevant biological categories from the leader-provided context (e.g., cell type markers, signaling pathways, functional states, disease-specific genes — whatever the user's goal requires).
Call browser_use MULTIPLE times, once per major biological category identified.
For each category, collect all well-established marker genes (typically 10-30+ per category, not just 3-5).
Sources: GeneCards, GO, UniProt, KEGG, Reactome, MSigDB, published marker gene lists, review articles.
A single browser_use call returning a handful of genes for an entire panel is INSUFFICIENT.
The number of biologically curated genes should scale with the gap between seed size and target N.
Do multiple rounds of lookup — breadth across ALL relevant categories AND depth within each.
3) Add Biologically Relevant Genes
For each candidate gene:
- check not already in seed panel
- ensure no redundancy
- maintain balanced biological coverage across categories
- categorize into a relevant biological category (from leader context, or inferred)
- after each batch of additions, check Completion Rule (ARI stability on training)
- if ARI drops sharply, try a different set; a modest drop for strong biological coverage is acceptable
- continue until all important biological genes are added or panel reaches size N
4) Consensus Fill (LAST RESORT ONLY — small gap filler)
Only if after exhaustive biological lookup, {seed + biological genes} < N:
- normalize scores per method (same scale, no method dominates)
- aggregate into a consensus table
- fill the small remaining gap by score priority, excluding genes already present
Deliverable: a gene × {method where it comes from, biological category, biological function, source/reference} table.
Note: Every accepted gene must be justified, assigned a biological category, and referenced with a source (seed/method score or website/literature) and a gene function description.
5) Benchmarking (MANDATORY)
5.0 Panel genes comparison
Create an UpSet plot for all N-size algorithmic panels to see overlap.
Use the full original dataset for evaluation.
5.1 Dataset
Benchmarking is performed on test datasets.
5.2 Metrics
For each subset compute (across test splits):
- all algorithmic N size panels
- final curated N size panel
- if curated N was not optimal per Completion Rule, also benchmark the optimal stable (<N) panel
- full gene set baseline
Compute:
- Leiden over-clustering on panel genes
- ARI, NMI between Leiden and true labels
- Silhouette Index using Leiden assignments
Plots:
- one figure per metric
- boxplots
- high-quality formatting
5.3 UMAP comparison
Compute UMAPs for:
- full genes (reference)
- each algorithmic N size panel
- final curated N size panel
- if needed, the optimal stable panel
Compare vs reference:
- qualitative
- quantitative (distance correlation / Procrustes-like metrics)
6) Summarizing
Report must include the full workflow (Steps 0 → 5) and at minimum, in a very well written pdf (ask reporter to make the pdf):
- Objective & context
- Dataset description (structure, labels, preprocessing status)
- Algorithmic methods run (HVG/DE/RF/scGeneFit/SpaPROS): what each optimizes (detailed)
- Sub-panel selection:
- ARI vs size curves per method
- UpSet plot of different panels (overlaps)
- selection decision (method + size) and why
- Consensus table construction:
- normalization choice
- aggregation rule
- resulting ranked list
- Curation & completion reasoning (step-by-step):
- per added gene: lookup → match to context → accept/reject
- redundancy checks + category balance
- all biological references
- Benchmarking results:
- UpSet plot comparing algorithmic panels and curated panel
- ARI/NMI/SI boxplots across test subsets
- UMAP comparisons + quantitative similarity metric
- interpretation of performance differences
Tables (MANDATORY)
- Recap table of final panel (all N genes):
| Gene | Methods where it appears | Biological Function | Relevance score |
|---|
- Per-category count recap table based on leader context.
Figures (MANDATORY)
The report should contain at least all of the following figures , and any other figures that you consider relevant:
- ARI vs size curves per method (See above Sub-panel selection)
- UpSet plot comparing algorithmic panels and curated panel (See above Benchmarking results)
- ARI/NMI/SI boxplots across test subsets (See above Benchmarking results)
- UMAP comparisons + quantitative similarity metric (See above Benchmarking results)
Guidelines for integrated notebook usage
Use the integrated_notebook toolset to create/manage/execute notebooks.
- Keep all related code in the same notebook
- Each notebook handles one specific analysis task
- Start each notebook with a markdown cell:
- After each code cell producing results, add a markdown cell explaining the result
- Save figures and also display them in notebook outputs
If memory becomes insufficient:
- close kernels using
manage_kernel
- reduce compute via stratified downsampling (preserve all cell types) and/or split heavy operations into separate cells
- document decisions explicitly (what was checked, what was changed, why)
[!IMPORTANT]
Do NOT lighten the computation or reduce the data to dodge a kernel crash — the point of the pipeline is to capture the full complexity of the dataset. If manage_kernel + splitting doesn't help, save the heavy code to a .py file and run it via python_interpreter, leaving a markdown cell in the notebook that points to the script path. This is a last resort — note it explicitly in report_analysis.md with a reason.
Visualization quality gate
We expect high-quality, publication-level figures.
After generating a figure:
- inspect via
observe_images
- if not good → replot
High-quality means:
- clear, readable
- labeled axes
- good color/contrast
- informative title (not too long)
If figure is not satisfactory → replot