{
  "skill": "alterlab-cellxgene",
  "evals": [
    {
      "id": "query-marker-genes-by-celltype",
      "prompt": "From the CZ CELLxGENE Census, pull human T cells and B cells but only the expression of CD4, CD8A, and CD19, and return it as an AnnData I can work with. Make sure I'm not double-counting cells.",
      "expected_output": "Invokes alterlab-cellxgene: opens the census with cellxgene_census.open_soma(), calls get_anndata with organism='Homo sapiens', a var_value_filter on feature_name for the three genes, and an obs_value_filter on cell_type that includes is_primary_data == True to avoid duplicates.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "is_primary_data" },
        { "type": "behavior", "value": "Uses cellxgene_census.get_anndata with var_value_filter on feature_name and an obs_value_filter, and includes is_primary_data == True." }
      ]
    },
    {
      "id": "explore-cell-types-in-tissue",
      "prompt": "Before I download anything, I just want to know what cell types are available in human brain in the Census and roughly how many cells of each. What's the lightweight way to check?",
      "expected_output": "Invokes alterlab-cellxgene: uses the explore-then-query pattern with cellxgene_census.get_obs (not get_anndata) on tissue_general == 'brain' and is_primary_data == True, selecting only column_names=['cell_type'] and doing value_counts to summarize without loading expression.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "get_obs" },
        { "type": "behavior", "value": "Queries metadata only via get_obs and value_counts before any expression download, filtering for primary data." }
      ]
    },
    {
      "id": "out-of-core-large-query",
      "prompt": "I need mean expression of FOXP2, TBR1, and SATB2 across every primary brain cell in the Census, but that's way more cells than my 48 GB machine can hold. How do I compute it without loading everything into memory?",
      "expected_output": "Invokes alterlab-cellxgene: uses out-of-core processing via axis_query() with soma.AxisQuery obs/var filters, iterates query.X('raw').tables() in batches accumulating sum and count incrementally to compute the mean, rather than get_anndata.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "axis_query" },
        { "type": "behavior", "value": "Recommends out-of-core iteration over X tables/batches instead of loading the full query into an AnnData." }
      ]
    },
    {
      "id": "reproducible-census-version",
      "prompt": "I'm writing the methods section of a paper that uses Census data. How do I make sure my exact same query is reproducible months from now when the Census has been updated?",
      "expected_output": "Invokes alterlab-cellxgene: advises pinning the release by passing census_version (e.g. open_soma(census_version='2025-11-08'), the current LTS build) so the analysis is reproducible across Census updates, and recording that version in the methods.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "census_version" },
        { "type": "behavior", "value": "Recommends explicitly pinning census_version for reproducibility rather than relying on the default 'stable' release." }
      ]
    },
    {
      "id": "near-miss-scanpy",
      "prompt": "I already loaded a lung dataset into an AnnData object. Normalize it, find highly variable genes, run PCA and UMAP, and cluster it with Leiden so I can see the cell populations.",
      "expected_output": "Should NOT trigger alterlab-cellxgene. This is downstream single-cell analysis (normalize, HVG, PCA, UMAP, Leiden) on an already-loaded dataset, which is alterlab-scanpy territory. cellxgene-census is for pulling reference data from the public atlas and defers analysis of your own dataset to scanpy.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "scanpy" }
      ]
    },
    {
      "id": "near-miss-anndata",
      "prompt": "I have three local .h5ad files on disk and I want to concatenate them along the cells axis with an inner join on genes and a batch label. No public data involved.",
      "expected_output": "Should NOT trigger alterlab-cellxgene. Concatenating local .h5ad files with ad.concat (axis, join, label) is alterlab-anndata's data-format territory. cellxgene-census is specifically for querying the public CZ CELLxGENE Census, not wrangling local files.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "anndata" }
      ]
    }
  ]
}

# Best Practices, Metadata Fields, and Troubleshooting

## Key Concepts and Best Practices

