Normalization and gene selection by analytical Pearson residuals

docs/api.rst

@@ -369,6 +369,24 @@ Collections of useful measurements for evaluating results.
    metrics.morans_i
 
 
+Experimental
+------------
+
+.. module:: scanpy.experimental
+.. currentmodule:: scanpy
+
+New methods that are in early development which are not (yet)
+integrated in Scanpy core.
+
+.. autosummary::
+   :toctree: generated/
+
+   experimental.pp.normalize_pearson_residuals
+   experimental.pp.normalize_pearson_residuals_pca
+   experimental.pp.highly_variable_genes
+   experimental.pp.recipe_pearson_residuals
+
+
 Classes
 -------
 

docs/references.rst

@@ -119,6 +119,10 @@ References
    *Laplacian Dynamics and Multiscale Modular Structure in Networks*
    `arXiv <https://arxiv.org/abs/0812.1770>`__.
 
+.. [Lause21] Lause *et al.* (2021)
+   *Analytic Pearson residuals for normalization of single-cell RNA-seq UMI data*,
+   `Genome Biology <https://doi.org/10.1186/s13059-021-02451-7>`__.
+
 .. [Leek12] Leek *et al.* (2012),
    *sva: Surrogate Variable Analysis. R package*
    `Bioconductor <https://doi.org/10.18129/B9.bioc.sva>`__.

docs/release-notes/1.9.0.rst

@@ -6,3 +6,13 @@
 - :func:`~scanpy.tl.filter_rank_genes_groups` now allows to filter with absolute values of log fold change :pr:`1649` :smaller:`S Rybakov`
 - :func:`~scanpy.pl.embedding_density` now allows more than 10 groups :pr:`1936` :smaller:`A Wolf`
 - :func:`~scanpy.logging.print_versions` now uses `session_info` :pr:`2089` :smaller:`P Angerer` :smaller:`I Virshup`
+
+.. rubric:: Experimental module
+
+- Added :mod:`scanpy.experimental` module!
+
+    - Added :func:`scanpy.experimental.pp.normalize_pearson_residuals` for Pearson Residuals normalization :pr:`1715` :smaller:`J Lause, G Palla, I Virshup`
+    - Added :func:`scanpy.experimental.pp.normalize_pearson_residuals_pca` for Pearson Residuals normalization and PCA :pr:`1715` :smaller:`J Lause, G Palla, I Virshup`
+    - Added :func:`scanpy.experimental.pp.highly_variable_genes` for HVG selection with Pearson Residuals :pr:`1715` :smaller:`J Lause, G Palla, I Virshup`
+    - Added :func:`scanpy.experimental.pp.normalize_pearson_residuals_pca` for Pearson Residuals normalization and dimensionality reduction with PCA :pr:`1715` :smaller:`J Lause, G Palla, I Virshup`
+    - Added :func:`scanpy.experimental.pp.recipe_pearson_residuals` for Pearson Residuals normalization, HVG selection and dimensionality reduction with PCA  :pr:`1715` :smaller:`J Lause, G Palla, I Virshup`

docs/tutorials.rst

@@ -97,11 +97,15 @@ See the `cell cycle`_ notebook.
 
 .. _cell cycle: https://nbviewer.jupyter.org/github/theislab/scanpy_usage/blob/master/180209_cell_cycle/cell_cycle.ipynb
 
+Normalization with Pearson Residuals
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Normalization of scRNA-seq data with Pearson Residuals, from [Lause21]_: :tutorial:`tutorial_pearson_residuals`
+
 .. image:: _static/img/tutorials/170522_visualizing_one_million_cells/tsne_1.3M.png
    :width: 120px
    :align: right
 
-
 Scaling Computations
 ~~~~~~~~~~~~~~~~~~~~
 

scanpy/__init__.py

@@ -14,7 +14,7 @@
     from . import tools as tl
     from . import preprocessing as pp
     from . import plotting as pl
-    from . import datasets, logging, queries, external, get, metrics
+    from . import datasets, logging, queries, external, get, metrics, experimental
 
     from anndata import AnnData, concat
     from anndata import (

scanpy/experimental/__init__.py

@@ -0,0 +1 @@
+from . import pp

scanpy/experimental/_docs.py

@@ -0,0 +1,79 @@
+"""Shared docstrings for experimental function parameters.
+"""
+
+doc_adata = """\
+adata
+    The annotated data matrix of shape `n_obs` × `n_vars`.
+    Rows correspond to cells and columns to genes.
+"""
+
+doc_dist_params = """\
+theta
+    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
+    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.
+"""
+
+doc_check_values = """\
+check_values
+    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.
+"""
+
+doc_layer = """\
+layer
+    Layer to use as input instead of `X`. If `None`, `X` is used.
+"""
+
+doc_subset = """\
+subset
+    Inplace subset to highly-variable genes if `True` otherwise merely indicate
+    highly variable genes.
+"""
+
+doc_genes_batch_chunk = """\
+n_top_genes
+    Number of highly-variable genes to keep. Mandatory if `flavor='seurat_v3'` or
+    `flavor='pearson_residuals'`.
+batch_key
+    If specified, 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
+    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.
+"""
+
+doc_pca_chunk = """\
+n_comps
+    Number of principal components to compute in the PCA step.
+random_state
+    Random seed for setting the initial states for the optimization in the PCA step.
+kwargs_pca
+    Dictionary of further keyword arguments passed on to `scanpy.pp.pca()`.
+"""
+
+doc_inplace = """\
+inplace
+    If `True`, update `adata` with results. Otherwise, return results. See below for
+    details of what is returned.
+"""
+
+doc_copy = """\
+copy
+    If `True`, the function runs on a copy of the input object and returns the
+    modified copy. Otherwise, the input object is modified direcly. Not compatible
+    with `inplace=False`.
+"""

scanpy/experimental/pp/__init__.py

@@ -0,0 +1,8 @@
+from scanpy.experimental.pp._normalization import (
+    normalize_pearson_residuals,
+    normalize_pearson_residuals_pca,
+)
+
+from scanpy.experimental.pp._highly_variable_genes import highly_variable_genes
+
+from scanpy.experimental.pp._recipes import recipe_pearson_residuals