PyPsupertime

PyPsupertime is a scalable python re-implementation of the R package psupertime for analysis of single-cell RNA sequencing data where the cell have an ordinal annotation (e.g. time series, or dosage). It can be used to identify a small subset of cells which contribute to the ordering and reconstruct a pseudotime.

The original methodology is published in Bioinformatics: https://doi.org/10.1093/bioinformatics/btac227

Getting Started

Install via pip:

pip install pypsupertime

This installs pypsupertime and its dependencies automatically

We recommend installing inside a virtualenv or pipenv environment.

Description

This package implements a modular API for preprocessing of single cell data, restructuring of the input data under different statistical assumptions, as well as creating and fitting memory-efficient supervised ordinal logistic models.

The central idea of the original work remained unchanged: Find a linear model that accurately predicts the ordinal labels while being restricted to a sparse set of features (genes), by searching along a path of regularization hyperparameters. To adress memory inefficiencies in the original work, the coordinate descent approach is replaced in favor of a stochistic gradient descent model with online fitting. Additionaly, new parametrized penalties are possible, allowing more control over the sparsity.

From a statistical perspective, the ordinal nature of the input data can be modeled under the cumulative proportional odds, forward continuation ratio and backward continuation ratio assumption.

All model and preprocessing classes fit implement scikit-learn estimators or transformers to fit seamlessly into its ecosystem and are wrapped by the Psupertime class at the core of this pacakge that allows all input and output data to be represented as anndata objects.

Please find a more detailed description documentation hosted by readthedocs here.

Basic Usage

The code below runs a psupertime analysis with default settings on a data set represented as AnnData object and stored in a .h5ad file. The data has a numeric ordinal cell annotation representing stored in the obs dataframe under the key "time".

from pypsupertime import Psupertime
p = Psupertime()
anndata = p.run("/path/to/data_sce.h5ad", "time")

Input Data: n_genes=24153, n_cells=992
Preprocessing: done. mode='all', n_genes=11305, n_cells=992
Grid Search CV: CPUs=4, n_folds=5
Regularization: done   
Refit on all data: done. accuracy=0.5195, n_genes=113
Total elapsed time:  0:01:26.141356

The code loads the single-cell data and perfoms default preprocessing, then runs a 5-fold cross validatied grid search along the default regularization path to identify the hyperparameter which results in the best-scoring sparse model, and finally refits the model with the ideal regularization on all data. Using that model, relevant genes are identified and the psupertime is predicted for all cells.

The following snippets show how to quickly inspect and evaluate the results. The regularization progress, model performance, and selected genes are shown here, but many more are available.

p.plot_grid_search(title="Grid Search")

p.plot_model_perf((adata.X, adata.obs.time), figsize=(6,5))

p.plot_identified_gene_coefficients(adata, n_top=20)

p.plot_labels_over_psupertime(adata, "time")

It is highly recommended to be aware of the preprocessing steps performed or perform key preprocessing manually. For a complete overview, look at the documentation.

Development Roadmap

• Hihgly variable gene selection
• (Re-)Add early stopping
• Add "layer"/"rep" parameter to Psupertime for specifying data in adata objects
• Extension of the pypsupertime.plots module with further analyses
• Extension of the Preprocessing to allow custom pipelines (see version 1.1.0)
• Integration into the scanpy project
• Unit Tests, when the code is stable enough

Changelog:

• Version 2.2.2:

  • Fixes learning rate calculation with "optimal" heuristic:
    • The initial learning rate had inverted sign, causing a step in the direction away from the gradient
    • The number of samples was added to early, essentially skipping the first iteration.
    • The learning rate did not decay as expected due to the sign error: For some sample sizes, the learning rate could increase for one epoch and change sign before decreasing again.
  • Changed from mean to sum as loss reduction method in CumulativePenaltySGD for checking early stopping condition
  • Changes grid search class to retrieve fitted model with tracked scores, instead of instantiating a new estimator.
  • Remove all but "optimal" learning rates for CumulativePenaltySGDModel (again)
  • Introduces min_iter parameter to train for a minimum number of epochs

• Version 2.2.1:

  • Bugfix: Attribute eta0 was never set in base class PsupertimeBaseModel
  • Fix optimal learning rate heuristic:
    • Models did not converge using the "optimal" learning rate method and were underfitted

• Version 2.2.0:

  • Adds adaptive and constant learning rates
  • Adds eta0 parameter to base class PsupertimeBaseModel
  • Breaks backward compatibility of learning_rate parameter

• Version 2.1.16:

  • Fixes critical training loop bug introduced after 2.1.9

• Version 2.1.15:

  • Changes early stopping behavior:
    • Always check if losses still improve based on training loss
    • if early_stopping is set to True, use a dedicated validation split instead
  • Verbosity levels 3 and 4 added: Output loss information

• Version 2.1.14:

  • Fixes dof trackig in ThresholdSGDModel

• Version 2.1.13:

  • Adds gamma hyperparameter for exponential learning rate decay to ThresholdSGDModel

• Version 2.1.12:

  • Fixes bug affecting score tracking in SGDModel

• Version 2.1.11:

  • Force track scores when refitting

• Version 2.1.10:

  • Fixes early stopping check

• Version 2.1.9:

  • Adds Pytorch dependency for new models to buildfile

• Version 2.1.8:

  • Bugfixes in score tracking
  • Adds parameter random_state to psupertime directly (had to be passed as estimator_param before)

• Version 2.1.7:

  • Refactoring of model classes and PsupertimeBaseModel to remove redundancies and improve readability
  • Bugfixes
  • Changes defaults
  • Tracks scores when refitting

• Version 2.0.0:

  • Changes undelying binary model from sklearn.SGDClassifier to pytorch model to gain more control over the loss and penalty calculation. The thresholds can now be trained without penalty. These changes remove/change some of the previous parameters and attributes in model classes

• Version 1.2.2:

  • Bugfixes
  • Removes inplace option in run(). It cannot currently be enforced in preprocessing. A copy of the input adata will always be created and the processed object returned.

• Version 1.2.0:

  • Returns an anndata object in Psupertime.run() if a filename is given

• Versions 1.1.1, 1.1.2, 1.1.3:

  • Bugfixes

• Version 1.1.0:

  • Add preprocessing_class parameter to enable using custom / no preprocessing
  • Adds heuristic for selecting the lowest regularization parameter when none is specified
  • Adds shorthand for selecting optimal regularization parameter at 1/2 standard error from the best score
  • Fix bug when using smooth with sparse matrices