### Always Filter for Primary Data
Unless analyzing duplicates, always include `is_primary_data == True` in queries to avoid counting cells multiple times:
```python
obs_value_filter="cell_type == 'B cell' and is_primary_data == True"
```

### Specify Census Version for Reproducibility
Pin the Census version in anything you will report on, and record it in your methods:
```python
census = cellxgene_census.open_soma(census_version="2025-11-08")  # an LTS build
```
`"stable"` resolves to the most recent **LTS** release (2025-11-08 as of 2026-09) and
`"latest"` to the weekly build, so both move under you. CZI commits to keeping LTS releases
publicly available for at least five years; the LTS series so far is 2023-05-15, 2023-07-25,
2023-12-15, 2024-07-01, 2025-01-30 and 2025-11-08. Cell counts differ between them — read
`census["census_info"]["summary"]` for the build you actually opened.

### Estimate Query Size Before Loading
For large queries, first check the number of cells to avoid memory issues:
```python
# Get cell count
metadata = cellxgene_census.get_obs(
    census, "homo_sapiens",
    value_filter="tissue_general == 'brain' and is_primary_data == True",
    column_names=["soma_joinid"]
)
n_cells = len(metadata)
print(f"Query will return {n_cells:,} cells")

# If too large (>100k), use out-of-core processing
```

### Use tissue_general for Broader Groupings
The `tissue_general` field provides coarser categories than `tissue`, useful for cross-tissue analyses:
```python
# Broader grouping
obs_value_filter="tissue_general == 'immune system'"

# Specific tissue
obs_value_filter="tissue == 'peripheral blood mononuclear cell'"
```

### Select Only Needed Columns
Minimize data transfer by specifying only required metadata columns:
```python
obs_column_names=["cell_type", "tissue_general", "disease"]  # Not all columns
```

### Check Dataset Presence for Gene-Specific Queries
When analyzing specific genes, verify which datasets measured them:
```python
presence = cellxgene_census.get_presence_matrix(
    census,
    "homo_sapiens",
    var_value_filter="feature_name in ['CD4', 'CD8A']"
)
```

### Two-Step Workflow: Explore Then Query
First explore metadata to understand available data, then query expression:
```python
# Step 1: Explore what's available
metadata = cellxgene_census.get_obs(
    census, "homo_sapiens",
    value_filter="disease == 'COVID-19' and is_primary_data == True",
    column_names=["cell_type", "tissue_general"]
)
print(metadata.value_counts())

# Step 2: Query based on findings
adata = cellxgene_census.get_anndata(
    census=census,
    organism="Homo sapiens",
    obs_value_filter="disease == 'COVID-19' and cell_type == 'T cell' and is_primary_data == True",
)
```

## Available Metadata Fields

### Cell Metadata (obs)
Key fields for filtering:
- `cell_type`, `cell_type_ontology_term_id`
- `tissue`, `tissue_general`, `tissue_ontology_term_id`
- `disease`, `disease_ontology_term_id`
- `assay`, `assay_ontology_term_id`
- `donor_id`, `sex`, `self_reported_ethnicity`
- `development_stage`, `development_stage_ontology_term_id`
- `dataset_id`
- `is_primary_data` (Boolean: True = unique cell)

### Gene Metadata (var)
- `feature_id` (Ensembl gene ID, e.g., "ENSG00000161798")
- `feature_name` (Gene symbol, e.g., "FOXP2")
- `feature_length` (Gene length in base pairs)

## Troubleshooting

### Query Returns Too Many Cells
- Add more specific filters to reduce scope
- Use `tissue` instead of `tissue_general` for finer granularity
- Filter by specific `dataset_id` if known
- Switch to out-of-core processing for large queries

### Memory Errors
- Reduce query scope with more restrictive filters
- Select fewer genes with `var_value_filter`
- Use out-of-core processing with `axis_query()`
- Process data in batches

### Duplicate Cells in Results
- Always include `is_primary_data == True` in filters
- Check if intentionally querying across multiple datasets

