Select Batch-Associated Features

Identifies genes that are associated with batch/cohort effects but NOT with the biological outcome, using PCA plus per-component ANOVA against supplied cohort labels. The returned genes act as proxies for the batch axis: they are useful as a confounding diagnostic in their own right, and as log-ratio denominators, where log(Signal) - log(BatchProxy) cancels batch-specific scaling factors (see construct_dual_stream_ratios()).

select_batch_associated_features(
  x,
  cohort,
  y = NULL,
  n_features = 20L,
  n_pcs = 50L,
  batch_pvalue_threshold = 0.1,
  biology_pvalue_threshold = 0.2,
  biology_penalty_quantile = 0.5
)

Arguments

x

Expression matrix (samples x genes). Should contain log-transformed expression values.

cohort

Factor of cohort/batch labels, one per row of x.

y

Binary outcome labels for biology filtering. Set to NULL to skip biology filtering (features selected purely on batch association).

n_features

Number of batch-associated genes to select (default 20).

n_pcs

Number of principal components to analyze (default 50).

batch_pvalue_threshold

P-value threshold for batch association. PCs with ANOVA p < this threshold are considered batch-associated (default 0.10).

biology_pvalue_threshold

P-value threshold for biology non-association. PCs with biology p > this threshold are considered biology-free (default 0.20).

biology_penalty_quantile

Quantile of the per-gene biology scores used to stabilise the batch/biology ratio (default 0.5, the median). This acts as a data-driven estimate of the biology-score noise floor; see Details. Larger values weaken the biology penalty, with 1.0 reducing the ranking to batch_score alone.

Value

A list with components:

features

Character vector of selected batch-associated gene names.

biology_epsilon

The stabilising constant actually added to the biology scores when computing ratio.

diagnostics

Data frame with per-PC statistics including variance explained, batch p-value/R^2, and biology p-value/R^2.

problematic_cohort

The cohort most associated with batch effects.

one_vs_rest

Data frame with per-cohort batch separation metrics.

pure_batch_pcs

Integer vector of PC indices identified as pure batch.

gene_scores

Data frame with per-gene batch and biology scores.

pca

The prcomp object from PCA on input data.

Details

The algorithm proceeds as follows:

1. PCA on combined expression matrix

2. Test each PC for batch association (ANOVA) and biology association (if y provided)

3. Identify "pure batch" PCs: batch-significant AND biology-nonsignificant

4. If no pure batch PCs found, fall back to LDA-based cohort discrimination

5. Extract genes with highest loadings on pure batch PCs

6. Score genes by ratio of batch to biology contribution

7. Return top n_features genes

Step 6 ranks by batch_score / (biology_score + epsilon). Both scores are variance-weighted means of absolute PCA loadings, so they share a scale that depends on the number of genes and the data at hand. epsilon is therefore set to the biology_penalty_quantile quantile of the observed biology scores rather than a fixed constant: a fixed tiny epsilon lets genes whose biology score is pure numerical noise (batch and biology both near zero) generate arbitrarily large ratios and outrank genuinely batch-driven genes. The median default assumes most genes are biologically null, making it a robust estimate of the biology-score noise floor.

Unlike select_ruv_features(), which estimates a latent unwanted-variation subspace from negative control genes, this function requires explicit cohort/batch labels and targets the batch axis directly.

See also

construct_dual_stream_ratios() for using the selected genes as ratio denominators, select_ruv_features() for the label-free alternative.

Examples

if (FALSE) { # \dontrun{
# Assume x is samples x genes, cohort and y are factors
result <- select_batch_associated_features(
  x = expression_matrix,
  cohort = cohort_labels,
  y = disease_labels,
  n_features = 20
)

# Use the batch-proxy genes as denominators
ratios <- construct_dual_stream_ratios(
  x = expression_matrix,
  numerators = signal_genes,
  denominators = result$features
)
} # }