Feature selection: choosing informative features

Feature selection decides which measured genes define PCA and the neighborhood graph. Highly variable genes (HVGs) help focus the analysis on differences among cells while reducing the computational cost of using every measured gene. The goal is to retain informative variation without letting broadly expressed or noisy genes dominate the comparison.

CyteArc’s features.hvgs models the relationship between mean expression and variance, then selects genes whose corrected variance is high relative to genes with similar abundance. This matters because a highly expressed gene can have a large raw variance without distinguishing cell types.

This is distinct from cell quality control. A gene can be measured correctly and still be excluded because it contributes broad technical or confounding variation to the graph.

Fit the mean-variance model¶

The PBMC store carries an analysis saved as docs_default. It used 500 genes and 15 PCs to keep this teaching example small. The current API defaults are 1,000 genes and 21 PCs. We first repeat the saved 500-gene selection to inspect how those genes were chosen.

Open the downloaded store and its saved analysis.

Fit and inspect the 500-gene selection used by the saved example.

The plot should retain genes above the fitted mean-variance trend across a useful expression range. A selection concentrated only among the most abundant genes can make library size or housekeeping programs dominate the graph.

Understand the default exclusions¶

CyteArc applies the following case-insensitive regular expression to gene names:

^MT-|^RPS|^RPL|^MRPS|^MRPL|^CCN|^HLA-|^H2-|^HIST|
^XIST$|^DDX3Y$|^USP9Y$|^EIF1AY$|^KDM5D$|^SRY$|^ZFY$|^UTY$|^TMSB4Y$|^NLGN4Y$

Matching starts at the beginning of the name. Each family is defined once in cytearc.features.gene_families; ribosomal covers RPS, RPL, MRPS, and MRPL. Count how many genes in this dataset fall into each family:

These families can dominate broad variation associated with sequencing depth, cell quality, cell state, or sex without representing the cell identities sought in a typical heterogeneity workflow. Mitochondrial, ribosomal, and sex-linked genes can still be biologically relevant to your question. Treat the default as a starting point and decide which programs the graph should capture before excluding them.

These are name patterns, not a complete classification of gene function. For example, ^CCN matches both CCN1 and CCND1, while ^HIST matches HIST1H2AC but not H2AC1. Check the actual names removed in your gene annotation before relying on a family label. Supplying blacklist= replaces the whole default pattern.

By default, max_cells is n_selected - 20. With the default strict upper bound, genes detected in at least that many selected cells are excluded as nearly ubiquitous. The cutoff is skipped when too few selected cells remain to satisfy the minimum detection requirement.

Clearing the blacklist keeps every gene name while retaining other HVG filters. Compare the same top_n with and without the default pattern:

Other overrides stay available for study-specific work:

# Replace the default with a study-specific, case-insensitive regex.
ds.features.hvgs(cell_selection, blacklist=r"^MT-|^RPS|^RPL", top_n=2000)

# Disable the nearly ubiquitous gene filter.
ds.features.hvgs(cell_selection, max_cells=np.inf, top_n=2000)

Repeating a selection with the same cells and settings reuses the saved result, including the diagnostics used for its plot. See Reuse results and compare analysis settings for details.

Compare feature-set size¶

The number of selected genes changes the PCA basis and can change neighbourhood structure. This comparison keeps all other graph choices fixed. The 500-gene branch comes directly from docs_default; only the 1,000-gene branch is new.

Normalize the selected genes and construct the new graph.

Build the layout and clustering for the new graph.

The cross-tabulation compares the two partitions cell by cell. Cluster numbers are arbitrary labels, so a changed number does not mean the population changed: ask whether the same cells stayed together. A row spread across several columns shows a group splitting when more genes are included; a column collecting several rows shows groups merging. The margins report cluster sizes, helping distinguish a small boundary change from a large reorganization.

In the saved comparison, the 500-gene cluster 3 divides mainly between 1,000-gene clusters 3 and 10; cluster 6 divides mainly between 6 and 9. Those are specific boundaries to review with markers. A larger set can recover weaker populations, but it can also include unwanted programs. The adjusted Rand index (ARI) measures agreement between partitions without requiring their cluster numbers to match. It does not identify which feature set is more biologically useful. Compare marker specificity and known biology instead of choosing the layout that appears most separated.

Install an externally chosen feature set¶

If your study calls for a specific gene program or an externally selected set, supply those features directly. A graph made from a small marker panel will emphasize the identities that panel was chosen to distinguish; it cannot provide an independent test of those identities. features.snapshot accepts either a Boolean mask aligned to the complete feature metadata order or integer indexes into that feature table. The indexes refer to genes, not cells. It records an immutable selection artifact, so downstream artifacts can trace which genes were used from the artifact itself. Exactly one input form is required; duplicate or out-of-range indexes, misaligned masks, and empty selections are rejected. Verify the mask length and selected count before building a graph:

CyteArc’s metadata helpers can construct these selections without a separate lookup routine. sift and multi_sift return Boolean masks for mask=; get_index_by returns integer feature-table indexes for feature_indexes=.

Both producer calls return an ArtifactRef. Pass that reference directly when continuing a branch:

# Keep the normalization fixed while comparing reductions.
normalized = ds.features.normalize(cell_selection, custom_features)

Retain or persist exact refs in the analysis record. To request the complete feature universe, use the canonical all-features producer:

# Use the complete feature universe for this assay.
all_features = ds.features.universe(from_assay="RNA")

all_features is an immutable all-true artifact for this exact assay axis.

In your own workflow, the standard pipeline uses 1,000 HVGs by default. Change the count when there is a reason to include more or fewer genes:

# Run the pipeline with the explicitly changed feature setting.
ds.pipeline.run(hvg_count=2000)

Other HVG settings are available through the pipeline’s params mapping. For example, if a study needs genes that the default blacklist removes:

# Run the pipeline with the explicitly changed feature setting.
ds.pipeline.run(params={"hvg": {"blacklist": ""}})

Use the explicit stage methods when you want to compare feature sets on the same frozen cells, as we did above. For a simple detection threshold instead of an HVG model, use features.detected(cell_selection, min_cells=...).

For scATAC-seq, prevalent peak selection is the analogous step. features.prevalent_peaks returns a feature-selection artifact whose scientific identity contains only its feature_summary input and top_n parameter. See Interpret cell populations with scATAC-seq; peak prevalence, not an RNA mean-variance model, defines the candidate accessibility features.

Common mistakes and limitations¶

• Assuming the default blacklist is a functional classification: it is a set of name patterns, so ^CCN also removes CCND1 and ^HIST misses names like H2AC1, and supplying blacklist= replaces the whole default rather than adding to it. Check which names your annotation actually drops, and clear or narrow the pattern when your question needs mitochondrial, ribosomal, or sex-linked programs.

• Treating the feature-set size as a neutral parameter: top_n changes the PCA basis and the neighbourhood graph, so clusters can split or merge when it changes. The ARI measures agreement between partitions, not biological usefulness; review changed boundaries with markers and known biology before preferring the more separated layout.

• Using an external marker panel as an independent test: a graph built from the panel emphasizes exactly the identities the panel was chosen to distinguish, so it cannot validate them. Keep a genome-wide HVG selection for that check.