### Gene Not Found
- Verify gene name spelling (case-sensitive)
- Try Ensembl ID with `feature_id` instead of `feature_name`
- Check dataset presence matrix to see if gene was measured
- Some genes may have been filtered during Census construction

### Version Inconsistencies
- Always specify `census_version` explicitly
- Use same version across all analyses
- Check release notes for version-specific changes

# CZ CELLxGENE Census Data Schema Reference

## Overview

The CZ CELLxGENE Census is a versioned collection of single-cell data built on the TileDB-SOMA framework. This reference documents the data structure, available metadata fields, and query syntax.

## High-Level Structure

The Census is organized as a `SOMACollection` with two main components:

### 1. census_info
Summary information including:
- **summary**: Build date, cell counts, dataset statistics
- **datasets**: All datasets from CELLxGENE Discover with metadata
- **summary_cell_counts**: Cell counts stratified by metadata categories

### 2. census_data
Organism-specific `SOMAExperiment` objects:
- **"homo_sapiens"**: Human single-cell data
- **"mus_musculus"**: Mouse single-cell data

## Data Structure Per Organism

Each organism experiment contains:

### obs (Cell Metadata)
Cell-level annotations stored as a `SOMADataFrame`. Access via:
```python
census["census_data"]["homo_sapiens"].obs
```

### ms["RNA"] (Measurement)
RNA measurement data including:
- **X**: Data matrices with layers:
  - `raw`: Raw count data
  - `normalized`: (if available) Normalized counts
- **var**: Gene metadata
- **feature_dataset_presence_matrix**: Sparse boolean array showing which genes were measured in each dataset

## Cell Metadata Fields (obs)

### Required/Core Fields

**Identity & Dataset:**
- `soma_joinid`: Unique integer identifier for joins
- `dataset_id`: Source dataset identifier
- `is_primary_data`: Boolean flag (True = unique cell, False = duplicate across datasets)

**Cell Type:**
- `cell_type`: Human-readable cell type name
- `cell_type_ontology_term_id`: Standardized ontology term (e.g., "CL:0000236")

**Tissue:**
- `tissue`: Specific tissue name
- `tissue_general`: Broader tissue category (useful for grouping)
- `tissue_ontology_term_id`: Standardized ontology term

**Assay:**
- `assay`: Sequencing technology used
- `assay_ontology_term_id`: Standardized ontology term

**Disease:**
- `disease`: Disease status or condition
- `disease_ontology_term_id`: Standardized ontology term

**Donor:**
- `donor_id`: Unique donor identifier
- `sex`: Biological sex (male, female, unknown)
- `self_reported_ethnicity`: Ethnicity information
- `development_stage`: Life stage (adult, child, embryonic, etc.)
- `development_stage_ontology_term_id`: Standardized ontology term

**Organism:**
- `organism`: Scientific name (Homo sapiens, Mus musculus)
- `organism_ontology_term_id`: Standardized ontology term

**Technical:**
- `suspension_type`: Sample preparation type (cell, nucleus, na)

## Gene Metadata Fields (var)

Access via:
```python
census["census_data"]["homo_sapiens"].ms["RNA"].var
```

**Available Fields:**
- `soma_joinid`: Unique integer identifier for joins
- `feature_id`: Ensembl gene ID (e.g., "ENSG00000161798")
- `feature_name`: Gene symbol (e.g., "FOXP2")
- `feature_length`: Gene length in base pairs

## Value Filter Syntax

Queries use Python-like expressions for filtering. The syntax is processed by TileDB-SOMA.

### Comparison Operators
- `==`: Equal to
- `!=`: Not equal to
- `<`, `>`, `<=`, `>=`: Numeric comparisons
- `in`: Membership test (e.g., `feature_id in ['ENSG00000161798', 'ENSG00000188229']`)

### Logical Operators
- `and`, `&`: Logical AND
- `or`, `|`: Logical OR

### Examples

**Single condition:**
```python
value_filter="cell_type == 'B cell'"
```

**Multiple conditions with AND:**
```python
value_filter="cell_type == 'B cell' and tissue_general == 'lung' and is_primary_data == True"
```

