scanpy.experimental.pp.recipe_pearson_residuals#

scanpy.experimental.pp.recipe_pearson_residuals(adata, *, theta=100, clip=None, n_top_genes=1000, batch=None, chunksize=1000, n_comps=50, rng=None, kwargs_pca=mappingproxy({}), check_values=True, use=None, out=False, inplace=True, batch_key=None)[source]#

Full pipeline for HVG selection and normalization by analytic Pearson residuals [Lause et al., 2021].

Applies gene selection based on Pearson residuals. On the resulting subset, Pearson residual normalization and PCA are performed.

Expects raw count input.

Parameters:
adata AnnData

The annotated data matrix of shape n_obs × n_vars. Rows correspond to cells and columns to genes.

theta float (default: 100)

The negative binomial overdispersion parameter theta for Pearson residuals. Higher values correspond to less overdispersion (var = mean + mean^2/theta), and theta=np.inf corresponds to a Poisson model.

clip float | None (default: None)

Determines if and how residuals are clipped:

  • If None, residuals are clipped to the interval [-sqrt(n_obs), sqrt(n_obs)], where n_obs is the number of cells in the dataset (default behavior).

  • If any scalar c, residuals are clipped to the interval [-c, c]. Set clip=np.inf for no clipping.

n_top_genes int (default: 1000)

Number of highly-variable genes to keep. Mandatory if flavor='seurat_v3' or flavor='pearson_residuals'.

batch AdRef | str | None (default: None)

If specified (e.g. A.obs['batch']), highly-variable genes are selected within each batch separately and merged. This simple process avoids the selection of batch-specific genes and acts as a lightweight batch correction method. Genes are first sorted by how many batches they are a HVG. If flavor='pearson_residuals', ties are broken by the median rank (across batches) based on within-batch residual variance.

chunksize int (default: 1000)

If flavor='pearson_residuals', this dertermines how many genes are processed at once while computing the residual variance. Choosing a smaller value will reduce the required memory.

n_comps int | None (default: 50)

Number of principal components to compute in the PCA step.

rng int | integer | Sequence[int] | SeedSequence | Generator | BitGenerator | None (default: None)

Random number generator for setting the initial states for the optimization in the PCA step.

kwargs_pca Mapping[str, Any] (default: mappingproxy({}))

Dictionary of further keyword arguments passed on to scanpy.pp.pca().

check_values bool (default: True)

If True, checks if counts in selected layer are integers as expected by this function, and return a warning if non-integers are found. Otherwise, proceed without checking. Setting this to False can speed up code for large datasets.

use LayerAcc | str | None (default: None)

Which count matrix to use as input. Given by an accessor, e.g. A.X or A.layers['counts']. strs are anndata.acc.AdAcc.resolve()d, e.g. 'layers.counts'.

out bool (default: False)

Where to write the gene selection and Pearson residual-based PCA results. If True, the gene selection and Pearson residual-based PCA results is returned instead of written. If False, it is written to adata (see below).

inplace bool (default: True)

Deprecated since version 1.13.0: Use out=True instead of inplace=False.

inplace=False is out=True.

batch_key str | None (default: None)

Deprecated since version 1.13.0: Use batch instead, e.g. batch=A.obs['b'] for batch_key='b'.

The name of a .obs column, i.e. batch_key='b' is batch=A.obs['b'].

Return type:

tuple[AnnData, DataFrame] | None

Returns:

If out=True, separately returns the Pearson residual-based PCA results (as AnnData) and the gene selection results (as DataFrame). Otherwise, updates adata with the following fields for gene selection results:

.var['highly_variable']bool

boolean indicator of highly-variable genes.

.var['means']float

means per gene.

.var['variances']float

variances per gene.

.var['residual_variances']float

Pearson residual variance per gene. Averaged in the case of multiple batches.

.var['highly_variable_rank']float

Rank of the gene according to residual variance, median rank in the case of multiple batches.

.var['highly_variable_nbatches']int

If batch is given, this denotes in how many batches genes are detected as HVG.

.var['highly_variable_intersection']bool

If batch is given, this denotes the genes that are highly variable in all batches.

The following fields contain Pearson residual-based PCA results and normalization settings:

.uns['pearson_residuals_normalization']['pearson_residuals_df']

The subset of highly variable genes, normalized by Pearson residuals.

.uns['pearson_residuals_normalization']['theta']

The used value of the overdisperion parameter theta.

.uns['pearson_residuals_normalization']['clip']

The used value of the clipping parameter.

By default, the PCA results are stored in these fields (see out of pca(), which can be passed via kwargs_pca):

.obsm['X_pca']

PCA representation of data after gene selection and Pearson residual normalization.

.varm['PCs']

The principal components containing the loadings. This will contain empty rows for the genes not selected during HVG selection.

.uns['pca']['variance_ratio']

Ratio of explained variance.

.uns['pca']['variance']

Explained variance, equivalent to the eigenvalues of the covariance matrix.