**Using IN for multiple values:**
```python
value_filter="tissue in ['lung', 'liver', 'kidney']"
```

**Complex condition:**
```python
value_filter="(cell_type == 'neuron' or cell_type == 'astrocyte') and disease != 'normal'"
```

**Filtering genes:**
```python
var_value_filter="feature_name in ['CD4', 'CD8A', 'CD19']"
```

## Data Inclusion Criteria

The Census includes all data from CZ CELLxGENE Discover meeting:

1. **Species**: Human (*Homo sapiens*) or mouse (*Mus musculus*)
2. **Technology**: Approved sequencing technologies for RNA
3. **Count Type**: Raw counts only (no processed/normalized-only data)
4. **Metadata**: Standardized following CELLxGENE schema
5. **Both spatial and non-spatial data**: Includes traditional and spatial transcriptomics

## Important Data Characteristics

### Duplicate Cells
Cells may appear across multiple datasets. Use `is_primary_data == True` to filter for unique cells in most analyses.

### Count Types
The Census includes:
- **Molecule counts**: From UMI-based methods
- **Full-gene sequencing read counts**: From non-UMI methods
These may need different normalization approaches.

### Versioning
Census releases are versioned (e.g., "2025-11-08", "stable", "latest"). "stable" is the newest LTS build and "latest" the weekly one, so pin an explicit date for reproducible analysis:
```python
census = cellxgene_census.open_soma(census_version="2025-11-08")
```

## Dataset Presence Matrix

Access which genes were measured in each dataset:
```python
presence_matrix = census["census_data"]["homo_sapiens"].ms["RNA"]["feature_dataset_presence_matrix"]
```

This sparse boolean matrix helps understand:
- Gene coverage across datasets
- Which datasets to include for specific gene analyses
- Technical batch effects related to gene coverage

## SOMA Object Types

Core TileDB-SOMA objects used:
- **DataFrame**: Tabular data (obs, var)
- **SparseNDArray**: Sparse matrices (X layers, presence matrix)
- **DenseNDArray**: Dense arrays (less common)
- **Collection**: Container for related objects
- **Experiment**: Top-level container for measurements
- **SOMAScene**: Spatial transcriptomics scenes
- **obs_spatial_presence**: Spatial data availability

# Additional Patterns and Pitfalls

This file holds patterns that go beyond the core workflows. For the standard
query recipes see `querying_expression.md` (open / explore / `get_anndata` /
`axis_query`) and `ml_and_scanpy.md` (PyTorch dataloader, scanpy, multi-dataset
concatenation). The best-practices checklist lives in
`best_practices_and_troubleshooting.md`.

## Incremental Variance (Welford's Algorithm)

When computing variance out-of-core over a query too large for RAM, accumulate
it in one pass with Welford's online algorithm. Vectorize per batch rather than
looping over individual values:

```python
import numpy as np

n = 0
mean = 0.0
M2 = 0.0

for batch in query.X("raw").tables():
    values = batch["soma_data"].to_numpy()
    # Chan et al. parallel/batch update
    b_n = len(values)
    if b_n == 0:
        continue
    b_mean = values.mean()
    b_M2 = ((values - b_mean) ** 2).sum()
    delta = b_mean - mean
    new_n = n + b_n
    mean += delta * b_n / new_n
    M2 += b_M2 + delta**2 * n * b_n / new_n
    n = new_n

variance = M2 / (n - 1) if n > 1 else float("nan")
```

Note: `X("raw")` returns only the *stored* (nonzero) entries, so a mean/variance
computed this way is over observed nonzero values. To get statistics over the
full dense matrix (including implicit zeros), divide sums by
`n_cells * n_genes` from the query shape instead of by the number of nonzero
entries.

## Filter by Ontology Term for Cross-Dataset Consistency

Free-text labels (`cell_type == 'B cell'`) can vary subtly across datasets;
ontology term IDs are canonical and more reliable for population-scale queries:

```python
# Equivalent to 'B cell' but stable across all datasets
obs_value_filter="cell_type_ontology_term_id == 'CL:0000236' and is_primary_data == True"
```

## Batch Processing Across Conditions

For systematic analyses sweeping one factor, loop the query and collect results.
Pin `census_version` once and reuse the open handle:

```python
tissues = ["lung", "liver", "kidney", "heart"]
results = {}

with cellxgene_census.open_soma(census_version="2025-11-08") as census:
    for tissue in tissues:
        adata = cellxgene_census.get_anndata(
            census=census,
            organism="Homo sapiens",
            obs_value_filter=f"tissue_general == '{tissue}' and is_primary_data == True",
        )
        results[tissue] = analyze(adata)
```

## Common Pitfalls to Avoid

1. **Forgetting `is_primary_data == True`** — silently counts duplicate cells that appear in multiple datasets.
2. **Loading before estimating** — run a `get_obs` count first; switch to `axis_query` out-of-core above ~100k cells.
3. **Skipping the context manager** — leak open SOMA/TileDB handles.
4. **Unpinned version** — results drift when the `stable` release rolls forward; pin `census_version`.
5. **Overly broad queries** — start focused, then widen.
6. **Ignoring the presence matrix** — a gene absent from a dataset reads as all-zero, not missing; check `get_presence_matrix` before interpreting zeros.
7. **Mixing count types** — UMI molecule counts and full-length read counts coexist in `raw` and may need different normalization.
8. **Using deprecated APIs** — ML loaders are now in `tiledbsoma_ml`, not `cellxgene_census.experimental.ml`; concatenate with `anndata.concat`, not `adata.concatenate`.

# Machine Learning, Scanpy, and Multi-Dataset Integration

Code patterns for training PyTorch models on Census data, integrating with scanpy, and combining multiple datasets/tissues.

## Machine Learning with PyTorch

The PyTorch integration now lives in the standalone **`tiledbsoma_ml`** package
(`pip install tiledbsoma-ml`), NOT in `cellxgene_census.experimental.ml` — that
prototype API has been superseded. Build an `ExperimentDataset` from an
`axis_query`, then wrap it with `experiment_dataloader`:

```python
import torch
import tiledbsoma as soma
from tiledbsoma_ml import ExperimentDataset, experiment_dataloader

with cellxgene_census.open_soma(census_version="2025-11-08") as census:
    experiment = census["census_data"]["homo_sapiens"]
    with experiment.axis_query(
        measurement_name="RNA",
        obs_query=soma.AxisQuery(
            value_filter="tissue_general == 'liver' and is_primary_data == True"
        ),
    ) as query:
        dataset = ExperimentDataset(
            query,
            layer_name="raw",
            obs_column_names=["cell_type"],
            batch_size=128,
            shuffle=True,
            seed=42,
        )
        dataloader = experiment_dataloader(dataset)

        # Each batch is a (X, obs) tuple: X is a NumPy array, obs a pandas DataFrame.
        for epoch in range(num_epochs):
            for X_batch, obs_batch in dataloader:
                X = torch.from_numpy(X_batch).float()
                labels = label_encoder.transform(obs_batch["cell_type"])

                outputs = model(X)
                loss = criterion(outputs, torch.from_numpy(labels))

                optimizer.zero_grad()
                loss.backward()
                optimizer.step()
```

**Train/test splitting** — call `random_split` on the dataset (positional
fractions, not `split=[...]`), then wrap each split:
```python
train_dataset, test_dataset = dataset.random_split(0.8, 0.2, seed=42)
train_loader = experiment_dataloader(train_dataset)
test_loader = experiment_dataloader(test_dataset)
```

## Integration with Scanpy

Seamlessly integrate Census data with scanpy workflows:

```python
import scanpy as sc

# Load data from Census
adata = cellxgene_census.get_anndata(
    census=census,
    organism="Homo sapiens",
    obs_value_filter="cell_type == 'neuron' and tissue_general == 'cortex' and is_primary_data == True",
)

# Standard scanpy workflow
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000)

# Dimensionality reduction
sc.pp.pca(adata, n_comps=50)
sc.pp.neighbors(adata)
sc.tl.umap(adata)

# Visualization
sc.pl.umap(adata, color=["cell_type", "tissue", "disease"])
```

## Multi-Dataset Integration

Prefer a single query with an `in` filter (Strategy 2) — it pulls a consistent
gene set in one pass. Only query separately and concatenate when you need to tag
or transform each slice differently:

```python
import anndata as ad

# Strategy 1: Query multiple tissues separately, then concatenate
tissues = ["lung", "liver", "kidney"]
adatas = []
for tissue in tissues:
    a = cellxgene_census.get_anndata(
        census=census,
        organism="Homo sapiens",
        obs_value_filter=f"tissue_general == '{tissue}' and is_primary_data == True",
    )
    adatas.append(a)

# Use anndata.concat (adata.concatenate() is deprecated). Inner join keeps
# only genes shared across all slices.
combined = ad.concat(adatas, join="inner", label="tissue", keys=tissues)

# Strategy 2 (preferred): one query, multiple tissues
adata = cellxgene_census.get_anndata(
    census=census,
    organism="Homo sapiens",
    obs_value_filter="tissue_general in ['lung', 'liver', 'kidney'] and is_primary_data == True",
)
```

## Worked Use Cases

### Use Case 1: Explore Cell Types in a Tissue
```python
with cellxgene_census.open_soma() as census:
    cells = cellxgene_census.get_obs(
        census, "homo_sapiens",
        value_filter="tissue_general == 'lung' and is_primary_data == True",
        column_names=["cell_type"]
    )
    print(cells["cell_type"].value_counts())
```

### Use Case 2: Query Marker Gene Expression
```python
with cellxgene_census.open_soma() as census:
    adata = cellxgene_census.get_anndata(
        census=census,
        organism="Homo sapiens",
        var_value_filter="feature_name in ['CD4', 'CD8A', 'CD19']",
        obs_value_filter="cell_type in ['T cell', 'B cell'] and is_primary_data == True",
    )
```

### Use Case 3: Train Cell Type Classifier
```python
import tiledbsoma as soma
from tiledbsoma_ml import ExperimentDataset, experiment_dataloader

with cellxgene_census.open_soma(census_version="2025-11-08") as census:
    experiment = census["census_data"]["homo_sapiens"]
    with experiment.axis_query(
        measurement_name="RNA",
        obs_query=soma.AxisQuery(
            value_filter="tissue_general == 'blood' and is_primary_data == True"
        ),
    ) as query:
        dataset = ExperimentDataset(
            query, layer_name="raw", obs_column_names=["cell_type"],
            batch_size=128, shuffle=True, seed=42,
        )
        dataloader = experiment_dataloader(dataset)
        for X_batch, obs_batch in dataloader:
            ...  # Training logic
```

### Use Case 4: Cross-Tissue Analysis
```python
with cellxgene_census.open_soma() as census:
    adata = cellxgene_census.get_anndata(
        census=census,
        organism="Homo sapiens",
        obs_value_filter="cell_type == 'macrophage' and tissue_general in ['lung', 'liver', 'brain'] and is_primary_data == True",
    )

    # Analyze macrophage differences across tissues
    sc.tl.rank_genes_groups(adata, groupby="tissue_general")
```

# Querying Expression Data

Code patterns for opening the Census, exploring metadata, and retrieving expression matrices at small/medium and large scales.

## Opening the Census

Always use the context manager to ensure proper resource cleanup:

```python
import cellxgene_census

# Open latest stable version
with cellxgene_census.open_soma() as census:
    summary = census["census_info"]["summary"].read().concat().to_pandas()

# Open a specific version for reproducibility (preferred for published work)
with cellxgene_census.open_soma(census_version="2025-11-08") as census:
    summary = census["census_info"]["summary"].read().concat().to_pandas()
```

**Key points:**
- Use context manager (`with` statement) for automatic cleanup
- Specify `census_version` for reproducible analyses
- Default opens latest "stable" release

## Exploring Census Information

Before querying expression data, explore available datasets and metadata.

**Access summary information:**
```python
# Get summary statistics
summary = census["census_info"]["summary"].read().concat().to_pandas()
print(f"Total cells: {summary['total_cell_count'][0]}")

# Get all datasets
datasets = census["census_info"]["datasets"].read().concat().to_pandas()

# Filter datasets by criteria
covid_datasets = datasets[datasets["disease"].str.contains("COVID", na=False)]
```

**Query cell metadata to understand available data:**
```python
# Get unique cell types in a tissue
cell_metadata = cellxgene_census.get_obs(
    census,
    "homo_sapiens",
    value_filter="tissue_general == 'brain' and is_primary_data == True",
    column_names=["cell_type"]
)
unique_cell_types = cell_metadata["cell_type"].unique()
print(f"Found {len(unique_cell_types)} cell types in brain")

# Count cells by tissue
tissue_counts = cell_metadata.groupby("tissue_general").size()
```

**Important:** Always filter for `is_primary_data == True` to avoid counting duplicate cells unless specifically analyzing duplicates.

## Small-to-Medium Scale Queries (`get_anndata`)

For queries returning < 100k cells that fit in memory, use `get_anndata()`:

```python
# Basic query with cell type and tissue filters
adata = cellxgene_census.get_anndata(
    census=census,
    organism="Homo sapiens",  # or "Mus musculus"
    obs_value_filter="cell_type == 'B cell' and tissue_general == 'lung' and is_primary_data == True",
    obs_column_names=["assay", "disease", "sex", "donor_id"],
)

# Query specific genes with multiple filters
adata = cellxgene_census.get_anndata(
    census=census,
    organism="Homo sapiens",
    var_value_filter="feature_name in ['CD4', 'CD8A', 'CD19', 'FOXP3']",
    obs_value_filter="cell_type == 'T cell' and disease == 'COVID-19' and is_primary_data == True",
    obs_column_names=["cell_type", "tissue_general", "donor_id"],
)
```

**Filter syntax:**
- Use `obs_value_filter` for cell filtering
- Use `var_value_filter` for gene filtering
- Combine conditions with `and`, `or`
- Use `in` for multiple values: `tissue in ['lung', 'liver']`
- Select only needed columns with `obs_column_names`

**Getting metadata separately:**
```python
# Query cell metadata
cell_metadata = cellxgene_census.get_obs(
    census, "homo_sapiens",
    value_filter="disease == 'COVID-19' and is_primary_data == True",
    column_names=["cell_type", "tissue_general", "donor_id"]
)

# Query gene metadata
gene_metadata = cellxgene_census.get_var(
    census, "homo_sapiens",
    value_filter="feature_name in ['CD4', 'CD8A']",
    column_names=["feature_id", "feature_name", "feature_length"]
)
```

## Large-Scale Queries (Out-of-Core Processing)

For queries exceeding available RAM, use `axis_query()` with iterative processing:

```python
import tiledbsoma as soma

# Create axis query
query = census["census_data"]["homo_sapiens"].axis_query(
    measurement_name="RNA",
    obs_query=soma.AxisQuery(
        value_filter="tissue_general == 'brain' and is_primary_data == True"
    ),
    var_query=soma.AxisQuery(
        value_filter="feature_name in ['FOXP2', 'TBR1', 'SATB2']"
    )
)

# Iterate through expression matrix in chunks
iterator = query.X("raw").tables()
for batch in iterator:
    # batch is a pyarrow.Table with columns:
    # - soma_data: expression value
    # - soma_dim_0: cell (obs) coordinate
    # - soma_dim_1: gene (var) coordinate
    process_batch(batch)
```

**Computing incremental statistics:**
```python
# Example: Calculate mean expression
n_observations = 0
sum_values = 0.0

iterator = query.X("raw").tables()
for batch in iterator:
    values = batch["soma_data"].to_numpy()
    n_observations += len(values)
    sum_values += values.sum()

mean_expression = sum_values / n_observations
```