Repository: theislab/scib (language: python)

Issue:
Use pre-commit
### Start using pre-commit

- [x] include working pre-commit config file
- [x] update syntax of all python files
- [x] include dependencies to 'dev' section of setup.cfg
- [x] include README section mentioning the installatio of scib with dev dependencies for developers

### Suggestions from @michalk8 :

> - [x] remove 3.7 CI, add 3.9 (unless you need to support it)
> - [x] use pre-commits (add e.g. black) instead of running raw flake8 (this PR was)
> - [x] add tests coverage
> - [x] not sure how much work, but macOS test CI would be nice to have

with emphasis on the first 2 points

### Notes:

+ This PR targets pre-commit integration and syntax updates. Separate points in the list above should be addressed in separate PRs.
+ We still need to support python 3.7, as this is the version used in scib-pipeline (expect users working with it already)


Discussion / hints:
I'd also recommend purging large files (most likely .ipynb) from the git repo history. When I cloned this repo, I downloaded like 300MiB of data, see also https://rtyley.github.io/bfg-repo-cleaner/.

Fix patch:
diff --git a/.codecov.yaml b/.codecov.yaml
new file mode 100644
index 00000000..829e56c0
--- /dev/null
+++ b/.codecov.yaml
@@ -0,0 +1,17 @@
+# Based on pydata/xarray
+codecov:
+    require_ci_to_pass: no
+
+coverage:
+    status:
+        project:
+            default:
+                # Require 1% coverage, i.e., always succeed
+                target: 1
+        patch: false
+        changes: false
+
+comment:
+    layout: diff, flags, files
+    behavior: once
+    require_base: no
diff --git a/.editorconfig b/.editorconfig
new file mode 100644
index 00000000..2fe0ce08
--- /dev/null
+++ b/.editorconfig
@@ -0,0 +1,12 @@
+root = true
+
+[*]
+indent_style = space
+indent_size = 4
+end_of_line = lf
+charset = utf-8
+trim_trailing_whitespace = true
+insert_final_newline = true
+
+[Makefile]
+indent_style = tab
diff --git a/.github/workflows/deployment.yml b/.github/workflows/deployment.yml
new file mode 100644
index 00000000..40dcd0da
--- /dev/null
+++ b/.github/workflows/deployment.yml
@@ -0,0 +1,38 @@
+name: Deployment
+
+on:
+    push:
+        branches: [main]
+        tags: [v*]
+
+jobs:
+    build:
+        if: ${{ github.event_name == 'push' && startsWith(github.ref, 'refs/tags') }}
+        runs-on: ubuntu-latest
+
+        steps:
+        -   uses: actions/checkout@v2
+            with:
+                fetch-depth: 0
+
+        -   name: Set up Python 3.9
+            uses: actions/setup-python@v2
+            with:
+                python-version: 3.9
+        -   name: Install pip dependencies
+            run: |
+                python -m pip install --upgrade pip
+                pip install build twine
+
+        -   name: Build a binary wheel and a source tarball
+            run: |
+                python -m build --sdist --wheel --outdir dist/
+                twine check dist/*
+
+        -   name: Publish on PyPI
+            uses: pypa/gh-action-pypi-publish@master
+            with:
+                user: __token__
+                password: ${{ secrets.PYPI_TOKEN }}
+                skip_existing: true
+                verbose: true
diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
deleted file mode 100644
index 43482b74..00000000
--- a/.github/workflows/python-package.yml
+++ /dev/null
@@ -1,42 +0,0 @@
-# This workflow will install Python dependencies, run tests and lint with a variety of Python versions
-# For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions
-
-name: scIB package via pip
-
-on: [push]
-
-jobs:
-  build:
-
-    runs-on: ubuntu-latest
-    strategy:
-      matrix:
-        python-version: ['3.7', '3.8']
-
-    steps:
-    - uses: actions/checkout@v2
-    - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v2
-      with:
-        python-version: ${{ matrix.python-version }}
-    - name: Install dependencies
-      run: |
-        python -m pip install --upgrade pip
-        python -m pip install flake8 pytest
-    - name: Lint with flake8
-      run: |
-        # stop the build if there are Python syntax errors or undefined names
-        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
-        # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
-        flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
-    - name: Install package
-      run: pip install .[test]
-    - name: Import package
-      run: |
-        pip list
-        cd tests
-        python -c 'import scib'
-    - name: Test with pytest
-      run: |
-        pytest --durations 0 -s
-
diff --git a/.gitignore b/.gitignore
index 053dabda..b5ec7c18 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,7 @@
-build/
 scratch/
 
+htmlcov/
+.coverage
 **.h5ad
 *test_output*
 
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
new file mode 100644
index 00000000..ad46b179
--- /dev/null
+++ b/.pre-commit-config.yaml
@@ -0,0 +1,76 @@
+fail_fast: false
+default_language_version:
+    python: python3
+default_stages:
+-   commit
+-   push
+minimum_pre_commit_version: 2.16.0
+repos:
+-   repo: https://github.com/psf/black
+    rev: 22.3.0
+    hooks:
+    -   id: black
+-   repo: https://github.com/asottile/blacken-docs
+    rev: v1.12.1
+    hooks:
+    -   id: blacken-docs
+-   repo: https://github.com/macisamuele/language-formatters-pre-commit-hooks
+    rev: v2.3.0
+    hooks:
+    -   id: pretty-format-yaml
+        args: [--autofix, --indent, '4']
+-   repo: https://github.com/pre-commit/mirrors-autopep8
+    rev: v1.6.0
+    hooks:
+    -   id: autopep8
+        args: [-i]
+-   repo: https://github.com/PyCQA/isort
+    rev: 5.10.1
+    hooks:
+    -   id: isort
+-   repo: https://github.com/asottile/yesqa
+    rev: v1.3.0
+    hooks:
+    -   id: yesqa
+        additional_dependencies:
+        -   flake8-tidy-imports
+        -   flake8-docstrings
+        -   flake8-rst-docstrings
+        -   flake8-comprehensions
+        -   flake8-bugbear
+        -   flake8-blind-except
+-   repo: https://github.com/pre-commit/pre-commit-hooks
+    rev: v4.1.0
+    hooks:
+    -   id: detect-private-key
+    -   id: check-ast
+    -   id: end-of-file-fixer
+    -   id: mixed-line-ending
+        args: [--fix=lf]
+    -   id: trailing-whitespace
+    -   id: check-case-conflict
+-   repo: https://github.com/myint/autoflake
+    rev: v1.4
+    hooks:
+    -   id: autoflake
+        args:
+        -   --in-place
+        -   --remove-all-unused-imports
+        -   --remove-unused-variable
+        -   --ignore-init-module-imports
+-   repo: https://github.com/PyCQA/flake8
+    rev: 4.0.1
+    hooks:
+    -   id: flake8
+        additional_dependencies:
+        -   flake8-tidy-imports
+        -   flake8-docstrings
+        -   flake8-rst-docstrings
+        -   flake8-comprehensions
+        -   flake8-bugbear
+        -   flake8-blind-except
+-   repo: https://github.com/asottile/pyupgrade
+    rev: v2.31.1
+    hooks:
+    -   id: pyupgrade
+        args: [--py3-plus, --py37-plus, --keep-runtime-typing]
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 9aae0277..02dc5666 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -7,16 +7,16 @@ version: 2
 
 # Set the version of Python and other tools you might need
 build:
-  os: ubuntu-20.04
-  tools:
-    python: "3.8"
-  apt_packages:
-    - r-base
+    os: ubuntu-20.04
+    tools:
+        python: '3.8'
+    apt_packages:
+    -   r-base
 
 # Build documentation in the docs/ directory with Sphinx
 sphinx:
-   configuration: docs/source/conf.py
-   fail_on_warning: false
+    configuration: docs/source/conf.py
+    fail_on_warning: false
 
 # If using Sphinx, optionally build your docs in additional formats such as PDF
 # formats:
@@ -24,8 +24,8 @@ sphinx:
 
 # Optionally declare the Python requirements required to build your docs
 python:
-  install:
-    - method: pip
-      path: .
-      extra_requirements:
-        - docs
+    install:
+    -   method: pip
+        path: .
+        extra_requirements:
+        -   docs
diff --git a/LICENSE.txt b/LICENSE.txt
index 3aa54c49..9cfe1282 100644
--- a/LICENSE.txt
+++ b/LICENSE.txt
@@ -1,5 +1,5 @@
 MIT License
-Copyright (c) 2018 
+Copyright (c) 2018
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
 in the Software without restriction, including without limitation the rights
@@ -14,4 +14,4 @@ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-SOFTWARE.
\ No newline at end of file
+SOFTWARE.
diff --git a/README.md b/README.md
index e766747d..99474c9a 100644
--- a/README.md
+++ b/README.md
@@ -1,19 +1,29 @@
+[![Stars](https://img.shields.io/github/stars/theislab/scib?logo=GitHub&color=yellow)](https://github.com/theislab/scib/stargazers)
+[![PyPI](https://img.shields.io/pypi/v/scib?logo=PyPI)](https://pypi.org/project/scib)
+[![PyPIDownloads](https://pepy.tech/badge/scib)](https://pepy.tech/project/scib)
+[![Build Status](https://github.com/theislab/scib/actions/workflows/test.yml/badge.svg)](https://github.com/theislab/scib/actions/workflows/test.yml)
+[![Documentation](https://readthedocs.org/projects/scib/badge/?version=latest)](https://scib.readthedocs.io/en/latest/?badge=latest)
+[![codecov](https://codecov.io/gh/theislab/scib/branch/main/graph/badge.svg?token=M1nuTpAxyS)](https://codecov.io/gh/theislab/scib)
+[![pre-commit](https://img.shields.io/badge/pre--commit-enabled-brightgreen?logo=pre-commit&logoColor=white)](https://github.com/pre-commit/pre-commit)
+
 # Benchmarking atlas-level data integration in single-cell genomics
 
 This repository contains the code for the `scib` package used in our benchmarking study for data integration tools.
 In [our study](https://doi.org/10.1038/s41592-021-01336-8), we benchmark 16 methods (see Tools) with 4 combinations of
-preprocessing steps leading to 68 methods combinations on 85 batches of gene expression and chromatin accessibility data.
+preprocessing steps leading to 68 methods combinations on 85 batches of gene expression and chromatin accessibility
+data.
 
 ![Workflow](https://raw.githubusercontent.com/theislab/scib/main/figure.png)
 
 ## Resources
 
-+ The git repository of the [`scib` package](https://github.com/theislab/scib) and its [documentation](https://scib.readthedocs.io/).
-+ The reusable pipeline we used in the study can be found in the
+- The git repository of the [`scib` package](https://github.com/theislab/scib) and
+  its [documentation](https://scib.readthedocs.io/).
+- The reusable pipeline we used in the study can be found in the
   separate [scib pipeline](https://github.com/theislab/scib-pipeline.git) repository. It is reproducible and automates
   the computation of preprocesssing combinations, integration methods and benchmarking metrics.
-+ On our [website](https://theislab.github.io/scib-reproducibility) we visualise the results of the study.
-+ For reproducibility and visualisation we have a dedicated
+- On our [website](https://theislab.github.io/scib-reproducibility) we visualise the results of the study.
+- For reproducibility and visualisation we have a dedicated
   repository: [scib-reproducibility](https://github.com/theislab/scib-reproducibility).
 
 ### Please cite:
@@ -23,22 +33,19 @@ Nat Methods 19, 41–50 (2022). [https://doi.org/10.1038/s41592-021-01336-8](htt
 
 ## Package: scib
 
-We created the python package called `scib` that uses `scanpy` to streamline the integration of single-cell datasets
-and evaluate the results.
-The package contains several modules for preprocessing an ``anndata`` object, running integration methods and
-evaluating the resulting using a number of metrics.
-For preprocessing, ``scib.preprocessing`` (or ``scib.pp``) contains functions for normalising, scaling or batch-aware
-selection of highly variable genes.
-Functions for the integration methods are in ``scib.integration`` or for short ``scib.ig`` and metrics are under
-``scib.metrics`` (or ``scib.me``).
+We created the python package called `scib` that uses `scanpy` to streamline the integration of single-cell datasets and
+evaluate the results. The package contains several modules for preprocessing an `anndata` object, running integration
+methods and evaluating the resulting using a number of metrics. For preprocessing, `scib.preprocessing` (or `scib.pp`)
+contains functions for normalising, scaling or batch-aware selection of highly variable genes. Functions for the
+integration methods are in `scib.integration` or for short `scib.ig` and metrics are under
+`scib.metrics` (or `scib.me`).
 
 The `scib` python package is available on [PyPI](https://pypi.org/) and can be installed through
 
-```
+```shell
 pip install scib
 ```
 
-
 Import `scib` in python:
 
 ```python
@@ -106,4 +113,21 @@ Tools that are compared include:
 - [scVI](https://github.com/YosefLab/scVI) 0.6.7
 - [Seurat v3](https://github.com/satijalab/seurat) 3.2.0 CCA (default) and RPCA
 - [TrVae](https://github.com/theislab/trvae) 0.0.1
-- [TrVaep](https://github.com/theislab/trvaep) 0.1.0
\ No newline at end of file
+- [TrVaep](https://github.com/theislab/trvaep) 0.1.0
+
+## Development
+
+For developing this package, please make sure to install additional dependencies so that you can use `pytest` and
+`pre-commit`.
+
+```shell
+pip install -e '.[test,dev]'
+```
+
+Please refer to the `setup.cfg` for more optional dependencies.
+
+Install `pre-commit` to the repository for running it automatically every time you commit in git.
+
+```shell
+pre-commit install
+```
diff --git a/docs/source/.gitignore b/docs/source/.gitignore
index b8798c54..8c6492ca 100644
--- a/docs/source/.gitignore
+++ b/docs/source/.gitignore
@@ -1,1 +1,1 @@
-api/
\ No newline at end of file
+api/
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 2ab4320a..b9770b37 100755
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -12,16 +12,19 @@
 #
 import os
 import sys
+
 import scib
 
-sys.path.insert(0, os.path.abspath('../..'))
+sys.path.insert(0, os.path.abspath("../.."))
 
 # -- Project information -----------------------------------------------------
 
-project = 'scib'
-copyright = '2021, Malte D. Luecken, Maren Buettner, Daniel C. Strobl, Michaela F. Mueller'
-author = 'Malte D. Luecken, Maren Buettner, Daniel C. Strobl, Michaela F. Mueller'
-github_url = 'https://github.com/theislab/scib'
+project = "scib"
+copyright = (
+    "2021, Malte D. Luecken, Maren Buettner, Daniel C. Strobl, Michaela F. Mueller"
+)
+author = "Malte D. Luecken, Maren Buettner, Daniel C. Strobl, Michaela F. Mueller"
+github_url = "https://github.com/theislab/scib"
 
 # The full version, including alpha/beta/rc tags
 release = scib.__version__
@@ -32,14 +35,14 @@
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
 extensions = [
-    'sphinx.ext.duration',
-    'sphinx.ext.doctest',
-    'sphinx.ext.autodoc',
-    'sphinx.ext.autosummary',
-    'sphinx.ext.intersphinx',
-    'sphinx_automodapi.automodapi',
-    'sphinx_automodapi.smart_resolver',
-    'myst_parser'
+    "sphinx.ext.duration",
+    "sphinx.ext.doctest",
+    "sphinx.ext.autodoc",
+    "sphinx.ext.autosummary",
+    "sphinx.ext.intersphinx",
+    "sphinx_automodapi.automodapi",
+    "sphinx_automodapi.smart_resolver",
+    "myst_parser",
 ]
 numpydoc_show_class_members = False
 
@@ -56,7 +59,7 @@
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
 #
-html_theme = 'sphinx_rtd_theme'
+html_theme = "sphinx_rtd_theme"
 
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
diff --git a/docs/source/installation.rst b/docs/source/installation.rst
index 3fa4cffd..f9dac71b 100644
--- a/docs/source/installation.rst
+++ b/docs/source/installation.rst
@@ -18,14 +18,14 @@ Installation with pip
 
 The ``scib`` python package is available on `PyPI <https://pypi.org/>`_ and can be installed through
 
-.. code-block:: python
+.. code-block:: bash
 
     pip install scib
 
 
 Alternatively, you can also install the package directly from GitHub directly via
 
-.. code-block:: python
+.. code-block:: bash
 
     pip install git+https://github.com/theislab/scib.git
 
@@ -62,7 +62,7 @@ installation of ``scib``. You can install optional Python dependencies via pip a
     Zsh often doesn't like square brackets. If you are a zsh user, use quotation marks around any statements containing
     square brackets. For example:
 
-    .. code-block:: python
+    .. code-block:: bash
 
         pip install 'scib[bbknn]'
 
diff --git a/docs/source/scib_integration.rst b/docs/source/scib_integration.rst
index 1631495e..567e1734 100644
--- a/docs/source/scib_integration.rst
+++ b/docs/source/scib_integration.rst
@@ -3,18 +3,18 @@ Integration
 
 Integration method functions require the preprocessed ``anndata`` object (here ``adata``) and the name of the batch column
 in ``adata.obs`` (here ``'batch'``).
-The methods can be called using the following, where ``<method>`` is the name of the integration method.
+The methods can be called using the following, where ``integration_method`` is the name of the integration method.
 
 .. code-block:: python
 
-    scib.ig.<method>(adata, batch='batch')
+    scib.ig.integration_method(adata, batch="batch")
 
 
 For example, in order to run Scanorama, on a dataset, call:
 
 .. code-block:: python
 
-    scib.ig.scanorama(adata, batch='batch')
+    scib.ig.scanorama(adata, batch="batch")
 
 .. warning::
 
@@ -22,7 +22,7 @@ For example, in order to run Scanorama, on a dataset, call:
 
     .. code-block:: python
 
-        scib.integration.run<method>(adata, batch='batch')
+        scib.integration.runIntegrationMethod(adata, batch="batch")
 
     Please use the snake_case naming without the ``run`` prefix.
 
@@ -32,8 +32,8 @@ For these, you need to additionally provide the corresponding label column of ``
 
 .. code-block:: python
 
-    scib.ig.scgen(adata, batch='batch', cell_type ='cell_type')
-    scib.ig.scanvi(adata, batch='batch', labels ='cell_type')
+    scib.ig.scgen(adata, batch="batch", cell_type="cell_type")
+    scib.ig.scanvi(adata, batch="batch", labels="cell_type")
 
 
 .. automodapi:: scib.integration
diff --git a/docs/source/scib_preprocessing.rst b/docs/source/scib_preprocessing.rst
index b990c1f8..eac5828c 100644
--- a/docs/source/scib_preprocessing.rst
+++ b/docs/source/scib_preprocessing.rst
@@ -19,4 +19,3 @@ The most relevant preprocessing steps are:
     :skip: readConos
     :skip: readSeurat
     :skip: saveSeurat
-
diff --git a/scib/__init__.py b/scib/__init__.py
index 83bac6b6..f2ce75ce 100644
--- a/scib/__init__.py
+++ b/scib/__init__.py
@@ -3,28 +3,27 @@
 except ImportError:  # for Python<3.8
     import importlib_metadata as metadata
 
-__version__ = metadata.version('scib')
+__version__ = metadata.version("scib")
 
-from . import integration, metrics, preprocessing
-from . import utils as utils
+from . import integration, metrics, preprocessing, utils
 from ._package_tools import rename_func
 from .metrics import clustering
 
 alias_func_map = {
-    'runScanorama': integration.scanorama,
-    'runTrVae': integration.trvae,
-    'runTrVaep': integration.trvaep,
-    'runScGen': integration.scgen,
-    'runScvi': integration.scvi,
-    'runScanvi': integration.scanvi,
-    'runMNN': integration.mnn,
-    'runBBKNN': integration.bbknn,
-    'runSaucie': integration.saucie,
-    'runCombat': integration.combat,
-    'runDESC': integration.desc,
-    'readConos': preprocessing.read_conos,
-    'readSeurat': preprocessing.read_seurat,
-    'saveSeurat': preprocessing.save_seurat,
+    "runScanorama": integration.scanorama,
+    "runTrVae": integration.trvae,
+    "runTrVaep": integration.trvaep,
+    "runScGen": integration.scgen,
+    "runScvi": integration.scvi,
+    "runScanvi": integration.scanvi,
+    "runMNN": integration.mnn,
+    "runBBKNN": integration.bbknn,
+    "runSaucie": integration.saucie,
+    "runCombat": integration.combat,
+    "runDESC": integration.desc,
+    "readConos": preprocessing.read_conos,
+    "readSeurat": preprocessing.read_seurat,
+    "saveSeurat": preprocessing.save_seurat,
 }
 
 for alias, func in alias_func_map.items():
diff --git a/scib/_package_tools.py b/scib/_package_tools.py
index 73c797dd..05e4509d 100644
--- a/scib/_package_tools.py
+++ b/scib/_package_tools.py
@@ -2,7 +2,7 @@
 import warnings
 from functools import wraps
 
-warnings.simplefilter('default')  # or 'always'
+warnings.simplefilter("default")  # or 'always'
 
 
 def wrap_func_naming(func, name):
@@ -15,7 +15,7 @@ def wrapper(*args, **kwargs):
         warnings.warn(
             f"Mixed case function naming is deprecated for '{name}'. "
             f"Please use '{func.__name__}' instead.",
-            DeprecationWarning
+            DeprecationWarning,
         )
         return func(*args, **kwargs)
 
diff --git a/scib/exceptions.py b/scib/exceptions.py
index 2f402441..e9a8e125 100644
--- a/scib/exceptions.py
+++ b/scib/exceptions.py
@@ -1,13 +1,13 @@
 class IntegrationMethodNotFound(ModuleNotFoundError):
-
     def __init__(self, exception):
-        self.message = f"\n{exception.name} is not installed by default, " \
-                       "please make sure you install it manually."
+        self.message = (
+            f"\n{exception.name} is not installed by default, "
+            "please make sure you install it manually."
+        )
         super().__init__(self.message)
 
 
 class RLibraryNotFound(ModuleNotFoundError):
-
     def __init__(self, exception):
         self.message = f"\nproblem loading library: {exception}"
         super().__init__(self.message)
diff --git a/scib/integration.py b/scib/integration.py
index d287a169..032dce86 100644
--- a/scib/integration.py
+++ b/scib/integration.py
@@ -12,11 +12,13 @@
 from . import utils
 from .exceptions import IntegrationMethodNotFound
 
-rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)  # Ignore R warning messages
+rpy2.rinterface_lib.callbacks.logger.setLevel(
+    logging.ERROR
+)  # Ignore R warning messages
 
 
 def scanorama(adata, batch, hvg=None, **kwargs):
-    """ Scanorama wrapper function
+    """Scanorama wrapper function
 
     Based on `scanorama <https://github.com/brianhie/scanorama>`_ version 1.7.0
 
@@ -37,7 +39,7 @@ def scanorama(adata, batch, hvg=None, **kwargs):
     corrected = anndata.AnnData.concatenate(
         *corrected, batch_key=batch, batch_categories=categories, index_unique=None
     )
-    corrected.obsm['X_emb'] = corrected.obsm['X_scanorama']
+    corrected.obsm["X_emb"] = corrected.obsm["X_scanorama"]
     # corrected.uns['emb']=True
 
     return corrected
@@ -62,18 +64,16 @@ def trvae(adata, batch, hvg=None):
     utils.check_sanity(adata, batch, hvg)
     n_batches = len(adata.obs[batch].cat.categories)
 
-    train_adata, valid_adata = trvae.utils.train_test_split(
-        adata,
-        train_frac=0.80
-    )
+    train_adata, valid_adata = trvae.utils.train_test_split(adata, train_frac=0.80)
 
     condition_encoder = trvae.utils.create_dictionary(
-        adata.obs[batch].cat.categories.tolist(), [])
+        adata.obs[batch].cat.categories.tolist(), []
+    )
 
     network = trvae.archs.trVAEMulti(
         x_dimension=train_adata.shape[1],
         n_conditions=n_batches,
-        output_activation='relu'
+        output_activation="relu",
     )
 
     network.train(
@@ -92,9 +92,9 @@ def trvae(adata, batch, hvg=None):
 
     network.get_corrected(adata, labels, return_z=False)
 
-    adata.obsm['X_emb'] = adata.obsm['mmd_latent']
-    del adata.obsm['mmd_latent']
-    adata.X = adata.obsm['reconstructed']
+    adata.obsm["X_emb"] = adata.obsm["mmd_latent"]
+    del adata.obsm["mmd_latent"]
+    adata.X = adata.obsm["reconstructed"]
 
     return adata
 
@@ -124,11 +124,17 @@ def trvaep(adata, batch, hvg=None):
     if issparse(adata.X):
         adata.X = adata.X.A
 
-    model = trvaep.CVAE(adata.n_vars, num_classes=n_batches,
-                        encoder_layer_sizes=[64, 32],
-                        decoder_layer_sizes=[32, 64], latent_dim=10,
-                        alpha=0.0001, use_mmd=True, beta=1,
-                        output_activation="ReLU")
+    model = trvaep.CVAE(
+        adata.n_vars,
+        num_classes=n_batches,
+        encoder_layer_sizes=[64, 32],
+        decoder_layer_sizes=[32, 64],
+        latent_dim=10,
+        alpha=0.0001,
+        use_mmd=True,
+        beta=1,
+        output_activation="ReLU",
+    )
 
     # Note: set seed for reproducibility of results
     trainer = trvaep.Trainer(model, adata, condition_key=batch, seed=42)
@@ -140,13 +146,13 @@ def trvaep(adata, batch, hvg=None):
 
     # Get latent representation
     latent_y = model.get_y(
-        adata.X, c=model.label_encoder.transform(
-            np.tile(np.array([main_batch]), len(adata))))
-    adata.obsm['X_emb'] = latent_y
+        adata.X,
+        c=model.label_encoder.transform(np.tile(np.array([main_batch]), len(adata))),
+    )
+    adata.obsm["X_emb"] = latent_y
 
     # Get reconstructed feature space:
-    data = model.predict(x=adata.X, y=adata.obs[batch].tolist(),
-                         target=main_batch)
+    data = model.predict(x=adata.X, y=adata.obs[batch].tolist(), target=main_batch)
     adata.X = data
 
     return adata
@@ -177,21 +183,10 @@ def scgen(adata, batch, cell_type, epochs=100, hvg=None, model_path=None, **kwar
         model_path = temp_dir.name
 
     # Fit the model
-    network = scgen.VAEArith(
-        x_dimension=adata.shape[1],
-        model_path=model_path
-    )
-    network.train(
-        train_data=adata,
-        n_epochs=epochs,
-        save=False
-    )
+    network = scgen.VAEArith(x_dimension=adata.shape[1], model_path=model_path)
+    network.train(train_data=adata, n_epochs=epochs, save=False)
     corrected_adata = scgen.batch_removal(
-        network,
-        adata,
-        batch_key=batch,
-        cell_label_key=cell_type,
-        **kwargs
+        network, adata, batch_key=batch, cell_label_key=cell_type, **kwargs
     )
 
     network.sess.close()
@@ -224,8 +219,10 @@ def scvi(adata, batch, hvg=None):
     utils.check_sanity(adata, batch, hvg)
 
     # Check for counts data layer
-    if 'counts' not in adata.layers:
-        raise TypeError('Adata does not contain a `counts` layer in `adata.layers[`counts`]`')
+    if "counts" not in adata.layers:
+        raise TypeError(
+            "Adata does not contain a `counts` layer in `adata.layers[`counts`]`"
+        )
 
     # Defaults from SCVI github tutorials scanpy_pbmc3k and harmonization
     n_epochs = np.min([round((20000 / adata.n_obs) * 400), 400])
@@ -234,20 +231,20 @@ def scvi(adata, batch, hvg=None):
     n_layers = 2
 
     net_adata = adata.copy()
-    net_adata.X = adata.layers['counts']
-    del net_adata.layers['counts']
+    net_adata.X = adata.layers["counts"]
+    del net_adata.layers["counts"]
     # Ensure that the raw counts are not accidentally used
     del net_adata.raw  # Note that this only works from anndata 0.7
 
     # Define batch indices
     le = LabelEncoder()
-    net_adata.obs['batch_indices'] = le.fit_transform(net_adata.obs[batch].values)
+    net_adata.obs["batch_indices"] = le.fit_transform(net_adata.obs[batch].values)
 
     net_adata = AnnDatasetFromAnnData(net_adata)
 
     vae = VAE(
         net_adata.nb_genes,
-        reconstruction_loss='nb',
+        reconstruction_loss="nb",
         n_batch=net_adata.n_batches,
         n_layers=n_layers,
         n_latent=n_latent,
@@ -263,10 +260,12 @@ def scvi(adata, batch, hvg=None):
 
     trainer.train(n_epochs=n_epochs, lr=1e-3)
 
-    full = trainer.create_posterior(trainer.model, net_adata, indices=np.arange(len(net_adata)))
+    full = trainer.create_posterior(
+        trainer.model, net_adata, indices=np.arange(len(net_adata))
+    )
     latent, _, _ = full.sequential().get_latent()
 
-    adata.obsm['X_emb'] = latent
+    adata.obsm["X_emb"] = latent
 
     return adata
 
@@ -295,25 +294,31 @@ def scanvi(adata, batch, labels):
 
     import numpy as np
 
-    if 'counts' not in adata.layers:
-        raise TypeError('Adata does not contain a `counts` layer in `adata.layers[`counts`]`')
+    if "counts" not in adata.layers:
+        raise TypeError(
+            "Adata does not contain a `counts` layer in `adata.layers[`counts`]`"
+        )
     # Check for counts data layer
 
     # STEP 1: prepare the data
     net_adata = adata.copy()
-    net_adata.X = adata.layers['counts']
-    del net_adata.layers['counts']
+    net_adata.X = adata.layers["counts"]
+    del net_adata.layers["counts"]
     # Ensure that the raw counts are not accidentally used
     del net_adata.raw  # Note that this only works from anndata 0.7
 
     # Define batch indices
     le = LabelEncoder()
-    net_adata.obs['batch_indices'] = le.fit_transform(net_adata.obs[batch].values)
-    net_adata.obs['labels'] = le.fit_transform(net_adata.obs[labels].values)
+    net_adata.obs["batch_indices"] = le.fit_transform(net_adata.obs[batch].values)
+    net_adata.obs["labels"] = le.fit_transform(net_adata.obs[labels].values)
 
     net_adata = AnnDatasetFromAnnData(net_adata)
 
-    print("scANVI dataset object with {} batches and {} cell types".format(net_adata.n_batches, net_adata.n_labels))
+    print(
+        "scANVI dataset object with {} batches and {} cell types".format(
+            net_adata.n_batches, net_adata.n_labels
+        )
+    )
 
     # if hvg is True:
     #    # this also corrects for different batches by default
@@ -321,7 +326,7 @@ def scanvi(adata, batch, labels):
 
     # # Defaults from SCVI github tutorials scanpy_pbmc3k and harmonization
     n_epochs_scVI = np.min([round((20000 / adata.n_obs) * 400), 400])  # 400
-    n_epochs_scANVI = int(np.min([10, np.max([2, round(n_epochs_scVI / 3.)])]))
+    n_epochs_scANVI = int(np.min([10, np.max([2, round(n_epochs_scVI / 3.0)])]))
     n_latent = 30
     n_hidden = 128
     n_layers = 2
@@ -330,7 +335,7 @@ def scanvi(adata, batch, labels):
 
     vae = VAE(
         net_adata.nb_genes,
-        reconstruction_loss='nb',
+        reconstruction_loss="nb",
         n_batch=net_adata.n_batches,
         n_latent=n_latent,
         n_hidden=n_hidden,
@@ -348,25 +353,35 @@ def scanvi(adata, batch, labels):
 
     # STEP 3: RUN scANVI
 
-    scanvi = SCANVI(net_adata.nb_genes, net_adata.n_batches, net_adata.n_labels,
-                    n_hidden=n_hidden, n_latent=n_latent, n_layers=n_layers, dispersion='gene',
-                    reconstruction_loss='nb')
+    scanvi = SCANVI(
+        net_adata.nb_genes,
+        net_adata.n_batches,
+        net_adata.n_labels,
+        n_hidden=n_hidden,
+        n_latent=n_latent,
+        n_layers=n_layers,
+        dispersion="gene",
+        reconstruction_loss="nb",
+    )
     scanvi.load_state_dict(trainer.model.state_dict(), strict=False)
 
     # use default parameter from semi-supervised trainer class
     trainer_scanvi = SemiSupervisedTrainer(scanvi, net_adata)
     # use all cells as labelled set
-    trainer_scanvi.labelled_set = trainer_scanvi.create_posterior(trainer_scanvi.model, net_adata,
-                                                                  indices=np.arange(len(net_adata)))
+    trainer_scanvi.labelled_set = trainer_scanvi.create_posterior(
+        trainer_scanvi.model, net_adata, indices=np.arange(len(net_adata))
+    )
     # put one cell in the unlabelled set
     trainer_scanvi.unlabelled_set = trainer_scanvi.create_posterior(indices=[0])
     trainer_scanvi.train(n_epochs=n_epochs_scANVI)
 
     # extract info from posterior
-    scanvi_full = trainer_scanvi.create_posterior(trainer_scanvi.model, net_adata, indices=np.arange(len(net_adata)))
+    scanvi_full = trainer_scanvi.create_posterior(
+        trainer_scanvi.model, net_adata, indices=np.arange(len(net_adata))
+    )
     latent, _, _ = scanvi_full.sequential().get_latent()
 
-    adata.obsm['X_emb'] = latent
+    adata.obsm["X_emb"] = latent
 
     return adata
 
@@ -422,17 +437,12 @@ def bbknn(adata, batch, hvg=None, **kwargs):
         raise IntegrationMethodNotFound(e)
 
     utils.check_sanity(adata, batch, hvg)
-    sc.pp.pca(adata, svd_solver='arpack')
+    sc.pp.pca(adata, svd_solver="arpack")
     if adata.n_obs < 1e5:
-        return bbknn.bbknn(
-            adata, batch_key=batch, copy=True, **kwargs)
+        return bbknn.bbknn(adata, batch_key=batch, copy=True, **kwargs)
     if adata.n_obs >= 1e5:
         return bbknn.bbknn(
-            adata,
-            batch_key=batch,
-            neighbors_within_batch=25,
-            copy=True,
-            **kwargs
+            adata, batch_key=batch, neighbors_within_batch=25, copy=True, **kwargs
         )
 
 
@@ -452,6 +462,7 @@ def saucie(adata, batch):
         raise IntegrationMethodNotFound(e)
 
     import sklearn.decomposition
+
     pca_op = sklearn.decomposition.PCA(100)
     if isinstance(adata.X, sp.sparse.csr_matrix):
         expr = adata.X.A
@@ -463,8 +474,8 @@ def saucie(adata, batch):
     loader_eval = SAUCIE.Loader(data, labels=adata.obs[batch].cat.codes, shuffle=False)
     saucie.train(loader_train, steps=5000)
     ret = adata.copy()
-    ret.obsm['X_emb'] = saucie.get_reconstruction(loader_eval)[0]
-    ret.X = pca_op.inverse_transform(ret.obsm['X_emb'])
+    ret.obsm["X_emb"] = saucie.get_reconstruction(loader_eval)[0]
+    ret.X = pca_op.inverse_transform(ret.obsm["X_emb"])
 
     return ret
 
@@ -511,20 +522,22 @@ def desc(adata, batch, res=0.8, ncores=None, tmp_dir=None, use_gpu=False):
 
     adata_out = desc.scale_bygroup(adata_out, groupby=batch, max_value=6)
 
-    adata_out = desc.train(adata_out,
-                           dims=[adata.shape[1], 128, 32],
-                           tol=0.001,
-                           n_neighbors=10,
-                           batch_size=256,
-                           louvain_resolution=res,
-                           save_encoder_weights=False,
-                           save_dir=tmp_dir,
-                           do_tsne=False,
-                           use_GPU=use_gpu,
-                           num_Cores=ncores,
-                           use_ae_weights=False,
-                           do_umap=False)
-
-    adata_out.obsm['X_emb'] = adata_out.obsm['X_Embeded_z' + str(res)]
+    adata_out = desc.train(
+        adata_out,
+        dims=[adata.shape[1], 128, 32],
+        tol=0.001,
+        n_neighbors=10,
+        batch_size=256,
+        louvain_resolution=res,
+        save_encoder_weights=False,
+        save_dir=tmp_dir,
+        do_tsne=False,
+        use_GPU=use_gpu,
+        num_Cores=ncores,
+        use_ae_weights=False,
+        do_umap=False,
+    )
+
+    adata_out.obsm["X_emb"] = adata_out.obsm["X_Embeded_z" + str(res)]
 
     return adata_out
diff --git a/scib/knn_graph/README.md b/scib/knn_graph/README.md
index eb1e465e..069beafd 100644
--- a/scib/knn_graph/README.md
+++ b/scib/knn_graph/README.md
@@ -1,8 +1,10 @@
 # knn-graph function written in C++
 
-The function `knn_graph.cpp` was written by Thomas Neumann and adapted by Maren Büttner for graph-LISI. 
-Manual compilation: 
+The function `knn_graph.cpp` was written by Thomas Neumann and adapted by Maren Büttner for graph-LISI.
+Manual compilation:
+
 ```
 g++ -std=c++11 -O3 knn_graph.cpp -o knn_graph.o
 ```
+
 A more modern version with `C++17` features is planned.
diff --git a/scib/knn_graph/knn_graph.cpp b/scib/knn_graph/knn_graph.cpp
index b20a40f9..a8e5f5d9 100644
--- a/scib/knn_graph/knn_graph.cpp
+++ b/scib/knn_graph/knn_graph.cpp
@@ -354,9 +354,9 @@ int main(int argc, char* argv[]) {
    string output_prefix = argv[2];
    // The number of neighbors we are interested in
    unsigned k = stoi(argv[3]); //convert input char to integer
-   
+
    unsigned n_chunks = stoi(argv[4]);
-   unsigned limit = matrix.getRowCount();  
+   unsigned limit = matrix.getRowCount();
    unsigned len_ch;
    if (n_chunks <= 1){
 	   n_chunks = 1;
@@ -370,12 +370,12 @@ int main(int argc, char* argv[]) {
 
    //ininitialize random number generator
    random_device rd; //used to obtain seed for random number engine
-   mt19937 gen(rd()); //standard merseen_twister_engine seeded with rd() 
-   uniform_int_distribution<> dis(1, 100); //uniform int distribution between 0 and 100  
+   mt19937 gen(rd()); //standard merseen_twister_engine seeded with rd()
+   uniform_int_distribution<> dis(1, 100); //uniform int distribution between 0 and 100
    int rand_res;
    //sanity check
-   //double sum = 0; 
-   
+   //double sum = 0;
+
    //variable declaration
    ofstream distances;
    ofstream indices;
@@ -389,7 +389,7 @@ int main(int argc, char* argv[]) {
    // write all neighbors and weights to two files
     dist = output_prefix + "_distances_" + to_string(n_ch) + ".txt";
     indi = output_prefix + "_indices_" + to_string(n_ch) + ".txt";
-   
+
     distances.open(dist, ios::out | ios::binary);
     indices.open(indi, ios::out | ios::binary);
     lower = n_ch * len_ch + 1;
@@ -397,16 +397,16 @@ int main(int argc, char* argv[]) {
     // don't run over upper limit
     if (upper > limit) {
      upper = limit;
-    } 
+    }
     for (unsigned row = lower; row <= upper; row++) {
       //cout << row << endl;
       // Ignore empty rows
       if (matrix.getRow(row).empty()) {
-	  //distances << row << endl; //add index to distances file to keep order 
+	  //distances << row << endl; //add index to distances file to keep order
 	  //indices << row << endl; //add index to indices file to keep order
 	  continue;
       }
-      // use subsampling 
+      // use subsampling
       rand_res = dis(gen); //generate random number
       if (rand_res> sub){ //skip for 100-sub percent of the data
           continue;
diff --git a/scib/metrics/__init__.py b/scib/metrics/__init__.py
index 1d8ec799..13454854 100644
--- a/scib/metrics/__init__.py
+++ b/scib/metrics/__init__.py
@@ -1,3 +1,13 @@
-from .metrics import *
-from .lisi import lisi_graph
-from .pcr import pc_regression, pcr
+from .ari import ari
+from .cell_cycle import cell_cycle
+from .clustering import opt_louvain
+from .graph_connectivity import graph_connectivity
+from .highly_variable_genes import hvg_overlap
+from .isolated_labels import isolated_labels
+from .kbet import kBET
+from .lisi import clisi_graph, ilisi_graph, lisi_graph
+from .metrics import metrics, metrics_all, metrics_fast, metrics_slim
+from .nmi import nmi
+from .pcr import pc_regression, pcr, pcr_comparison
+from .silhouette import silhouette, silhouette_batch
+from .trajectory import trajectory_conservation
diff --git a/scib/metrics/ari.py b/scib/metrics/ari.py
index 61fbd546..3e1cc3ef 100644
--- a/scib/metrics/ari.py
+++ b/scib/metrics/ari.py
@@ -29,10 +29,10 @@ def ari(adata, group1, group2, implementation=None):
 
     if len(group1) != len(group2):
         raise ValueError(
-            f'different lengths in group1 ({len(group1)}) and group2 ({len(group2)})'
+            f"different lengths in group1 ({len(group1)}) and group2 ({len(group2)})"
         )
 
-    if implementation == 'sklearn':
+    if implementation == "sklearn":
         return adjusted_rand_score(group1, group2)
 
     def binom_sum(x, k=2):
diff --git a/scib/metrics/cell_cycle.py b/scib/metrics/cell_cycle.py
index 02cfb9e8..648634f7 100644
--- a/scib/metrics/cell_cycle.py
+++ b/scib/metrics/cell_cycle.py
@@ -7,16 +7,16 @@
 
 
 def cell_cycle(
-        adata_pre,
-        adata_post,
-        batch_key,
-        embed=None,
-        agg_func=np.mean,
-        organism='mouse',
-        n_comps=50,
-        verbose=False,
-        recompute_cc=True,
-        precompute_pcr_key=None
+    adata_pre,
+    adata_post,
+    batch_key,
+    embed=None,
+    agg_func=np.mean,
+    organism="mouse",
+    n_comps=50,
+    verbose=False,
+    recompute_cc=True,
+    precompute_pcr_key=None,
 ):
     """Cell cycle conservation score
 
@@ -51,7 +51,7 @@ def cell_cycle(
     check_adata(adata_pre)
     check_adata(adata_post)
 
-    if embed == 'X_pca':
+    if embed == "X_pca":
         embed = None
 
     batches = adata_pre.obs[batch_key].unique()
@@ -59,11 +59,14 @@ def cell_cycle(
     scores_before = []
     scores_after = []
 
-    recompute_cc = recompute_cc \
-                   or 'S_score' not in adata_pre.obs_keys() \
-                   or 'G2M_score' not in adata_pre.obs_keys()
-    recompute_pcr = precompute_pcr_key is None \
-                    or precompute_pcr_key not in adata_pre.uns_keys()
+    recompute_cc = (
+        recompute_cc
+        or "S_score" not in adata_pre.obs_keys()
+        or "G2M_score" not in adata_pre.obs_keys()
+    )
+    recompute_pcr = (
+        precompute_pcr_key is None or precompute_pcr_key not in adata_pre.uns_keys()
+    )
 
     for batch in batches:
         before, after = get_pcr_before_after(
@@ -77,7 +80,7 @@ def cell_cycle(
             recompute_pcr=recompute_pcr,
             pcr_key=precompute_pcr_key,
             n_comps=n_comps,
-            verbose=verbose
+            verbose=verbose,
         )
 
         # scale result
@@ -93,7 +96,9 @@ def cell_cycle(
             score = 0
 
         if verbose:
-            print(f"batch: {batch}\t before: {before}\t after: {after}\t score: {score}")
+            print(
+                f"batch: {batch}\t before: {before}\t after: {after}\t score: {score}"
+            )
 
         scores_before.append(before)
         scores_after.append(after)
@@ -102,24 +107,24 @@ def cell_cycle(
     if agg_func is None:
         return pd.DataFrame(
             [batches, scores_before, scores_after, scores_final],
-            columns=['batch', 'before', 'after', 'score']
+            columns=["batch", "before", "after", "score"],
         )
     else:
         return agg_func(scores_final)
 
 
 def get_pcr_before_after(
-        adata_pre,
-        adata_post,
-        batch_key,
-        batch,
-        embed,
-        organism,
-        recompute_cc,
-        recompute_pcr,
-        pcr_key,
-        n_comps,
-        verbose
+    adata_pre,
+    adata_post,
+    batch_key,
+    batch,
+    embed,
+    organism,
+    recompute_cc,
+    recompute_pcr,
+    pcr_key,
+    n_comps,
+    verbose,
 ):
     """
     Principle component regression value on cell cycle scores for one batch
@@ -151,8 +156,8 @@ def get_pcr_before_after(
     if raw_sub.shape[0] != int_sub.shape[0]:
         raise ValueError(
             f'batch "{batch}" of batch_key "{batch_key}" has unequal number of '
-            f'entries before and after integration.\n'
-            f'before: {raw_sub.shape[0]} after: {int_sub.shape[0]}'
+            f"entries before and after integration.\n"
+            f"before: {raw_sub.shape[0]} after: {int_sub.shape[0]}"
         )
 
     # compute cc scores if necessary
@@ -162,27 +167,19 @@ def get_pcr_before_after(
         score_cell_cycle(raw_sub, organism=organism)
 
     # regression variable
-    covariate = raw_sub.obs[['S_score', 'G2M_score']]
+    covariate = raw_sub.obs[["S_score", "G2M_score"]]
 
     # PCR on adata before integration
     if recompute_pcr:
         before = pc_regression(
-            raw_sub.X,
-            covariate,
-            pca_var=None,
-            n_comps=n_comps,
-            verbose=verbose
+            raw_sub.X, covariate, pca_var=None, n_comps=n_comps, verbose=verbose
         )
     else:
         before = pd.Series(raw_sub.uns[pcr_key])
 
     # PCR on adata after integration
     after = pc_regression(
-        int_sub,
-        covariate,
-        pca_var=None,
-        n_comps=n_comps,
-        verbose=verbose
+        int_sub, covariate, pca_var=None, n_comps=n_comps, verbose=verbose
     )
 
     return before, after
diff --git a/scib/metrics/clustering.py b/scib/metrics/clustering.py
index c6acb3ed..cef1dade 100644
--- a/scib/metrics/clustering.py
+++ b/scib/metrics/clustering.py
@@ -7,17 +7,17 @@
 
 
 def opt_louvain(
-        adata,
-        label_key,
-        cluster_key,
-        function=None,
-        resolutions=None,
-        use_rep=None,
-        inplace=True,
-        plot=False,
-        force=True,
-        verbose=True,
-        **kwargs
+    adata,
+    label_key,
+    cluster_key,
+    function=None,
+    resolutions=None,
+    use_rep=None,
+    inplace=True,
+    plot=False,
+    force=True,
+    verbose=True,
+    **kwargs,
 ):
     """Optimised Louvain clustering
 
@@ -44,19 +44,23 @@ def opt_louvain(
     """
 
     if verbose:
-        print('Clustering...')
+        print("Clustering...")
 
     if function is None:
         function = nmi
 
     if cluster_key in adata.obs.columns:
         if force:
-            print(f"Warning: cluster key {cluster_key} already exists "
-                  "in adata.obs and will be overwritten")
+            print(
+                f"Warning: cluster key {cluster_key} already exists "
+                "in adata.obs and will be overwritten"
+            )
         else:
-            raise ValueError(f"cluster key {cluster_key} already exists in " +
-                             "adata, please remove the key or choose a different name." +
-                             "If you want to force overwriting the key, specify `force=True`")
+            raise ValueError(
+                f"cluster key {cluster_key} already exists in "
+                + "adata, please remove the key or choose a different name."
+                + "If you want to force overwriting the key, specify `force=True`"
+            )
 
     if resolutions is None:
         n = 20
@@ -68,17 +72,17 @@ def opt_louvain(
     score_all = []
 
     try:
-        adata.uns['neighbors']
+        adata.uns["neighbors"]
     except KeyError:
         if verbose:
-            print('computing neighbours for opt_cluster')
+            print("computing neighbours for opt_cluster")
         sc.pp.neighbors(adata, use_rep=use_rep)
 
     for res in resolutions:
         sc.tl.louvain(adata, resolution=res, key_added=cluster_key)
         score = function(adata, label_key, cluster_key, **kwargs)
         if verbose:
-            print(f'resolution: {res}, {function.__name__}: {score}')
+            print(f"resolution: {res}, {function.__name__}: {score}")
         score_all.append(score)
         if score_max < score:
             score_max = score
@@ -87,14 +91,18 @@ def opt_louvain(
         del adata.obs[cluster_key]
 
     if verbose:
-        print(f'optimised clustering against {label_key}')
-        print(f'optimal cluster resolution: {res_max}')
-        print(f'optimal score: {score_max}')
+        print(f"optimised clustering against {label_key}")
+        print(f"optimal cluster resolution: {res_max}")
+        print(f"optimal score: {score_max}")
 
-    score_all = pd.DataFrame(zip(resolutions, score_all), columns=('resolution', 'score'))
+    score_all = pd.DataFrame(
+        zip(resolutions, score_all), columns=("resolution", "score")
+    )
     if plot:
         # score vs. resolution profile
-        sns.lineplot(data=score_all, x='resolution', y='score').set_title('Optimal cluster resolution profile')
+        sns.lineplot(data=score_all, x="resolution", y="score").set_title(
+            "Optimal cluster resolution profile"
+        )
         plt.show()
 
     if inplace:
diff --git a/scib/metrics/graph_connectivity.py b/scib/metrics/graph_connectivity.py
index 8aee881e..fbdc24d3 100644
--- a/scib/metrics/graph_connectivity.py
+++ b/scib/metrics/graph_connectivity.py
@@ -19,9 +19,9 @@ def graph_connectivity(adata, label_key):
     :param adata: adata with computed neighborhood graph
     :param label_key: name in adata.obs containing the cell identity labels
     """
-    if 'neighbors' not in adata.uns:
+    if "neighbors" not in adata.uns:
         raise KeyError(
-            'Please compute the neighborhood graph before running this function!'
+            "Please compute the neighborhood graph before running this function!"
         )
 
     clust_res = []
@@ -29,8 +29,7 @@ def graph_connectivity(adata, label_key):
     for label in adata.obs[label_key].cat.categories:
         adata_sub = adata[adata.obs[label_key].isin([label])]
         _, labels = connected_components(
-            adata_sub.obsp['connectivities'],
-            connection='strong'
+            adata_sub.obsp["connectivities"], connection="strong"
         )
         tab = pd.value_counts(labels)
         clust_res.append(tab.max() / sum(tab))
diff --git a/scib/metrics/highly_variable_genes.py b/scib/metrics/highly_variable_genes.py
index 86668413..65144077 100644
--- a/scib/metrics/highly_variable_genes.py
+++ b/scib/metrics/highly_variable_genes.py
@@ -7,7 +7,7 @@
 def precompute_hvg_batch(adata, batch, features, n_hvg=500, save_hvg=False):
     """
     Compute HVGs per batch
-    
+
     :param adata: anndata object
     :param batch: key in adata.obs
     :param features: features to subset to
@@ -22,13 +22,15 @@ def precompute_hvg_batch(adata, batch, features, n_hvg=500, save_hvg=False):
         sc.pp.filter_genes(i, min_cells=1)
         n_hvg_tmp = np.minimum(n_hvg, int(0.5 * i.n_vars))
         if n_hvg_tmp < n_hvg:
-            print(i.obs[batch][0] + ' has less than the specified number of genes')
-            print('Number of genes: ' + str(i.n_vars))
-        hvg = sc.pp.highly_variable_genes(i, flavor='cell_ranger', n_top_genes=n_hvg_tmp, inplace=False)
-        hvg_dir[i.obs[batch][0]] = i.var.index[hvg['highly_variable']]
+            print(i.obs[batch][0] + " has less than the specified number of genes")
+            print("Number of genes: " + str(i.n_vars))
+        hvg = sc.pp.highly_variable_genes(
+            i, flavor="cell_ranger", n_top_genes=n_hvg_tmp, inplace=False
+        )
+        hvg_dir[i.obs[batch][0]] = i.var.index[hvg["highly_variable"]]
 
     if save_hvg:
-        adata.uns['hvg_before'] = hvg_dir
+        adata.uns["hvg_before"] = hvg_dir
     else:
         return hvg_dir
 
@@ -36,9 +38,9 @@ def precompute_hvg_batch(adata, batch, features, n_hvg=500, save_hvg=False):
 def hvg_overlap(adata_pre, adata_post, batch, n_hvg=500, verbose=False):
     """Highly variable gene overlap
 
-    Metric that computes the average percentage of overlapping 
+    Metric that computes the average percentage of overlapping
     highly variable genes per batch pre post integration.
-    
+
     :param adata_pre: Unintegrated anndata object
     :param adata_post: Integrated anndata object
     :param batch: Batch variable
@@ -64,12 +66,9 @@ def hvg_overlap(adata_pre, adata_post, batch, n_hvg=500, verbose=False):
 
         tmp_pre = hvg_pre_list[batch_var]
         hvg_post = sc.pp.highly_variable_genes(
-            ad_post,
-            flavor='cell_ranger',
-            n_top_genes=n_hvg_tmp,
-            inplace=False
+            ad_post, flavor="cell_ranger", n_top_genes=n_hvg_tmp, inplace=False
         )
-        tmp_post = ad_post.var.index[hvg_post['highly_variable']]
+        tmp_post = ad_post.var.index[hvg_post["highly_variable"]]
         n_hvg_real = np.minimum(len(tmp_pre), len(tmp_post))
         overlap.append((len(set(tmp_pre).intersection(set(tmp_post)))) / n_hvg_real)
     return np.mean(overlap)
diff --git a/scib/metrics/isolated_labels.py b/scib/metrics/isolated_labels.py
index db324e42..e62f796e 100644
--- a/scib/metrics/isolated_labels.py
+++ b/scib/metrics/isolated_labels.py
@@ -6,14 +6,14 @@
 
 
 def isolated_labels(
-        adata,
-        label_key,
-        batch_key,
-        embed,
-        cluster=True,
-        iso_threshold=None,
-        return_all=False,
-        verbose=True
+    adata,
+    label_key,
+    batch_key,
+    embed,
+    cluster=True,
+    iso_threshold=None,
+    return_all=False,
+    verbose=True,
 ):
     """Isolated label score
 
@@ -39,21 +39,12 @@ def isolated_labels(
     """
     scores = {}
     isolated_labels = get_isolated_labels(
-        adata,
-        label_key,
-        batch_key,
-        iso_threshold,
-        verbose
+        adata, label_key, batch_key, iso_threshold, verbose
     )
 
     for label in isolated_labels:
         score = score_isolated_label(
-            adata,
-            label_key,
-            label,
-            embed,
-            cluster,
-            verbose=verbose
+            adata, label_key, label, embed, cluster, verbose=verbose
         )
         scores[label] = score
     scores = pd.Series(scores)
@@ -65,13 +56,13 @@ def isolated_labels(
 
 
 def score_isolated_label(
-        adata,
-        label_key,
-        isolated_label,
-        embed,
-        cluster=True,
-        iso_label_key='iso_label',
-        verbose=False
+    adata,
+    label_key,
+    isolated_label,
+    embed,
+    cluster=True,
+    iso_label_key="iso_label",
+    verbose=False,
 ):
     """
     Compute label score for a single label
@@ -116,9 +107,11 @@ def max_f1(adata, label_key, cluster_key, label, argmax=False):
             use_rep=embed,
             function=max_f1,
             verbose=False,
-            inplace=True
+            inplace=True,
+        )
+        score = max_f1(
+            adata_tmp, label_key, iso_label_key, isolated_label, argmax=False
         )
-        score = max_f1(adata_tmp, label_key, iso_label_key, isolated_label, argmax=False)
     else:
         # AWS score between label
         adata_tmp.obs[iso_label_key] = adata_tmp.obs[label_key] == isolated_label
diff --git a/scib/metrics/kbet.py b/scib/metrics/kbet.py
index d0570832..ea05b3ea 100644
--- a/scib/metrics/kbet.py
+++ b/scib/metrics/kbet.py
@@ -4,6 +4,7 @@
 import numpy as np
 import pandas as pd
 import rpy2.rinterface_lib.callbacks
+import rpy2.rinterface_lib.embedded
 import rpy2.robjects as ro
 import scanpy as sc
 import scipy.sparse
@@ -12,18 +13,20 @@
 from ..utils import check_adata, check_batch
 from .utils import NeighborsError, diffusion_conn, diffusion_nn
 
-rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)  # Ignore R warning messages
+rpy2.rinterface_lib.callbacks.logger.setLevel(
+    logging.ERROR
+)  # Ignore R warning messages
 
 
 def kBET(
-        adata,
-        batch_key,
-        label_key,
-        scaled=True,
-        embed='X_pca',
-        type_=None,
-        return_df=False,
-        verbose=False
+    adata,
+    batch_key,
+    label_key,
+    scaled=True,
+    embed="X_pca",
+    type_=None,
+    return_df=False,
+    verbose=False,
 ):
     """kBET score
 
@@ -54,25 +57,27 @@ def kBET(
 
     # compute connectivities for non-knn type data integrations
     # and increase neighborhoods for knn type data integrations
-    if type_ != 'knn':
+    if type_ != "knn":
         adata_tmp = sc.pp.neighbors(adata, n_neighbors=50, use_rep=embed, copy=True)
     else:
         # check if pre-computed neighbours are stored in input file
         adata_tmp = adata.copy()
-        if 'diffusion_connectivities' not in adata.uns['neighbors']:
+        if "diffusion_connectivities" not in adata.uns["neighbors"]:
             if verbose:
-                print(f"Compute: Diffusion neighbours.")
+                print("Compute diffusion neighbours")
             adata_tmp = diffusion_conn(adata, min_k=50, copy=True)
-        adata_tmp.obsp['connectivities'] = adata_tmp.uns['neighbors']['diffusion_connectivities']
+        adata_tmp.obsp["connectivities"] = adata_tmp.uns["neighbors"][
+            "diffusion_connectivities"
+        ]
 
     if verbose:
         print(f"batch: {batch_key}")
 
     # set upper bound for k0
-    size_max = 2 ** 31 - 1
+    size_max = 2**31 - 1
 
     # prepare call of kBET per cluster
-    kBET_scores = {'cluster': [], 'kBET': []}
+    kBET_scores = {"cluster": [], "kBET": []}
     for clus in adata_tmp.obs[label_key].unique():
 
         # subset by label
@@ -80,40 +85,41 @@ def kBET(
 
         # check if neighborhood size too small or only one batch in subset
         if np.logical_or(
-                adata_sub.n_obs < 10,
-                len(adata_sub.obs[batch_key].cat.categories) == 1
+            adata_sub.n_obs < 10, len(adata_sub.obs[batch_key].cat.categories) == 1
         ):
             print(f"{clus} consists of a single batch or is too small. Skip.")
             score = np.nan
         else:
-            quarter_mean = np.floor(np.mean(adata_sub.obs[batch_key].value_counts()) / 4).astype('int')
+            quarter_mean = np.floor(
+                np.mean(adata_sub.obs[batch_key].value_counts()) / 4
+            ).astype("int")
             k0 = np.min([70, np.max([10, quarter_mean])])
             # check k0 for reasonability
             if k0 * adata_sub.n_obs >= size_max:
-                k0 = np.floor(size_max / adata_sub.n_obs).astype('int')
+                k0 = np.floor(size_max / adata_sub.n_obs).astype("int")
 
             matrix = np.zeros(shape=(adata_sub.n_obs, k0 + 1))
 
             if verbose:
                 print(f"Use {k0} nearest neighbors.")
             n_comp, labs = scipy.sparse.csgraph.connected_components(
-                adata_sub.obsp['connectivities'],
-                connection='strong'
+                adata_sub.obsp["connectivities"], connection="strong"
             )
 
             if n_comp == 1:  # a single component to compute kBET on
                 try:
-                    nn_index_tmp = diffusion_nn(adata_sub, k=k0).astype('float')
+                    nn_index_tmp = diffusion_nn(adata_sub, k=k0).astype("float")
                     # call kBET
                     score = kBET_single(
                         matrix=matrix,
                         batch=adata_sub.obs[batch_key],
-                        knn=nn_index_tmp + 1,  # nn_index in python is 0-based and 1-based in R
+                        knn=nn_index_tmp
+                        + 1,  # nn_index in python is 0-based and 1-based in R
                         verbose=verbose,
-                        k0=k0
+                        k0=k0,
                     )
                 except NeighborsError:
-                    print('Not enough neighbours')
+                    print("Not enough neighbours")
                     score = 1  # i.e. 100% rejection
 
             else:
@@ -133,23 +139,26 @@ def kBET(
                     nn_index_tmp[:] = np.nan
 
                     try:
-                        nn_index_tmp[idx_nonan] = diffusion_nn(adata_sub_sub, k=k0).astype('float')
+                        nn_index_tmp[idx_nonan] = diffusion_nn(
+                            adata_sub_sub, k=k0
+                        ).astype("float")
                         # call kBET
                         score = kBET_single(
                             matrix=matrix,
                             batch=adata_sub.obs[batch_key],
-                            knn=nn_index_tmp + 1,  # nn_index in python is 0-based and 1-based in R
+                            knn=nn_index_tmp
+                            + 1,  # nn_index in python is 0-based and 1-based in R
                             verbose=verbose,
-                            k0=k0
+                            k0=k0,
                         )
                     except NeighborsError:
-                        print('Not enough neighbours')
+                        print("Not enough neighbours")
                         score = 1  # i.e. 100% rejection
                 else:  # if there are too many too small connected components, set kBET score to 1
                     score = 1  # i.e. 100% rejection
 
-        kBET_scores['cluster'].append(clus)
-        kBET_scores['kBET'].append(score)
+        kBET_scores["cluster"].append(clus)
+        kBET_scores["kBET"].append(score)
 
     kBET_scores = pd.DataFrame.from_dict(kBET_scores)
     kBET_scores = kBET_scores.reset_index(drop=True)
@@ -157,17 +166,11 @@ def kBET(
     if return_df:
         return kBET_scores
 
-    final_score = np.nanmean(kBET_scores['kBET'])
+    final_score = np.nanmean(kBET_scores["kBET"])
     return 1 - final_score if scaled else final_score
 
 
-def kBET_single(
-        matrix,
-        batch,
-        k0=10,
-        knn=None,
-        verbose=False
-):
+def kBET_single(matrix, batch, k0=10, knn=None, verbose=False):
     """Single kBET run
 
     Compute k-nearest neighbour batch effect test (kBET) score as described in
@@ -179,21 +182,21 @@ def kBET_single(
     """
     try:
         ro.r("library(kBET)")
-    except Exception as ex:
+    except rpy2.rinterface_lib.embedded.RRuntimeError as ex:
         RLibraryNotFound(ex)
 
     anndata2ri.activate()
 
     if verbose:
         print("importing expression matrix")
-    ro.globalenv['data_mtrx'] = matrix
-    ro.globalenv['batch'] = batch
+    ro.globalenv["data_mtrx"] = matrix
+    ro.globalenv["batch"] = batch
 
     if verbose:
         print("kBET estimation")
 
-    ro.globalenv['knn_graph'] = knn
-    ro.globalenv['k0'] = k0
+    ro.globalenv["knn_graph"] = knn
+    ro.globalenv["k0"] = k0
     ro.r(
         "batch.estimate <- kBET("
         "  data_mtrx,"
@@ -210,7 +213,8 @@ def kBET_single(
 
     try:
         score = ro.r("batch.estimate$summary$kBET.observed")[0]
-    except rpy2.rinterface_lib.embedded.RRuntimeError:
+    except rpy2.rinterface_lib.embedded.RRuntimeError as ex:
+        print(f"Error computing kBET: {ex}\nSetting value to np.nan")
         score = np.nan
 
     anndata2ri.deactivate()
diff --git a/scib/metrics/lisi.py b/scib/metrics/lisi.py
index 70f685a0..82e23306 100644
--- a/scib/metrics/lisi.py
+++ b/scib/metrics/lisi.py
@@ -10,6 +10,7 @@
 import numpy as np
 import pandas as pd
 import rpy2.rinterface_lib.callbacks
+import rpy2.rinterface_lib.embedded
 import rpy2.robjects as ro
 import scanpy as sc
 import scipy.sparse
@@ -17,19 +18,17 @@
 from scipy.io import mmwrite
 
 import scib
+
 from ..exceptions import RLibraryNotFound
 from ..utils import check_adata, check_batch
 
-rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)  # Ignore R warning messages
+rpy2.rinterface_lib.callbacks.logger.setLevel(
+    logging.ERROR
+)  # Ignore R warning messages
 
 
 # Graph LISI (analoguous to lisi function)
-def lisi_graph(
-        adata,
-        batch_key,
-        label_key,
-        **kwargs
-):
+def lisi_graph(adata, batch_key, label_key, **kwargs):
     """cLISI and iLISI scores
 
     This is a reimplementation of the LISI (Local Inverse Simpson’s Index) metrics
@@ -49,14 +48,14 @@ def lisi_graph(
 
 
 def ilisi_graph(
-        adata,
-        batch_key,
-        k0=90,
-        type_=None,
-        subsample=None,
-        scale=True,
-        n_cores=1,
-        verbose=False
+    adata,
+    batch_key,
+    k0=90,
+    type_=None,
+    subsample=None,
+    scale=True,
+    n_cores=1,
+    verbose=False,
 ):
     """Integration LISI (iLISI) score
 
@@ -71,7 +70,7 @@ def ilisi_graph(
     :param k0: number of nearest neighbors to compute lisi score
         Please note that the initial neighborhood size that is
         used to compute shortest paths is 15.
-    :param type_: type of data integration, either knn, full or embed
+    :param `type_`: type of data integration, either knn, full or embed
     :param subsample: Percentage of observations (integer between 0 and 100)
         to which lisi scoring should be subsampled
     :param scale: scale output values between 0 and 1 (True/False)
@@ -90,7 +89,7 @@ def ilisi_graph(
         perplexity=None,
         subsample=subsample,
         n_cores=n_cores,
-        verbose=verbose
+        verbose=verbose,
     )
 
     # iLISI: nbatches good, 1 bad
@@ -104,15 +103,15 @@ def ilisi_graph(
 
 
 def clisi_graph(
-        adata,
-        batch_key,  # TODO: remove
-        label_key,
-        k0=90,
-        type_=None,
-        subsample=None,
-        scale=True,
-        n_cores=1,
-        verbose=False
+    adata,
+    batch_key,  # TODO: remove
+    label_key,
+    k0=90,
+    type_=None,
+    subsample=None,
+    scale=True,
+    n_cores=1,
+    verbose=False,
 ):
     """Cell-type LISI (cLISI) score
 
@@ -128,7 +127,7 @@ def clisi_graph(
     :param k0: number of nearest neighbors to compute lisi score
         Please note that the initial neighborhood size that is
         used to compute shortest paths is 15.
-    :param type_: type of data integration, either knn, full or embed
+    :param `type_`: type of data integration, either knn, full or embed
     :param subsample: Percentage of observations (integer between 0 and 100)
         to which lisi scoring should be subsampled
     :param scale: scale output values between 0 and 1 (True/False)
@@ -149,7 +148,7 @@ def clisi_graph(
         perplexity=None,
         subsample=subsample,
         n_cores=n_cores,
-        verbose=verbose
+        verbose=verbose,
     )
 
     # cLISI: 1 good, nlabs bad
@@ -163,13 +162,12 @@ def clisi_graph(
 
 
 def recompute_knn(adata, type_):
-    """Recompute neighbours
-    """
-    if type_ == 'embed':
-        return sc.pp.neighbors(adata, n_neighbors=15, use_rep='X_emb', copy=True)
-    elif type_ == 'full':
-        if 'X_pca' not in adata.obsm.keys():
-            sc.pp.pca(adata, svd_solver='arpack')
+    """Recompute neighbours"""
+    if type_ == "embed":
+        return sc.pp.neighbors(adata, n_neighbors=15, use_rep="X_emb", copy=True)
+    elif type_ == "full":
+        if "X_pca" not in adata.obsm.keys():
+            sc.pp.pca(adata, svd_solver="arpack")
         return sc.pp.neighbors(adata, n_neighbors=15, copy=True)
     else:
         # if knn - do not compute a new neighbourhood graph (it exists already)
@@ -177,13 +175,13 @@ def recompute_knn(adata, type_):
 
 
 def lisi_graph_py(
-        adata,
-        batch_key,
-        n_neighbors=90,
-        perplexity=None,
-        subsample=None,
-        n_cores=1,
-        verbose=False
+    adata,
+    batch_key,
+    n_neighbors=90,
+    perplexity=None,
+    subsample=None,
+    n_cores=1,
+    verbose=False,
 ):
     """
     Function to prepare call of compute_simpson_index
@@ -194,9 +192,10 @@ def lisi_graph_py(
     # use no more than the available cores
     n_cores = max(1, min(n_cores, mp.cpu_count()))
 
-    if 'neighbors' not in adata.uns:
-        raise AttributeError(f"key 'neighbors' not found. Please make sure that a " +
-                             "kNN graph has been computed")
+    if "neighbors" not in adata.uns:
+        raise AttributeError(
+            "Key 'neighbors' not found. Please make sure that a kNN graph has been computed"
+        )
     elif verbose:
         print("using precomputed kNN graph")
 
@@ -215,37 +214,37 @@ def lisi_graph_py(
     subset = 100  # default, no subsampling
     if subsample is not None:
         subset = subsample  # do not use subsampling
-        if isinstance(subsample, int) == False:  # need to set as integer
+        if isinstance(subsample, int) is False:  # need to set as integer
             subset = int(subsample)
 
     # run LISI in python
     if verbose:
         print("Compute knn on shortest paths")
 
-        # set connectivities to 3e-308 if they are lower than 3e-308 (because cpp can't handle double values smaller than that).
-    connectivities = adata.obsp['connectivities']  # csr matrix format
+    # set connectivities to 3e-308 if they are lower than 3e-308 (because cpp can't handle double values smaller than that).
+    connectivities = adata.obsp["connectivities"]  # csr matrix format
     large_enough = connectivities.data >= 3e-308
     if verbose:
-        n_too_small = np.sum(large_enough == False)
+        n_too_small = np.sum(large_enough is False)
         if n_too_small:
-            print(f"{n_too_small} connectivities are smaller than 3e-308 and will be set to 3e-308")
-            print(connectivities.data[large_enough == False])
-    connectivities.data[large_enough == False] = 3e-308
+            print(
+                f"{n_too_small} connectivities are smaller than 3e-308 and will be set to 3e-308"
+            )
+            print(connectivities.data[large_enough is False])
+    connectivities.data[large_enough is False] = 3e-308
 
     # temporary file
     tmpdir = tempfile.TemporaryDirectory(prefix="lisi_")
-    prefix = tmpdir.name + '/graph_lisi'
-    mtx_file_path = prefix + '_input.mtx'
+    prefix = tmpdir.name + "/graph_lisi"
+    mtx_file_path = prefix + "_input.mtx"
 
-    mmwrite(
-        mtx_file_path,
-        connectivities,
-        symmetry='general'
-    )
+    mmwrite(mtx_file_path, connectivities, symmetry="general")
     # call knn-graph computation in Cpp
 
     root = pathlib.Path(scib.__file__).parent  # get current root directory
-    cpp_file_path = root / 'knn_graph/knn_graph.o'  # create POSIX path to file to execute compiled cpp-code
+    cpp_file_path = (
+        root / "knn_graph/knn_graph.o"
+    )  # create POSIX path to file to execute compiled cpp-code
     # comment: POSIX path needs to be converted to string - done below with 'as_posix()'
     # create evenly split chunks if n_obs is divisible by n_chunks (doesn't really make sense on 2nd thought)
     args_int = [
@@ -254,15 +253,14 @@ def lisi_graph_py(
         prefix,
         str(n_neighbors),
         str(n_cores),  # number of splits
-        str(subset)
+        str(subset),
     ]
     if verbose:
         print(f'call {" ".join(args_int)}')
     try:
         subprocess.run(args_int)
-    except Exception as e:
-        print(e)
-        print("Couldn't compute LISI, returning NaN")
+    except rpy2.rinterface_lib.embedded.RRuntimeError as ex:
+        print(f"Error computing LISI kNN graph {ex}\nSetting value to np.nan")
         return np.nan
 
     if verbose:
@@ -284,8 +282,8 @@ def lisi_graph_py(
                 itertools.repeat(n_batches),
                 itertools.repeat(n_neighbors),
                 itertools.repeat(perplexity),
-                chunk_no
-            )
+                chunk_no,
+            ),
         )
         pool.close()
         pool.join()
@@ -298,7 +296,7 @@ def lisi_graph_py(
             batch_labels=batch,
             n_batches=n_batches,
             perplexity=perplexity,
-            n_neighbors=n_neighbors
+            n_neighbors=n_neighbors,
         )
 
     tmpdir.cleanup()
@@ -308,13 +306,9 @@ def lisi_graph_py(
 
 # LISI core functions
 
+
 def compute_simpson_index(
-        D=None,
-        knn_idx=None,
-        batch_labels=None,
-        n_batches=None,
-        perplexity=15,
-        tol=1e-5
+    D=None, knn_idx=None, batch_labels=None, n_batches=None, perplexity=15, tol=1e-5
 ):
     """
     Simpson index of batch labels subset by group.
@@ -345,16 +339,16 @@ def compute_simpson_index(
         Hdiff = H - logU
         tries = 0
         # first get neighbor probabilities
-        while (np.logical_and(np.abs(Hdiff) > tol, tries < 50)):
-            if (Hdiff > 0):
+        while np.logical_and(np.abs(Hdiff) > tol, tries < 50):
+            if Hdiff > 0:
                 betamin = beta
-                if (betamax == np.inf):
+                if betamax == np.inf:
                     beta *= 2
                 else:
                     beta = (beta + betamax) / 2
             else:
                 betamax = beta
-                if (betamin == -np.inf):
+                if betamin == -np.inf:
                     beta /= 2
                 else:
                     beta = (beta + betamin) / 2
@@ -363,12 +357,12 @@ def compute_simpson_index(
             Hdiff = H - logU
             tries += 1
 
-        if (H == 0):
+        if H == 0:
             simpson[i] = -1
             continue
 
             # then compute Simpson's Index
-        non_nan_knn = knn_idx[i][np.invert(np.isnan(knn_idx[i]))].astype('int')
+        non_nan_knn = knn_idx[i][np.invert(np.isnan(knn_idx[i]))].astype("int")
         batch = batch_labels[non_nan_knn]
         # convertToOneHot omits all nan entries.
         # Therefore, we run into errors in np.matmul.
@@ -383,13 +377,13 @@ def compute_simpson_index(
 
 
 def compute_simpson_index_graph(
-        file_prefix=None,
-        batch_labels=None,
-        n_batches=None,
-        n_neighbors=90,
-        perplexity=30,
-        chunk_no=0,
-        tol=1e-5
+    file_prefix=None,
+    batch_labels=None,
+    n_batches=None,
+    n_neighbors=90,
+    perplexity=30,
+    chunk_no=0,
+    tol=1e-5,
 ):
     """
     Simpson index of batch labels subset by group.
@@ -403,8 +397,8 @@ def compute_simpson_index_graph(
     :param tol: a tolerance for testing effective neighborhood size
     :returns: the simpson index for the neighborhood of each cell
     """
-    index_file = file_prefix + '_indices_' + str(chunk_no) + '.txt'
-    distance_file = file_prefix + '_distances_' + str(chunk_no) + '.txt'
+    index_file = file_prefix + "_indices_" + str(chunk_no) + ".txt"
+    distance_file = file_prefix + "_distances_" + str(chunk_no) + ".txt"
 
     # initialize
     P = np.zeros(n_neighbors)
@@ -417,14 +411,14 @@ def compute_simpson_index_graph(
         return lists
 
     # read distances and indices with nan value handling
-    indices = pd.read_table(index_file, index_col=0, header=None, sep=',')
+    indices = pd.read_table(index_file, index_col=0, header=None, sep=",")
     indices = indices.T
 
-    distances = pd.read_table(distance_file, index_col=0, header=None, sep=',')
+    distances = pd.read_table(distance_file, index_col=0, header=None, sep=",")
     distances = distances.T
 
     # get cell ids
-    chunk_ids = indices.columns.values.astype('int')
+    chunk_ids = indices.columns.values.astype("int")
 
     # define result vector
     simpson = np.zeros(len(chunk_ids))
@@ -437,14 +431,14 @@ def compute_simpson_index_graph(
 
         if get_col.isnull().sum() > 0:
             # not enough neighbors
-            print(f'Chunk {chunk_id} does not have enough neighbors. Skipping...')
+            print(f"Chunk {chunk_id} does not have enough neighbors. Skipping...")
             simpson[i] = 1  # np.nan #set nan for testing
             continue
 
-        knn_idx = get_col.astype('int') - 1  # get 0-based indexing
+        knn_idx = get_col.astype("int") - 1  # get 0-based indexing
 
         # read line i from distances matrix
-        D_act = distances[chunk_id].values.astype('float')
+        D_act = distances[chunk_id].values.astype("float")
 
         # start lisi estimation
         beta = 1
@@ -490,9 +484,9 @@ def Hbeta(D_row, beta):
     """
     Helper function for simpson index computation
     """
-    P = np.exp(- D_row * beta)
+    P = np.exp(-D_row * beta)
     sumP = np.nansum(P)
-    if (sumP == 0):
+    if sumP == 0:
         H = 0
         P = np.zeros(len(D_row))
     else:
@@ -538,16 +532,9 @@ def convert_to_one_hot(vector, num_classes=None):
 
 # Deprecated functions
 
+
 @deprecated
-def lisi(
-        adata,
-        batch_key,
-        label_key,
-        k0=90,
-        type_=None,
-        scale=True,
-        verbose=False
-):
+def lisi(adata, batch_key, label_key, k0=90, type_=None, scale=True, verbose=False):
     """Compute iLISI and cLISI scores
 
     This is a reimplementation of the LISI (Local Inverse Simpson’s Index) metrics
@@ -567,18 +554,20 @@ def lisi(
     #    if verbose:
     #        print("recompute kNN graph with {k0} nearest neighbors.")
     # recompute neighbours
-    if (type_ == 'embed'):
-        adata_tmp = sc.pp.neighbors(adata, n_neighbors=k0, use_rep='X_emb', copy=True)
-    elif (type_ == 'full'):
-        if 'X_pca' not in adata.obsm.keys():
-            sc.pp.pca(adata, svd_solver='arpack')
+    if type_ == "embed":
+        adata_tmp = sc.pp.neighbors(adata, n_neighbors=k0, use_rep="X_emb", copy=True)
+    elif type_ == "full":
+        if "X_pca" not in adata.obsm.keys():
+            sc.pp.pca(adata, svd_solver="arpack")
         adata_tmp = sc.pp.neighbors(adata, n_neighbors=k0, copy=True)
     else:
         adata_tmp = adata.copy()
     # if knn - do not compute a new neighbourhood graph (it exists already)
 
     # lisi_score = lisi_knn(adata=adata, batch_key=batch_key, label_key=label_key, verbose=verbose)
-    lisi_score = lisi_knn_py(adata=adata_tmp, batch_key=batch_key, label_key=label_key, verbose=verbose)
+    lisi_score = lisi_knn_py(
+        adata=adata_tmp, batch_key=batch_key, label_key=label_key, verbose=verbose
+    )
 
     # iLISI: nbatches good, 1 bad
     ilisi_score = np.nanmedian(lisi_score[batch_key])
@@ -594,43 +583,36 @@ def lisi(
 
 
 @deprecated
-def lisi_knn_py(
-        adata,
-        batch_key,
-        label_key,
-        perplexity=None,
-        verbose=False
-):
+def lisi_knn_py(adata, batch_key, label_key, perplexity=None, verbose=False):
     """
     Compute LISI score on kNN graph provided in the adata object. By default, perplexity
     is chosen as 1/3 * number of nearest neighbours in the knn-graph.
     """
 
-    if 'neighbors' not in adata.uns:
-        raise AttributeError(f"key 'neighbors' not found. Please make sure that a " +
-                             "kNN graph has been computed")
+    if "neighbors" not in adata.uns:
+        raise AttributeError(
+            "key 'neighbors' not found. Please make sure that a kNN graph has been computed"
+        )
     elif verbose:
         print("using precomputed kNN graph")
 
     # get knn index matrix
     if verbose:
         print("Convert nearest neighbor matrix and distances for LISI.")
-    dist_mat = scipy.sparse.find(adata.obsp['distances'])
+    dist_mat = scipy.sparse.find(adata.obsp["distances"])
     # get number of nearest neighbours parameter
-    if 'params' not in adata.uns['neighbors']:
+    if "params" not in adata.uns["neighbors"]:
         # estimate the number of nearest neighbors as the median
         # of the distance matrix
         _, e = np.unique(dist_mat[0], return_counts=True)
         n_nn = np.nanmedian(e)
-        n_nn = n_nn.astype('int')
+        n_nn = n_nn.astype("int")
     else:
-        n_nn = adata.uns['neighbors']['params']['n_neighbors'] - 1
+        n_nn = adata.uns["neighbors"]["params"]["n_neighbors"] - 1
     # initialise index and fill it with NaN values
-    nn_index = np.empty(shape=(adata.obsp['distances'].shape[0],
-                               n_nn))
+    nn_index = np.empty(shape=(adata.obsp["distances"].shape[0], n_nn))
     nn_index[:] = np.NaN
-    nn_dists = np.empty(shape=(adata.obsp['distances'].shape[0],
-                               n_nn))
+    nn_dists = np.empty(shape=(adata.obsp["distances"].shape[0], n_nn))
     nn_dists[:] = np.NaN
     index_out = []
     for cell_id in np.arange(np.min(dist_mat[0]), np.max(dist_mat[0]) + 1):
@@ -638,7 +620,9 @@ def lisi_knn_py(
         num_idx = get_idx.sum()
         # in case that get_idx contains more than n_nn neighbours, cut away the outlying ones
         fin_idx = np.min([num_idx, n_nn])
-        nn_index[cell_id, :fin_idx] = dist_mat[1][get_idx][np.argsort(dist_mat[2][get_idx])][:fin_idx]
+        nn_index[cell_id, :fin_idx] = dist_mat[1][get_idx][
+            np.argsort(dist_mat[2][get_idx])
+        ][:fin_idx]
         nn_dists[cell_id, :fin_idx] = np.sort(dist_mat[2][get_idx])[:fin_idx]
         if num_idx < n_nn:
             index_out.append(cell_id)
@@ -665,18 +649,20 @@ def lisi_knn_py(
     if verbose:
         print("LISI score estimation")
 
-    simpson_estimate_batch = compute_simpson_index(D=nn_dists,
-                                                   knn_idx=nn_index,
-                                                   batch_labels=batch,
-                                                   n_batches=n_batches,
-                                                   perplexity=perplexity,
-                                                   )
-    simpson_estimate_label = compute_simpson_index(D=nn_dists,
-                                                   knn_idx=nn_index,
-                                                   batch_labels=label,
-                                                   n_batches=n_labels,
-                                                   perplexity=perplexity
-                                                   )
+    simpson_estimate_batch = compute_simpson_index(
+        D=nn_dists,
+        knn_idx=nn_index,
+        batch_labels=batch,
+        n_batches=n_batches,
+        perplexity=perplexity,
+    )
+    simpson_estimate_label = compute_simpson_index(
+        D=nn_dists,
+        knn_idx=nn_index,
+        batch_labels=label,
+        n_batches=n_labels,
+        perplexity=perplexity,
+    )
     simpson_est_batch = 1 / simpson_estimate_batch
     simpson_est_label = 1 / simpson_estimate_label
     # extract results
@@ -691,24 +677,20 @@ def scale_lisi(ilisi_score, clisi_score, nbatches):
     # scale iLISI score to 0 bad 1 good
     ilisi_score = (ilisi_score - 1) / (nbatches - 1)
     # scale clisi score to 0 bad 1 good
-    clisi_score = (nbatches - clisi_score) / (nbatches - 1)  # Scaled incorrectly by n_batches
+    clisi_score = (nbatches - clisi_score) / (
+        nbatches - 1
+    )  # Scaled incorrectly by n_batches
     return ilisi_score, clisi_score
 
 
 @deprecated
-def lisi_knn(
-        adata,
-        batch_key,
-        label_key,
-        perplexity=None,
-        verbose=False
-):
+def lisi_knn(adata, batch_key, label_key, perplexity=None, verbose=False):
     """
     Compute LISI score on kNN graph provided in the adata object. By default, perplexity
     is chosen as 1/3 * number of nearest neighbours in the knn-graph.
     """
 
-    if 'neighbors' not in adata.uns:
+    if "neighbors" not in adata.uns:
         raise AttributeError(
             "key 'neighbors' not found. Please make sure that a "
             "kNN graph has been computed"
@@ -719,20 +701,18 @@ def lisi_knn(
     # get knn index matrix
     if verbose:
         print("Convert nearest neighbor matrix and distances for LISI.")
-    dist_mat = scipy.sparse.find(adata.obsp['distances'])
+    dist_mat = scipy.sparse.find(adata.obsp["distances"])
     # get number of nearest neighbours parameter
-    if 'params' not in adata.uns['neighbors']:
+    if "params" not in adata.uns["neighbors"]:
         # estimate the number of nearest neighbors as the median
         # of the distance matrix
         _, e = np.unique(dist_mat[0], return_counts=True)
         n_nn = np.nanmin(e)
-        n_nn = n_nn.astype('int')
+        n_nn = n_nn.astype("int")
     else:
-        n_nn = adata.uns['neighbors']['params']['n_neighbors'] - 1
-    nn_index = np.empty(shape=(adata.obsp['distances'].shape[0],
-                               n_nn))
-    nn_dists = np.empty(shape=(adata.obsp['distances'].shape[0],
-                               n_nn))
+        n_nn = adata.uns["neighbors"]["params"]["n_neighbors"] - 1
+    nn_index = np.empty(shape=(adata.obsp["distances"].shape[0], n_nn))
+    nn_dists = np.empty(shape=(adata.obsp["distances"].shape[0], n_nn))
     index_out = []
     for cell_id in np.arange(np.min(dist_mat[0]), np.max(dist_mat[0]) + 1):
         get_idx = dist_mat[0] == cell_id
@@ -740,7 +720,9 @@ def lisi_knn(
         # in case that get_idx contains more than n_nn neighbours, cut away the outlying ones
         # potential enhancement: handle case where less than n_nn neighbours are reported
         if num_idx >= n_nn:
-            nn_index[cell_id, :] = dist_mat[1][get_idx][np.argsort(dist_mat[2][get_idx])][:n_nn]
+            nn_index[cell_id, :] = dist_mat[1][get_idx][
+                np.argsort(dist_mat[2][get_idx])
+            ][:n_nn]
             nn_dists[cell_id, :] = np.sort(dist_mat[2][get_idx])[:n_nn]
         else:
             index_out.append(cell_id)
@@ -753,7 +735,9 @@ def lisi_knn(
         # nn_dists = np.delete(nn_dists, index_out, 0)
         # nn_index = np.delete(nn_index, index_out, 0)
         if verbose:
-            print(f"{out_cells} had less than {n_nn} neighbors and were omitted in LISI score.")
+            print(
+                f"{out_cells} had less than {n_nn} neighbors and were omitted in LISI score."
+            )
 
     if perplexity is None:
         # use LISI default
@@ -762,27 +746,29 @@ def lisi_knn(
     # run LISI in R
     try:
         ro.r("library(lisi)")
-    except Exception as ex:
+    except rpy2.rinterface_lib.embedded.RRuntimeError as ex:
         RLibraryNotFound(ex)
 
     anndata2ri.activate()
 
     if verbose:
         print("importing knn-graph")
-    ro.globalenv['nn_indx'] = nn_index.astype('int').T
-    ro.globalenv['nn_dst'] = nn_dists.T
-    ro.globalenv['perplexity'] = perplexity
-    ro.globalenv['batch'] = adata.obs[batch_key].cat.codes.values
-    ro.globalenv['n_batches'] = len(np.unique(adata.obs[batch_key]))
-    ro.globalenv['label'] = adata.obs[label_key].cat.codes.values
-    ro.globalenv['n_labels'] = len(np.unique(adata.obs[label_key]))
+    ro.globalenv["nn_indx"] = nn_index.astype("int").T
+    ro.globalenv["nn_dst"] = nn_dists.T
+    ro.globalenv["perplexity"] = perplexity
+    ro.globalenv["batch"] = adata.obs[batch_key].cat.codes.values
+    ro.globalenv["n_batches"] = len(np.unique(adata.obs[batch_key]))
+    ro.globalenv["label"] = adata.obs[label_key].cat.codes.values
+    ro.globalenv["n_labels"] = len(np.unique(adata.obs[label_key]))
 
     if verbose:
         print("LISI score estimation")
-    simpson_estimate_batch = ro.r(
-        f"simpson.estimate_batch <- compute_simpson_index(nn_dst, nn_indx, batch, n_batches, perplexity)")  # batch_label_keys)")
-    simpson_estimate_label = ro.r(
-        f"simpson.estimate_label <- compute_simpson_index(nn_dst, nn_indx, label, n_labels, perplexity)")  # batch_label_keys)")
+    ro.r(
+        "simpson.estimate_batch <- compute_simpson_index(nn_dst, nn_indx, batch, n_batches, perplexity)"
+    )  # batch_label_keys)")
+    ro.r(
+        "simpson.estimate_label <- compute_simpson_index(nn_dst, nn_indx, label, n_labels, perplexity)"
+    )  # batch_label_keys)")
     simpson_est_batch = 1 / np.squeeze(ro.r("simpson.estimate_batch"))
     simpson_est_label = 1 / np.squeeze(ro.r("simpson.estimate_label"))
 
@@ -796,13 +782,7 @@ def lisi_knn(
 
 
 @deprecated
-def lisi_matrix(
-        adata,
-        batch_key,
-        label_key,
-        matrix=None,
-        verbose=False
-):
+def lisi_matrix(adata, batch_key, label_key, matrix=None, verbose=False):
     """
     Computes the LISI scores for a given data matrix in adata.X. The scoring function of the
     LISI R package is called with default parameters. This function takes a data matrix and
@@ -819,25 +799,27 @@ def lisi_matrix(
     # run LISI in R
     try:
         ro.r("library(lisi)")
-    except Exception as ex:
+    except rpy2.rinterface_lib.embedded.RRuntimeError as ex:
         RLibraryNotFound(ex)
 
     anndata2ri.activate()
 
     if verbose:
         print("importing expression matrix")
-    ro.globalenv['data_mtrx'] = matrix
+    ro.globalenv["data_mtrx"] = matrix
 
     if verbose:
         print(f"covariates: {batch_key} and {label_key}")
     metadata = adata.obs[[batch_key, label_key]]
-    ro.globalenv['metadata'] = metadata
+    ro.globalenv["metadata"] = metadata
     batch_label_keys = ro.StrVector([batch_key, label_key])
-    ro.globalenv['batch_label_keys'] = batch_label_keys
+    ro.globalenv["batch_label_keys"] = batch_label_keys
 
     if verbose:
         print("LISI score estimation")
-    lisi_estimate = ro.r(f"lisi.estimate <- compute_lisi(data_mtrx, metadata, batch_label_keys)")  # batch_label_keys)")
+    lisi_estimate = ro.r(
+        "lisi.estimate <- compute_lisi(data_mtrx, metadata, batch_label_keys)"
+    )  # batch_label_keys)")
     anndata2ri.deactivate()
 
     return lisi_estimate
diff --git a/scib/metrics/metrics.py b/scib/metrics/metrics.py
index 9fae5e95..393cf26b 100755
--- a/scib/metrics/metrics.py
+++ b/scib/metrics/metrics.py
@@ -1,6 +1,5 @@
 import numpy as np
 import pandas as pd
-from deprecated import deprecated
 
 from ..utils import check_adata, check_batch
 from .ari import ari
@@ -17,13 +16,7 @@
 from .trajectory import trajectory_conservation
 
 
-def metrics_fast(
-        adata,
-        adata_int,
-        batch_key,
-        label_key,
-        **kwargs
-):
+def metrics_fast(adata, adata_int, batch_key, label_key, **kwargs):
     """Only metrics with minimal preprocessing and runtime
 
 
@@ -58,17 +51,11 @@ def metrics_fast(
         hvg_score_=True,
         graph_conn_=True,
         pcr_=True,
-        **kwargs
+        **kwargs,
     )
 
 
-def metrics_slim(
-        adata,
-        adata_int,
-        batch_key,
-        label_key,
-        **kwargs
-):
+def metrics_slim(adata, adata_int, batch_key, label_key, **kwargs):
     """All metrics apart from kBET and LISI scores
 
     :Biological conservation:
@@ -117,17 +104,11 @@ def metrics_slim(
         nmi_=True,
         ari_=True,
         cell_cycle_=True,
-        **kwargs
+        **kwargs,
     )
 
 
-def metrics_all(
-        adata,
-        adata_int,
-        batch_key,
-        label_key,
-        **kwargs
-):
+def metrics_all(adata, adata_int, batch_key, label_key, **kwargs):
     """All metrics
 
     :Biological conservation:
@@ -184,42 +165,42 @@ def metrics_all(
         kBET_=True,
         ilisi_=True,
         clisi_=True,
-        **kwargs
+        **kwargs,
     )
 
 
 def metrics(
-        adata,
-        adata_int,
-        batch_key,
-        label_key,
-        embed='X_pca',
-        cluster_key='cluster',
-        cluster_nmi=None,
-        ari_=False,
-        nmi_=False,
-        nmi_method='arithmetic',
-        nmi_dir=None,
-        silhouette_=False,
-        si_metric='euclidean',
-        pcr_=False,
-        cell_cycle_=False,
-        organism='mouse',
-        hvg_score_=False,
-        isolated_labels_=False,  # backwards compatibility
-        isolated_labels_f1_=False,
-        isolated_labels_asw_=False,
-        n_isolated=None,
-        graph_conn_=False,
-        trajectory_=False,
-        kBET_=False,
-        lisi_graph_=False,
-        ilisi_=False,
-        clisi_=False,
-        subsample=0.5,
-        n_cores=1,
-        type_=None,
-        verbose=False,
+    adata,
+    adata_int,
+    batch_key,
+    label_key,
+    embed="X_pca",
+    cluster_key="cluster",
+    cluster_nmi=None,
+    ari_=False,
+    nmi_=False,
+    nmi_method="arithmetic",
+    nmi_dir=None,
+    silhouette_=False,
+    si_metric="euclidean",
+    pcr_=False,
+    cell_cycle_=False,
+    organism="mouse",
+    hvg_score_=False,
+    isolated_labels_=False,  # backwards compatibility
+    isolated_labels_f1_=False,
+    isolated_labels_asw_=False,
+    n_isolated=None,
+    graph_conn_=False,
+    trajectory_=False,
+    kBET_=False,
+    lisi_graph_=False,
+    ilisi_=False,
+    clisi_=False,
+    subsample=0.5,
+    n_cores=1,
+    type_=None,
+    verbose=False,
 ):
     """Master metrics function
 
@@ -318,42 +299,35 @@ def metrics(
             plot=False,
             verbose=verbose,
             inplace=True,
-            force=True
+            force=True,
         )
         if cluster_nmi is not None:
             nmi_all.to_csv(cluster_nmi, header=False)
-            print(f'saved clustering NMI values to {cluster_nmi}')
+            print(f"saved clustering NMI values to {cluster_nmi}")
 
     if nmi_:
-        print('NMI...')
+        print("NMI...")
         nmi_score = nmi(
             adata_int,
             group1=cluster_key,
             group2=label_key,
             method=nmi_method,
-            nmi_dir=nmi_dir
+            nmi_dir=nmi_dir,
         )
     else:
         nmi_score = np.nan
 
     if ari_:
-        print('ARI...')
-        ari_score = ari(
-            adata_int,
-            group1=cluster_key,
-            group2=label_key
-        )
+        print("ARI...")
+        ari_score = ari(adata_int, group1=cluster_key, group2=label_key)
     else:
         ari_score = np.nan
 
     if silhouette_:
-        print('Silhouette score...')
+        print("Silhouette score...")
         # global silhouette coefficient
         asw_label = silhouette(
-            adata_int,
-            group_key=label_key,
-            embed=embed,
-            metric=si_metric
+            adata_int, group_key=label_key, embed=embed, metric=si_metric
         )
         # silhouette coefficient per batch
         asw_batch = silhouette_batch(
@@ -363,33 +337,29 @@ def metrics(
             embed=embed,
             metric=si_metric,
             return_all=False,
-            verbose=False
+            verbose=False,
         )
     else:
         asw_label = np.nan
         asw_batch = np.nan
 
     if pcr_:
-        print('PC regression...')
+        print("PC regression...")
         pcr_score = pcr_comparison(
-            adata,
-            adata_int,
-            embed=embed,
-            covariate=batch_key,
-            verbose=verbose
+            adata, adata_int, embed=embed, covariate=batch_key, verbose=verbose
         )
     else:
         pcr_score = np.nan
 
     if cell_cycle_:
-        print('cell cycle effect...')
+        print("cell cycle effect...")
         cc_score = cell_cycle(
             adata,
             adata_int,
             batch_key=batch_key,
             embed=embed,
             agg_func=np.mean,
-            organism=organism
+            organism=organism,
         )
     else:
         cc_score = np.nan
@@ -403,36 +373,37 @@ def metrics(
             embed=embed,
             cluster=True,
             iso_threshold=n_isolated,
-            verbose=False
+            verbose=False,
         )
     else:
         il_score_f1 = np.nan
 
     if isolated_labels_asw_ or isolated_labels_:
         print("Isolated labels ASW...")
-        il_score_asw = isolated_labels(
-            adata_int,
-            label_key=label_key,
-            batch_key=batch_key,
-            embed=embed,
-            cluster=False,
-            iso_threshold=n_isolated,
-            verbose=False
-        ) if silhouette_ else np.nan
+        il_score_asw = (
+            isolated_labels(
+                adata_int,
+                label_key=label_key,
+                batch_key=batch_key,
+                embed=embed,
+                cluster=False,
+                iso_threshold=n_isolated,
+                verbose=False,
+            )
+            if silhouette_
+            else np.nan
+        )
     else:
         il_score_asw = np.nan
 
     if graph_conn_:
-        print('Graph connectivity...')
-        graph_conn_score = graph_connectivity(
-            adata_int,
-            label_key=label_key
-        )
+        print("Graph connectivity...")
+        graph_conn_score = graph_connectivity(adata_int, label_key=label_key)
     else:
         graph_conn_score = np.nan
 
     if kBET_:
-        print('kBET...')
+        print("kBET...")
         kbet_score = kBET(
             adata_int,
             batch_key=batch_key,
@@ -440,7 +411,7 @@ def metrics(
             type_=type_,
             embed=embed,
             scaled=True,
-            verbose=verbose
+            verbose=verbose,
         )
     else:
         kbet_score = np.nan
@@ -450,7 +421,7 @@ def metrics(
         ilisi_ = True
 
     if clisi_:
-        print('cLISI score...')
+        print("cLISI score...")
         clisi = clisi_graph(
             adata_int,
             batch_key=batch_key,
@@ -459,13 +430,13 @@ def metrics(
             subsample=subsample * 100,
             scale=True,
             n_cores=n_cores,
-            verbose=verbose
+            verbose=verbose,
         )
     else:
         clisi = np.nan
 
     if ilisi_:
-        print('iLISI score...')
+        print("iLISI score...")
         ilisi = ilisi_graph(
             adata_int,
             batch_key=batch_key,
@@ -473,7 +444,7 @@ def metrics(
             subsample=subsample * 100,
             scale=True,
             n_cores=n_cores,
-            verbose=verbose
+            verbose=verbose,
         )
     else:
         ilisi = np.nan
@@ -484,7 +455,7 @@ def metrics(
         hvg_score = np.nan
 
     if trajectory_:
-        print('Trajectory conservation score...')
+        print("Trajectory conservation score...")
         trajectory_score = trajectory_conservation(
             adata,
             adata_int,
@@ -495,40 +466,20 @@ def metrics(
         trajectory_score = np.nan
 
     results = {
-        'NMI_cluster/label': nmi_score,
-        'ARI_cluster/label': ari_score,
-        'ASW_label': asw_label,
-        'ASW_label/batch': asw_batch,
-        'PCR_batch': pcr_score,
-        'cell_cycle_conservation': cc_score,
-        'isolated_label_F1': il_score_f1,
-        'isolated_label_silhouette': il_score_asw,
-        'graph_conn': graph_conn_score,
-        'kBET': kbet_score,
-        'iLISI': ilisi,
-        'cLISI': clisi,
-        'hvg_overlap': hvg_score,
-        'trajectory': trajectory_score
+        "NMI_cluster/label": nmi_score,
+        "ARI_cluster/label": ari_score,
+        "ASW_label": asw_label,
+        "ASW_label/batch": asw_batch,
+        "PCR_batch": pcr_score,
+        "cell_cycle_conservation": cc_score,
+        "isolated_label_F1": il_score_f1,
+        "isolated_label_silhouette": il_score_asw,
+        "graph_conn": graph_conn_score,
+        "kBET": kbet_score,
+        "iLISI": ilisi,
+        "cLISI": clisi,
+        "hvg_overlap": hvg_score,
+        "trajectory": trajectory_score,
     }
 
-    return pd.DataFrame.from_dict(results, orient='index')
-
-
-@deprecated
-def measureTM(*args, **kwargs):
-    """
-    :param *args: function to be tested for time and memory
-    :param **kwargs: list of function parameters
-    :returns: (memory (MB), time (s), list of *args function outputs)
-    """
-    import cProfile
-    from pstats import Stats
-
-    import memory_profiler
-
-    prof = cProfile.Profile()
-    out = memory_profiler.memory_usage((prof.runcall, args, kwargs), retval=True)
-    mem = np.max(out[0]) - out[0][0]
-    print(f'memory usage:{round(mem, 0)} MB')
-    print(f'runtime: {round(Stats(prof).total_tt, 0)} s')
-    return mem, Stats(prof).total_tt, out[1:]
+    return pd.DataFrame.from_dict(results, orient="index")
diff --git a/scib/metrics/nmi.py b/scib/metrics/nmi.py
index d1c6e448..13afae2b 100644
--- a/scib/metrics/nmi.py
+++ b/scib/metrics/nmi.py
@@ -35,11 +35,11 @@ def nmi(adata, group1, group2, method="arithmetic", nmi_dir=None):
 
     if len(group1) != len(group2):
         raise ValueError(
-            f'different lengths in group1 ({len(group1)}) and group2 ({len(group2)})'
+            f"different lengths in group1 ({len(group1)}) and group2 ({len(group2)})"
         )
 
     # choose method
-    if method in ['max', 'min', 'geometric', 'arithmetic']:
+    if method in ["max", "min", "geometric", "arithmetic"]:
         nmi_value = normalized_mutual_info_score(group1, group2, average_method=method)
     elif method == "Lancichinetti":
         nmi_value = nmi_Lanc(group1, group2, nmi_dir=nmi_dir)
@@ -73,7 +73,8 @@ def onmi(group1, group2, nmi_dir=None, verbose=True):
     nmi_call = subprocess.Popen(
         [nmi_dir + "onmi", group1_file, group2_file],
         stdout=subprocess.PIPE,
-        stderr=subprocess.STDOUT)
+        stderr=subprocess.STDOUT,
+    )
 
     stdout, stderr = nmi_call.communicate()
     if stderr:
@@ -83,7 +84,7 @@ def onmi(group1, group2, nmi_dir=None, verbose=True):
     if verbose:
         print(nmi_out)
 
-    nmi_split = [x.strip().split('\t') for x in nmi_out.split('\n')]
+    nmi_split = [x.strip().split("\t") for x in nmi_out.split("\n")]
     nmi_max = float(nmi_split[0][1])
 
     # remove temporary files
@@ -102,7 +103,8 @@ def nmi_Lanc(group1, group2, nmi_dir="external/mutual3/", verbose=True):
 
     if nmi_dir is None:
         raise FileNotFoundError(
-            "Please provide the directory of the compiled C code from https://sites.google.com/site/andrealancichinetti/mutual3.tar.gz")
+            "Please provide the directory of the compiled C code from https://sites.google.com/site/andrealancichinetti/mutual3.tar.gz"
+        )
 
     group1_file = write_tmp_labels(group1, to_int=False)
     group2_file = write_tmp_labels(group2, to_int=False)
@@ -110,17 +112,18 @@ def nmi_Lanc(group1, group2, nmi_dir="external/mutual3/", verbose=True):
     nmi_call = subprocess.Popen(
         [nmi_dir + "mutual", group1_file, group2_file],
         stdout=subprocess.PIPE,
-        stderr=subprocess.STDOUT)
+        stderr=subprocess.STDOUT,
+    )
 
     stdout, stderr = nmi_call.communicate()
     if stderr:
         print(stderr)
     nmi_out = stdout.decode().strip()
 
-    return float(nmi_out.split('\t')[1])
+    return float(nmi_out.split("\t")[1])
 
 
-def write_tmp_labels(group_assignments, to_int=False, delim='\n'):
+def write_tmp_labels(group_assignments, to_int=False, delim="\n"):
     """
     write the values of a specific obs column into a temporary file in text format
     needed for external C NMI implementations (onmi and nmi_Lanc functions), because they require files as input
@@ -136,15 +139,15 @@ def write_tmp_labels(group_assignments, to_int=False, delim='\n'):
         for label in set(group_assignments):
             label_map[label] = i
             i += 1
-        labels = delim.join([str(label_map[name]) for name in group_assignments])
+        delim.join([str(label_map[name]) for name in group_assignments])
     else:
-        labels = delim.join([str(name) for name in group_assignments])
+        delim.join([str(name) for name in group_assignments])
 
     clusters = {label: [] for label in set(group_assignments)}
     for i, label in enumerate(group_assignments):
         clusters[label].append(str(i))
 
-    output = '\n'.join([' '.join(c) for c in clusters.values()])
+    output = "\n".join([" ".join(c) for c in clusters.values()])
     output = str.encode(output)
 
     # write to file
diff --git a/scib/metrics/pcr.py b/scib/metrics/pcr.py
index da7a0db3..50a57d15 100644
--- a/scib/metrics/pcr.py
+++ b/scib/metrics/pcr.py
@@ -8,13 +8,7 @@
 
 
 def pcr_comparison(
-        adata_pre,
-        adata_post,
-        covariate,
-        embed=None,
-        n_comps=50,
-        scale=True,
-        verbose=False
+    adata_pre, adata_post, covariate, embed=None, n_comps=50, scale=True, verbose=False
 ):
     """Principal component regression score
 
@@ -37,7 +31,7 @@ def pcr_comparison(
         Difference of variance contribution of PCR (scaled between 0 and 1 by default)
     """
 
-    if embed == 'X_pca':
+    if embed == "X_pca":
         embed = None
 
     pcr_before = pcr(
@@ -45,7 +39,7 @@ def pcr_comparison(
         covariate=covariate,
         recompute_pca=True,
         n_comps=n_comps,
-        verbose=verbose
+        verbose=verbose,
     )
 
     pcr_after = pcr(
@@ -54,7 +48,7 @@ def pcr_comparison(
         embed=embed,
         recompute_pca=True,
         n_comps=n_comps,
-        verbose=verbose
+        verbose=verbose,
     )
 
     if scale:
@@ -70,14 +64,7 @@ def pcr_comparison(
         return pcr_after - pcr_before
 
 
-def pcr(
-        adata,
-        covariate,
-        embed=None,
-        n_comps=50,
-        recompute_pca=True,
-        verbose=False
-):
+def pcr(adata, covariate, embed=None, n_comps=50, recompute_pca=True, verbose=False):
     """Principal component regression for anndata object
 
     Wraps :func:`~scib.metrics.pc_regression` while checking whether to:
@@ -110,10 +97,12 @@ def pcr(
         return pc_regression(adata.obsm[embed], covariate_values, n_comps=n_comps)
 
     # use existing PCA computation
-    elif (recompute_pca == False) and ('X_pca' in adata.obsm) and ('pca' in adata.uns):
+    elif (recompute_pca is False) and ("X_pca" in adata.obsm) and ("pca" in adata.uns):
         if verbose:
             print("using existing PCA")
-        return pc_regression(adata.obsm['X_pca'], covariate_values, pca_var=adata.uns['pca']['variance'])
+        return pc_regression(
+            adata.obsm["X_pca"], covariate_values, pca_var=adata.uns["pca"]["variance"]
+        )
 
     # recompute PCA
     else:
@@ -123,12 +112,7 @@ def pcr(
 
 
 def pc_regression(
-        data,
-        covariate,
-        pca_var=None,
-        n_comps=50,
-        svd_solver='arpack',
-        verbose=False
+    data, covariate, pca_var=None, n_comps=50, svd_solver="arpack", verbose=False
 ):
     """Principal component regression
 
@@ -136,7 +120,7 @@ def pc_regression(
 
     .. math::
 
-        Var(C|B) = \\sum^G_{i=1} Var(C|PC_i) \cdot R^2(PC_i|B)
+        Var(C|B) = \\sum^G_{i=1} Var(C|PC_i) \\cdot R^2(PC_i|B)
 
     for :math:`G` principal components (:math:`PC_i`), where :math:`Var(C|PC_i)` is the variance of the data matrix
     :math:`C` explained by the i-th principal component, and :math:`R^2(PC_i|B)` is the :math:`R^2` of the i-th
@@ -159,7 +143,9 @@ def pc_regression(
     if isinstance(data, (np.ndarray, sparse.csr_matrix, sparse.csc_matrix)):
         matrix = data
     else:
-        raise TypeError(f'invalid type: {data.__class__} is not a numpy array or sparse matrix')
+        raise TypeError(
+            f"invalid type: {data.__class__} is not a numpy array or sparse matrix"
+        )
 
     # perform PCA if no variance contributions are given
     if pca_var is None:
@@ -168,12 +154,18 @@ def pc_regression(
             n_comps = min(matrix.shape)
 
         if n_comps == min(matrix.shape):
-            svd_solver = 'full'
+            svd_solver = "full"
 
         if verbose:
             print("compute PCA")
-        pca = sc.tl.pca(matrix, n_comps=n_comps, use_highly_variable=False,
-                        return_info=True, svd_solver=svd_solver, copy=True)
+        pca = sc.tl.pca(
+            matrix,
+            n_comps=n_comps,
+            use_highly_variable=False,
+            return_info=True,
+            svd_solver=svd_solver,
+            copy=True,
+        )
         X_pca = pca[0].copy()
         pca_var = pca[3].copy()
         del pca
@@ -181,7 +173,7 @@ def pc_regression(
         X_pca = matrix
         n_comps = matrix.shape[1]
 
-    ## PC Regression
+    # PC Regression
     if verbose:
         print("fit regression on PCs")
 
diff --git a/scib/metrics/silhouette.py b/scib/metrics/silhouette.py
index 37a313c3..53e9d675 100644
--- a/scib/metrics/silhouette.py
+++ b/scib/metrics/silhouette.py
@@ -2,13 +2,7 @@
 from sklearn.metrics.cluster import silhouette_samples, silhouette_score
 
 
-def silhouette(
-        adata,
-        group_key,
-        embed,
-        metric='euclidean',
-        scale=True
-):
+def silhouette(adata, group_key, embed, metric="euclidean", scale=True):
     """Average silhouette width (ASW)
 
     Wrapper for sklearn silhouette function values range from [-1, 1] with
@@ -25,11 +19,9 @@ def silhouette(
     """
     if embed not in adata.obsm.keys():
         print(adata.obsm.keys())
-        raise KeyError(f'{embed} not in obsm')
+        raise KeyError(f"{embed} not in obsm")
     asw = silhouette_score(
-        X=adata.obsm[embed],
-        labels=adata.obs[group_key],
-        metric=metric
+        X=adata.obsm[embed], labels=adata.obs[group_key], metric=metric
     )
     if scale:
         asw = (asw + 1) / 2
@@ -37,14 +29,14 @@ def silhouette(
 
 
 def silhouette_batch(
-        adata,
-        batch_key,
-        group_key,
-        embed,
-        metric='euclidean',
-        return_all=False,
-        scale=True,
-        verbose=True
+    adata,
+    batch_key,
+    group_key,
+    embed,
+    metric="euclidean",
+    return_all=False,
+    scale=True,
+    verbose=True,
 ):
     """Batch ASW
 
@@ -88,9 +80,9 @@ def silhouette_batch(
     """
     if embed not in adata.obsm.keys():
         print(adata.obsm.keys())
-        raise KeyError(f'{embed} not in obsm')
+        raise KeyError(f"{embed} not in obsm")
 
-    sil_all = pd.DataFrame(columns=['group', 'silhouette_score'])
+    sil_all = pd.DataFrame(columns=["group", "silhouette_score"])
 
     for group in adata.obs[group_key].unique():
         adata_group = adata[adata.obs[group_key] == group]
@@ -100,9 +92,7 @@ def silhouette_batch(
             continue
 
         sil_per_group = silhouette_samples(
-            adata_group.obsm[embed],
-            adata_group.obs[batch_key],
-            metric=metric
+            adata_group.obsm[embed], adata_group.obs[batch_key], metric=metric
         )
 
         # take only absolute value
@@ -113,18 +103,20 @@ def silhouette_batch(
             sil_per_group = [1 - i for i in sil_per_group]
 
         sil_all = sil_all.append(
-            pd.DataFrame({
-                'group': [group] * len(sil_per_group),
-                'silhouette_score': sil_per_group
-            })
+            pd.DataFrame(
+                {
+                    "group": [group] * len(sil_per_group),
+                    "silhouette_score": sil_per_group,
+                }
+            )
         )
 
     sil_all = sil_all.reset_index(drop=True)
-    sil_means = sil_all.groupby('group').mean()
-    asw = sil_means['silhouette_score'].mean()
+    sil_means = sil_all.groupby("group").mean()
+    asw = sil_means["silhouette_score"].mean()
 
     if verbose:
-        print(f'mean silhouette per group: {sil_means}')
+        print(f"mean silhouette per group: {sil_means}")
 
     if return_all:
         return asw, sil_means, sil_all
diff --git a/scib/metrics/trajectory.py b/scib/metrics/trajectory.py
index 02ca92f9..5e13559d 100644
--- a/scib/metrics/trajectory.py
+++ b/scib/metrics/trajectory.py
@@ -8,11 +8,7 @@
 
 
 def trajectory_conservation(
-        adata_pre,
-        adata_post,
-        label_key,
-        pseudotime_key="dpt_pseudotime",
-        batch_key=None
+    adata_pre, adata_post, label_key, pseudotime_key="dpt_pseudotime", batch_key=None
 ):
     """Trajectory conservation score
 
@@ -38,52 +34,56 @@ def trajectory_conservation(
     adata_pre_ti = adata_pre[cell_subset]
     adata_post_ti = adata_post[cell_subset]
     try:
-        iroot, adata_post_ti2 = get_root(adata_pre_ti, adata_post_ti, label_key, pseudotime_key)
+        iroot, adata_post_ti2 = get_root(
+            adata_pre_ti, adata_post_ti, label_key, pseudotime_key
+        )
     except RootCellError:
-        print('No root cell found, setting trajectory conservation metric to 0.')
+        print("No root cell found, setting trajectory conservation metric to 0.")
         return 0  # failure to find root cell means no TI conservation
 
-    adata_post_ti2.uns['iroot'] = iroot
+    adata_post_ti2.uns["iroot"] = iroot
 
     sc.tl.dpt(adata_post_ti2)  # stored in 'dpt_pseudotime'
-    adata_post_ti2.obs.loc[adata_post_ti2.obs['dpt_pseudotime'] > 1, 'dpt_pseudotime'] = 0
-    adata_post_ti.obs['dpt_pseudotime'] = 0
-    adata_post_ti.obs['dpt_pseudotime'] = adata_post_ti2.obs['dpt_pseudotime']
-    adata_post_ti.obs['dpt_pseudotime'].fillna(0, inplace=True)
-
-    if batch_key == None:
+    adata_post_ti2.obs.loc[
+        adata_post_ti2.obs["dpt_pseudotime"] > 1, "dpt_pseudotime"
+    ] = 0
+    adata_post_ti.obs["dpt_pseudotime"] = 0
+    adata_post_ti.obs["dpt_pseudotime"] = adata_post_ti2.obs["dpt_pseudotime"]
+    adata_post_ti.obs["dpt_pseudotime"].fillna(0, inplace=True)
+
+    if batch_key is None:
         pseudotime_before = adata_pre_ti.obs[pseudotime_key]
-        pseudotime_after = adata_post_ti.obs['dpt_pseudotime']
-        correlation = pseudotime_before.corr(pseudotime_after, 'spearman')
+        pseudotime_after = adata_post_ti.obs["dpt_pseudotime"]
+        correlation = pseudotime_before.corr(pseudotime_after, "spearman")
         return (correlation + 1) / 2  # scaled
     else:
         check_batch(batch_key, adata_pre.obs)
         check_batch(batch_key, adata_post.obs)
 
         # check if batches match
-        if not np.array_equal(adata_post_ti.obs[batch_key], adata_pre_ti.obs[batch_key]):
+        if not np.array_equal(
+            adata_post_ti.obs[batch_key], adata_pre_ti.obs[batch_key]
+        ):
             raise ValueError(
-                'Batch columns do not match\n'
+                "Batch columns do not match\n"
                 f"adata_post_ti.obs['batch']:\n {adata_post_ti.obs[batch_key]}\n"
                 f"adata_pre_ti.obs['batch']:\n {adata_pre_ti.obs[batch_key]}\n"
             )
 
         corr = pd.Series()
         for i in adata_pre_ti.obs[batch_key].unique():
-            pseudotime_before = adata_pre_ti.obs[adata_pre_ti.obs[batch_key] == i][pseudotime_key]
-            pseudotime_after = adata_post_ti.obs[adata_post_ti.obs[batch_key] == i]['dpt_pseudotime']
-            corr[i] = pseudotime_before.corr(pseudotime_after, 'spearman')
+            pseudotime_before = adata_pre_ti.obs[adata_pre_ti.obs[batch_key] == i][
+                pseudotime_key
+            ]
+            pseudotime_after = adata_post_ti.obs[adata_post_ti.obs[batch_key] == i][
+                "dpt_pseudotime"
+            ]
+            corr[i] = pseudotime_before.corr(pseudotime_after, "spearman")
 
         return (corr.mean() + 1) / 2  # scaled
 
 
-def get_root(
-        adata_pre,
-        adata_post,
-        ct_key,
-        pseudotime_key="dpt_pseudotime",
-        dpt_dim=3
-):
+def get_root(adata_pre, adata_post, ct_key, pseudotime_key="dpt_pseudotime", dpt_dim=3):
     """Determine root cell for integrated adata based on unintegrated adata
 
     :param adata_pre: unintegrated adata
@@ -93,16 +93,19 @@ def get_root(
         Column can contain empty entries, the dataset will be subset to the cells with scores.
     :param dpt_dim: number of diffmap dimensions used to determine root
     """
-    n_components, adata_post.obs['neighborhood'] = connected_components(
-        csgraph=adata_post.obsp['connectivities'],
-        directed=False,
-        return_labels=True
+    n_components, adata_post.obs["neighborhood"] = connected_components(
+        csgraph=adata_post.obsp["connectivities"], directed=False, return_labels=True
     )
 
     start_clust = adata_pre.obs.groupby([ct_key]).mean()[pseudotime_key].idxmin()
     min_dpt = adata_pre.obs[adata_pre.obs[ct_key] == start_clust].index
-    which_max_neigh = adata_post.obs['neighborhood'] == adata_post.obs['neighborhood'].value_counts().idxmax()
-    min_dpt = [value for value in min_dpt if value in adata_post.obs[which_max_neigh].index]
+    which_max_neigh = (
+        adata_post.obs["neighborhood"]
+        == adata_post.obs["neighborhood"].value_counts().idxmax()
+    )
+    min_dpt = [
+        value for value in min_dpt if value in adata_post.obs[which_max_neigh].index
+    ]
 
     adata_post_ti = adata_post[which_max_neigh]
 
@@ -120,7 +123,7 @@ def get_root(
 
         # count opt cell
         if len(diffmap_min_dpt) == 0:
-            raise RootCellError('No root cell in largest component')
+            raise RootCellError("No root cell in largest component")
 
         # choose optimum function
         if len(diffmap_min_dpt) > 0 and diffmap_min_dpt.mean() < diffmap_mean:
diff --git a/scib/metrics/utils.py b/scib/metrics/utils.py
index a60606fb..f3af7327 100644
--- a/scib/metrics/utils.py
+++ b/scib/metrics/utils.py
@@ -4,6 +4,7 @@
 
 # Errors
 
+
 class RootCellError(Exception):
     def __init__(self, message):
         self.message = message
@@ -16,6 +17,7 @@ def __init__(self, message):
 
 # Diffusion
 
+
 def diffusion_conn(adata, min_k=50, copy=True, max_iterations=26):
     """
     Diffusion for connectivites matrix extension
@@ -31,19 +33,18 @@ def diffusion_conn(adata, min_k=50, copy=True, max_iterations=26):
        with the diffusion-enhanced connectivities matrix is in
        `adata.uns["neighbors"]["conectivities"]`
     """
-    if 'neighbors' not in adata.uns:
+    if "neighbors" not in adata.uns:
         raise ValueError(
-            '`neighbors` not in adata object. '
-            'Please compute a neighbourhood graph!'
+            "`neighbors` not in adata object. " "Please compute a neighbourhood graph!"
         )
 
-    if 'connectivities' not in adata.obsp:
+    if "connectivities" not in adata.obsp:
         raise ValueError(
-            '`connectivities` not in `adata.obsp`. '
-            'Please pass an object with connectivities computed!'
+            "`connectivities` not in `adata.obsp`. "
+            "Please pass an object with connectivities computed!"
         )
 
-    T = adata.obsp['connectivities']
+    T = adata.obsp["connectivities"]
 
     # Normalize T with max row sum
     # Note: This keeps the matrix symmetric and ensures |M| doesn't keep growing
@@ -52,8 +53,9 @@ def diffusion_conn(adata, min_k=50, copy=True, max_iterations=26):
     M = T
 
     # Check for disconnected component
-    n_comp, labs = sparse.csgraph.connected_components(adata.obsp['connectivities'],
-                                                       connection='strong')
+    n_comp, labs = sparse.csgraph.connected_components(
+        adata.obsp["connectivities"], connection="strong"
+    )
 
     if n_comp > 1:
         tab = pd.value_counts(labs)
@@ -64,25 +66,27 @@ def diffusion_conn(adata, min_k=50, copy=True, max_iterations=26):
 
     T_agg = T
     i = 2
-    while ((M[large_comp_mask, :][:, large_comp_mask] > 0).sum(1).min() < min_k) and (i < max_iterations):
-        print(f'Adding diffusion to step {i}')
+    while ((M[large_comp_mask, :][:, large_comp_mask] > 0).sum(1).min() < min_k) and (
+        i < max_iterations
+    ):
+        print(f"Adding diffusion to step {i}")
         T_agg *= T
         M += T_agg
         i += 1
 
     if (M[large_comp_mask, :][:, large_comp_mask] > 0).sum(1).min() < min_k:
         raise ValueError(
-            'could not create diffusion connectivities matrix'
-            f'with at least {min_k} non-zero entries in'
-            f'{max_iterations} iterations.\n Please increase the'
-            'value of max_iterations or reduce k_min.\n'
+            "could not create diffusion connectivities matrix"
+            f"with at least {min_k} non-zero entries in"
+            f"{max_iterations} iterations.\n Please increase the"
+            "value of max_iterations or reduce k_min.\n"
         )
 
     M.setdiag(0)
 
     if copy:
         adata_tmp = adata.copy()
-        adata_tmp.uns['neighbors'].update({'diffusion_connectivities': M})
+        adata_tmp.uns["neighbors"].update({"diffusion_connectivities": M})
         return adata_tmp
 
     else:
@@ -99,34 +103,39 @@ def diffusion_nn(adata, k, max_iterations=26):
     Return:
        `k_indices` a numpy.ndarray of the indices of the k-nearest neighbors.
     """
-    if 'neighbors' not in adata.uns:
-        raise ValueError('`neighbors` not in adata object. '
-                         'Please compute a neighbourhood graph!')
+    if "neighbors" not in adata.uns:
+        raise ValueError(
+            "`neighbors` not in adata object. " "Please compute a neighbourhood graph!"
+        )
 
-    if 'connectivities' not in adata.obsp:
-        raise ValueError('`connectivities` not in `adata.obsp`. '
-                         'Please pass an object with connectivities computed!')
+    if "connectivities" not in adata.obsp:
+        raise ValueError(
+            "`connectivities` not in `adata.obsp`. "
+            "Please pass an object with connectivities computed!"
+        )
 
-    T = adata.obsp['connectivities']
+    T = adata.obsp["connectivities"]
 
     # Row-normalize T
     T = sparse.diags(1 / T.sum(1).A.ravel()) * T
 
-    T_agg = T ** 3
-    M = T + T ** 2 + T_agg
+    T_agg = T**3
+    M = T + T**2 + T_agg
     i = 4
 
     while ((M > 0).sum(1).min() < (k + 1)) and (i < max_iterations):
         # note: k+1 is used as diag is non-zero (self-loops)
-        print(f'Adding diffusion to step {i}')
+        print(f"Adding diffusion to step {i}")
         T_agg *= T
         M += T_agg
         i += 1
 
     if (M > 0).sum(1).min() < (k + 1):
-        raise NeighborsError(f'could not find {k} nearest neighbors in {max_iterations}'
-                             'diffusion steps.\n Please increase max_iterations or reduce'
-                             ' k.\n')
+        raise NeighborsError(
+            f"could not find {k} nearest neighbors in {max_iterations}"
+            "diffusion steps.\n Please increase max_iterations or reduce"
+            " k.\n"
+        )
 
     M.setdiag(0)
     k_indices = np.argpartition(M.A, -k, axis=1)[:, -k:]
@@ -136,16 +145,19 @@ def diffusion_nn(adata, k, max_iterations=26):
 
 # Not used
 
+
 def get_hvg_indices(adata, verbose=True):
     if "highly_variable" not in adata.var.columns:
         if verbose:
-            print(f"No highly variable genes computed, continuing with full matrix {adata.shape}")
+            print(
+                f"No highly variable genes computed, continuing with full matrix {adata.shape}"
+            )
         return np.array(range(adata.n_vars))
-    return np.where((adata.var["highly_variable"] == True))[0]
+    return np.where(adata.var["highly_variable"] is True)[0]
 
 
 def select_hvg(adata, select=True):
-    if select and 'highly_variable' in adata.var:
-        return adata[:, adata.var['highly_variable']].copy()
+    if select and "highly_variable" in adata.var:
+        return adata[:, adata.var["highly_variable"]].copy()
     else:
         return adata
diff --git a/scib/preprocessing.py b/scib/preprocessing.py
index ea68a6d3..4c46d24c 100644
--- a/scib/preprocessing.py
+++ b/scib/preprocessing.py
@@ -1,12 +1,10 @@
-# TODO: move util functions e.g. reader functions elsewhere
-
 import logging
 import tempfile
 
 import anndata2ri
 import numpy as np
-# rpy2 for running R code
-import rpy2.rinterface_lib.callbacks
+import rpy2.rinterface_lib.callbacks  # rpy2 for running R code
+import rpy2.rinterface_lib.embedded
 import rpy2.robjects as ro
 import scanpy as sc
 import seaborn
@@ -14,15 +12,15 @@
 from matplotlib import pyplot as plt
 from scipy import sparse
 
-# access to other methods of this module
-from . import utils
+from . import utils  # TODO: move util fcns (eg reader) elsewhere
 from .exceptions import RLibraryNotFound
 
-rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)  # Ignore R warning messages
-seaborn.set_context('talk')
+# Ignore R warning messages
+rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)
+seaborn.set_context("talk")
 
 
-def summarize_counts(adata, count_matrix=None, mt_gene_regex='^MT-'):
+def summarize_counts(adata, count_matrix=None, mt_gene_regex="^MT-"):
     """Summarise counts of the given count matrix
 
     This function is useful for quality control.
@@ -40,11 +38,11 @@ def summarize_counts(adata, count_matrix=None, mt_gene_regex='^MT-'):
 
     if count_matrix is None:
         count_matrix = adata.X
-    adata.obs['n_counts'] = count_matrix.sum(1)
-    adata.obs['log_counts'] = np.log(adata.obs['n_counts'])
-    adata.obs['n_genes'] = (count_matrix > 0).sum(1)
+    adata.obs["n_counts"] = count_matrix.sum(1)
+    adata.obs["log_counts"] = np.log(adata.obs["n_counts"])
+    adata.obs["n_genes"] = (count_matrix > 0).sum(1)
 
-    if mt_gene_regex != None:
+    if mt_gene_regex is not None:
         # for each cell compute fraction of counts in mito genes vs. all genes
         mito_genes = adata.var_names.str.match(mt_gene_regex)
         mt_sum = np.sum(adata[:, mito_genes].X, axis=1)
@@ -53,7 +51,7 @@ def summarize_counts(adata, count_matrix=None, mt_gene_regex='^MT-'):
         if sparse.issparse(adata.X):
             mt_sum = mt_sum.A1
             total_sum = total_sum.A1
-        adata.obs['mito_frac'] = mt_sum / total_sum
+        adata.obs["mito_frac"] = mt_sum / total_sum
 
         # mt_gene_mask = [gene.startswith('mt-') for gene in adata.var_names]
         # mt_count = count_matrix[:, mt_gene_mask].sum(1)
@@ -62,17 +60,17 @@ def summarize_counts(adata, count_matrix=None, mt_gene_regex='^MT-'):
         # adata.obs['mt_frac'] = mt_count/adata.obs['n_counts']
 
 
-### Quality Control
+# Quality Control
 def plot_qc(
-        adata,
-        color=None,
-        bins=60,
-        legend_loc='right margin',
-        histogram=True,
-        count_threshold=(0, np.inf),
-        count_filter_threshold=(0, np.inf),
-        gene_threshold=(0, np.inf),
-        gene_filter_threshold=(0, np.inf)
+    adata,
+    color=None,
+    bins=60,
+    legend_loc="right margin",
+    histogram=True,
+    count_threshold=(0, np.inf),
+    count_filter_threshold=(0, np.inf),
+    gene_threshold=(0, np.inf),
+    gene_filter_threshold=(0, np.inf),
 ):
     """Create QC Plots
 
@@ -105,7 +103,7 @@ def plot_qc(
         title=color,
         gene_threshold=gene_filter_threshold[0],
         count_threshold=count_filter_threshold[0],
-        legend_loc=legend_loc
+        legend_loc=legend_loc,
     )
 
     if not histogram:
@@ -116,12 +114,12 @@ def plot_qc(
         # count filtering
         plot_count_filter(
             adata,
-            obs_col='n_counts',
+            obs_col="n_counts",
             bins=bins,
             lower=count_threshold[0],
             filter_lower=count_filter_threshold[0],
             upper=count_threshold[1],
-            filter_upper=count_filter_threshold[1]
+            filter_upper=count_filter_threshold[1],
         )
 
     if gene_filter_threshold != (0, np.inf):
@@ -129,85 +127,101 @@ def plot_qc(
         # gene filtering
         plot_count_filter(
             adata,
-            obs_col='n_genes',
+            obs_col="n_genes",
             bins=bins,
             lower=gene_threshold[0],
             filter_lower=gene_filter_threshold[0],
             upper=gene_threshold[1],
-            filter_upper=gene_filter_threshold[1]
+            filter_upper=gene_filter_threshold[1],
         )
 
 
 def plot_scatter(
-        adata,
-        count_threshold=0,
-        gene_threshold=0,
-        color=None,
-        title='Scatter plot of count depth vs gene counts',
-        lab_size=15,
-        tick_size=11,
-        legend_loc='right margin',
-        palette=None
+    adata,
+    count_threshold=0,
+    gene_threshold=0,
+    color=None,
+    title="Scatter plot of count depth vs gene counts",
+    lab_size=15,
+    tick_size=11,
+    legend_loc="right margin",
+    palette=None,
 ):
     utils.check_adata(adata)
     if color:
         utils.check_batch(color, adata.obs)
 
-    ax = sc.pl.scatter(adata, 'n_counts', 'n_genes', color=color, show=False,
-                       legend_fontweight=50, legend_loc=legend_loc, palette=palette)
+    ax = sc.pl.scatter(
+        adata,
+        "n_counts",
+        "n_genes",
+        color=color,
+        show=False,
+        legend_fontweight=50,
+        legend_loc=legend_loc,
+        palette=palette,
+    )
     ax.set_title(title, fontsize=lab_size)
     ax.set_xlabel("Count depth", fontsize=lab_size)
     ax.set_ylabel("Number of genes", fontsize=lab_size)
     ax.tick_params(labelsize=tick_size)
 
     if gene_threshold > 0:
-        ax.axhline(gene_threshold, 0, 1, color='red')
+        ax.axhline(gene_threshold, 0, 1, color="red")
     if count_threshold > 0:
-        ax.axvline(count_threshold, 0, 1, color='red')
+        ax.axvline(count_threshold, 0, 1, color="red")
 
     fig = plt.gcf()
     cbar_ax = fig.axes[-1]
     cbar_ax.tick_params(labelsize=tick_size)
-    f1 = ax.get_figure()
+    ax.get_figure()
     plt.show()
 
 
-def plot_count_filter(adata, obs_col='n_counts', bins=60, lower=0, upper=np.inf, filter_lower=0, filter_upper=np.inf):
+def plot_count_filter(
+    adata,
+    obs_col="n_counts",
+    bins=60,
+    lower=0,
+    upper=np.inf,
+    filter_lower=0,
+    filter_upper=np.inf,
+):
     plot_data = adata.obs[obs_col]
 
     sns.distplot(plot_data, kde=False, bins=bins)
 
     if lower > 0:
-        plt.axvline(lower, linestyle='--', color='g')
+        plt.axvline(lower, linestyle="--", color="g")
     if filter_lower > 0:
-        plt.axvline(filter_lower, linestyle='-', color='r')
+        plt.axvline(filter_lower, linestyle="-", color="r")
     if not np.isinf(upper):
-        plt.axvline(upper, linestyle='--', color='g')
+        plt.axvline(upper, linestyle="--", color="g")
     if not np.isinf(upper):
-        plt.axvline(filter_upper, linestyle='-', color='r')
-    plt.ticklabel_format(style='sci', axis='x', scilimits=(0, 2))
+        plt.axvline(filter_upper, linestyle="-", color="r")
+    plt.ticklabel_format(style="sci", axis="x", scilimits=(0, 2))
     plt.show()
 
     # determine lower bound of total, look at points below lower bound
     if filter_lower > 0:
         print(f"lower threshold: {filter_lower}")
         sns.distplot(plot_data[plot_data < lower], kde=False, bins=bins)
-        plt.axvline(filter_lower, linestyle='-', color='r')
-        plt.axvline(lower, linestyle='--', color='g')
-        plt.ticklabel_format(style='sci', axis='x', scilimits=(0, 2))
+        plt.axvline(filter_lower, linestyle="-", color="r")
+        plt.axvline(lower, linestyle="--", color="g")
+        plt.ticklabel_format(style="sci", axis="x", scilimits=(0, 2))
         plt.show()
 
     # determine upper bound of total
     if not np.isinf(filter_upper) and not np.isinf(upper):
         print(f"upper threshold: {filter_upper}")
         sns.distplot(plot_data[plot_data > upper], kde=False, bins=bins)
-        plt.axvline(filter_upper, linestyle='-', color='r')
-        plt.axvline(upper, linestyle='--', color='g')
-        plt.ticklabel_format(style='sci', axis='x', scilimits=(0, 2))
+        plt.axvline(filter_upper, linestyle="-", color="r")
+        plt.axvline(upper, linestyle="--", color="g")
+        plt.ticklabel_format(style="sci", axis="x", scilimits=(0, 2))
         plt.show()
 
 
-### Normalisation
+# Normalisation
 def normalize(adata, sparsify=True, precluster=True, min_mean=0.1, log=True):
     """Normalise counts using the ``scran`` normalisation method
 
@@ -226,13 +240,17 @@ def normalize(adata, sparsify=True, precluster=True, min_mean=0.1, log=True):
 
     # Check for 0 count cells
     if np.any(adata.X.sum(axis=1) == 0):
-        raise ValueError('found 0 count cells in the AnnData object.'
-                         ' Please filter these from your dataset.')
+        raise ValueError(
+            "found 0 count cells in the AnnData object."
+            " Please filter these from your dataset."
+        )
 
     # Check for 0 count genes
     if np.any(adata.X.sum(axis=0) == 0):
-        raise ValueError('found 0 count genes in the AnnData object.'
-                         ' Please filter these from your dataset.')
+        raise ValueError(
+            "found 0 count genes in the AnnData object."
+            " Please filter these from your dataset."
+        )
 
     if sparsify:
         # massive speedup when working with sparse matrix
@@ -240,8 +258,8 @@ def normalize(adata, sparsify=True, precluster=True, min_mean=0.1, log=True):
             adata.X = sparse.csr_matrix(adata.X)
 
     try:
-        ro.r('library(scran)')
-    except Exception as ex:
+        ro.r("library(scran)")
+    except rpy2.rinterface_lib.embedded.RRuntimeError as ex:
         RLibraryNotFound(ex)
 
     anndata2ri.activate()
@@ -255,42 +273,44 @@ def normalize(adata, sparsify=True, precluster=True, min_mean=0.1, log=True):
     if sparse.issparse(X):
         is_sparse = True
 
-        if X.nnz > 2 ** 31 - 1:
+        if X.nnz > np.power(2, 31) - 1:
             X = X.tocoo()
         else:
             X = X.tocsc()
 
-    ro.globalenv['data_mat'] = X
+    ro.globalenv["data_mat"] = X
 
     if precluster:
         # Preliminary clustering for differentiated normalisation
         adata_pp = adata.copy()
         sc.pp.normalize_per_cell(adata_pp, counts_per_cell_after=1e6)
         sc.pp.log1p(adata_pp)
-        sc.pp.pca(adata_pp, n_comps=15, svd_solver='arpack')
+        sc.pp.pca(adata_pp, n_comps=15, svd_solver="arpack")
         sc.pp.neighbors(adata_pp)
-        sc.tl.louvain(adata_pp, key_added='groups', resolution=0.5)
+        sc.tl.louvain(adata_pp, key_added="groups", resolution=0.5)
 
-        ro.globalenv['input_groups'] = adata_pp.obs['groups']
+        ro.globalenv["input_groups"] = adata_pp.obs["groups"]
         size_factors = ro.r(
-            'sizeFactors('
-            '   computeSumFactors('
-            '       SingleCellExperiment(list(counts=data_mat)),'
-            '       clusters = input_groups,'
-            f'       min.mean = {min_mean}'
-            '   )'
-            ')'
+            "sizeFactors("
+            "   computeSumFactors("
+            "       SingleCellExperiment(list(counts=data_mat)),"
+            "       clusters = input_groups,"
+            f"       min.mean = {min_mean}"
+            "   )"
+            ")"
         )
 
         del adata_pp
 
     else:
-        size_factors = ro.r('sizeFactors(computeSumFactors(SingleCellExperiment('
-                            f'list(counts=data_mat)), min.mean = {min_mean}))')
+        size_factors = ro.r(
+            "sizeFactors(computeSumFactors(SingleCellExperiment("
+            f"list(counts=data_mat)), min.mean = {min_mean}))"
+        )
 
     # modify adata
-    adata.obs['size_factors'] = size_factors
-    adata.X /= adata.obs['size_factors'].values[:, None]
+    adata.obs["size_factors"] = size_factors
+    adata.X /= adata.obs["size_factors"].values[:, None]
     if log:
         print("Note! Performing log1p-transformation after normalization.")
         sc.pp.log1p(adata)
@@ -304,11 +324,13 @@ def normalize(adata, sparsify=True, precluster=True, min_mean=0.1, log=True):
     adata.raw = adata  # Store the full data set in 'raw' as log-normalised data for statistical testing
 
     # Free memory in R
-    ro.r('rm(list=ls())')
-    ro.r('lapply(names(sessionInfo()$loadedOnly), require, character.only = TRUE)')
-    ro.r('invisible(lapply(paste0("package:", names(sessionInfo()$otherPkgs)), '
-         'detach, character.only=TRUE, unload=TRUE))')
-    ro.r('gc()')
+    ro.r("rm(list=ls())")
+    ro.r("lapply(names(sessionInfo()$loadedOnly), require, character.only = TRUE)")
+    ro.r(
+        'invisible(lapply(paste0("package:", names(sessionInfo()$otherPkgs)), '
+        "detach, character.only=TRUE, unload=TRUE))"
+    )
+    ro.r("gc()")
 
     anndata2ri.deactivate()
 
@@ -353,15 +375,15 @@ def scale_batch(adata, batch):
 
 
 def hvg_intersect(
-        adata,
-        batch,
-        target_genes=2000,
-        flavor='cell_ranger',
-        n_bins=20,
-        adataOut=False,
-        n_stop=8000,
-        min_genes=500,
-        step_size=1000
+    adata,
+    batch,
+    target_genes=2000,
+    flavor="cell_ranger",
+    n_bins=20,
+    adataOut=False,
+    n_stop=8000,
+    min_genes=500,
+    step_size=1000,
 ):
     """Highly variable gene selection
 
@@ -386,14 +408,12 @@ def hvg_intersect(
     hvg_res = []
 
     for i in split:
-        sc.pp.filter_genes(i, min_cells=1)  # remove genes unexpressed (otherwise hvg might break)
+        sc.pp.filter_genes(
+            i, min_cells=1
+        )  # remove genes unexpressed (otherwise hvg might break)
         hvg_res.append(
             sc.pp.highly_variable_genes(
-                i,
-                flavor=flavor,
-                n_top_genes=n_hvg,
-                n_bins=n_bins,
-                inplace=False
+                i, flavor=flavor, n_top_genes=n_hvg, n_bins=n_bins, inplace=False
             )
         )
 
@@ -401,11 +421,11 @@ def hvg_intersect(
         genes = []
 
         for i in range(len(split)):
-            dispersion_norm = hvg_res[i]['dispersions_norm']
+            dispersion_norm = hvg_res[i]["dispersions_norm"]
             dispersion_norm = dispersion_norm[~np.isnan(dispersion_norm)]
             dispersion_norm[::-1].sort()
             disp_cut_off = dispersion_norm[n_hvg - 1]
-            gene_subset = np.nan_to_num(hvg_res[i]['dispersions_norm']) >= disp_cut_off
+            gene_subset = np.nan_to_num(hvg_res[i]["dispersions_norm"]) >= disp_cut_off
 
             genes.append(set(split[i].var[gene_subset].index))
 
@@ -415,9 +435,11 @@ def hvg_intersect(
         else:
             if n_hvg > n_stop:
                 if len(intersect) < min_genes:
-                    raise Exception(f'Only {len(intersect)} HVGs were found in the intersection.\n'
-                                    f'This is fewer than {min_genes} HVGs set as the minimum.\n'
-                                    'Consider raising `n_stop` or reducing `min_genes`.')
+                    raise Exception(
+                        f"Only {len(intersect)} HVGs were found in the intersection.\n"
+                        f"This is fewer than {min_genes} HVGs set as the minimum.\n"
+                        "Consider raising `n_stop` or reducing `min_genes`."
+                    )
                 break
             n_hvg = int(n_hvg + step_size)
 
@@ -427,13 +449,20 @@ def hvg_intersect(
     return list(intersect)
 
 
-def hvg_batch(adata, batch_key=None, target_genes=2000, flavor='cell_ranger', n_bins=20, adataOut=False):
+def hvg_batch(
+    adata,
+    batch_key=None,
+    target_genes=2000,
+    flavor="cell_ranger",
+    n_bins=20,
+    adataOut=False,
+):
     """Batch-aware highly variable gene selection
 
     Method to select HVGs based on mean dispersions of genes that are highly
     variable genes in all batches. Using a the top target_genes per batch by
-    average normalize dispersion. If target genes still hasn't been reached, 
-    then HVGs in all but one batches are used to fill up. This is continued 
+    average normalize dispersion. If target genes still hasn't been reached,
+    then HVGs in all but one batches are used to fill up. This is continued
     until HVGs in a single batch are considered.
 
     :param adata: ``anndata`` object
@@ -458,12 +487,13 @@ def hvg_batch(adata, batch_key=None, target_genes=2000, flavor='cell_ranger', n_
         flavor=flavor,
         n_top_genes=target_genes,
         n_bins=n_bins,
-        batch_key=batch_key
+        batch_key=batch_key,
     )
 
-    nbatch1_dispersions = adata_hvg.var['dispersions_norm'][
-        adata_hvg.var.highly_variable_nbatches > len(adata_hvg.obs[batch_key].cat.categories) - 1
-        ]
+    nbatch1_dispersions = adata_hvg.var["dispersions_norm"][
+        adata_hvg.var.highly_variable_nbatches
+        > len(adata_hvg.obs[batch_key].cat.categories) - 1
+    ]
 
     nbatch1_dispersions.sort_values(ascending=False, inplace=True)
 
@@ -472,28 +502,33 @@ def hvg_batch(adata, batch_key=None, target_genes=2000, flavor='cell_ranger', n_
 
     else:
         enough = False
-        print(f'Using {len(nbatch1_dispersions)} HVGs from full intersect set')
+        print(f"Using {len(nbatch1_dispersions)} HVGs from full intersect set")
         hvg = nbatch1_dispersions.index[:]
         not_n_batches = 1
 
         while not enough:
             target_genes_diff = target_genes - len(hvg)
 
-            tmp_dispersions = adata_hvg.var['dispersions_norm'][adata_hvg.var.highly_variable_nbatches ==
-                                                                (n_batches - not_n_batches)]
+            tmp_dispersions = adata_hvg.var["dispersions_norm"][
+                adata_hvg.var.highly_variable_nbatches == (n_batches - not_n_batches)
+            ]
 
             if len(tmp_dispersions) < target_genes_diff:
-                print(f'Using {len(tmp_dispersions)} HVGs from n_batch-{not_n_batches} set')
+                print(
+                    f"Using {len(tmp_dispersions)} HVGs from n_batch-{not_n_batches} set"
+                )
                 hvg = hvg.append(tmp_dispersions.index)
                 not_n_batches += 1
 
             else:
-                print(f'Using {target_genes_diff} HVGs from n_batch-{not_n_batches} set')
+                print(
+                    f"Using {target_genes_diff} HVGs from n_batch-{not_n_batches} set"
+                )
                 tmp_dispersions.sort_values(ascending=False, inplace=True)
                 hvg = hvg.append(tmp_dispersions.index[:target_genes_diff])
                 enough = True
 
-    print(f'Using {len(hvg)} HVGs')
+    print(f"Using {len(hvg)} HVGs")
 
     if not adataOut:
         del adata_hvg
@@ -502,19 +537,19 @@ def hvg_batch(adata, batch_key=None, target_genes=2000, flavor='cell_ranger', n_
         return adata_hvg[:, hvg].copy()
 
 
-### Feature Reduction
+# Feature Reduction
 def reduce_data(
-        adata,
-        batch_key=None,
-        flavor='cell_ranger',
-        n_top_genes=2000,
-        n_bins=20,
-        pca=True,
-        pca_comps=50,
-        overwrite_hvg=True,
-        neighbors=True,
-        use_rep='X_pca',
-        umap=True
+    adata,
+    batch_key=None,
+    flavor="cell_ranger",
+    n_top_genes=2000,
+    n_bins=20,
+    pca=True,
+    pca_comps=50,
+    overwrite_hvg=True,
+    neighbors=True,
+    use_rep="X_pca",
+    umap=True,
 ):
     """Apply feature selection and dimensionality reduction steps.
 
@@ -548,34 +583,35 @@ def reduce_data(
 
         overwrite_hvg = False
 
-        ## quick fix: HVG doesn't work on dense matrix
+        # quick fix: HVG doesn't work on dense matrix
         if not sparse.issparse(adata.X):
             adata.X = sparse.csr_matrix(adata.X)
 
         if batch_key is not None:
-            hvg_list = hvg_batch(adata, batch_key=batch_key, target_genes=n_top_genes, n_bins=n_bins)
-            adata.var['highly_variable'] = np.in1d(adata.var_names, hvg_list)
+            hvg_list = hvg_batch(
+                adata, batch_key=batch_key, target_genes=n_top_genes, n_bins=n_bins
+            )
+            adata.var["highly_variable"] = np.in1d(adata.var_names, hvg_list)
 
         else:
             print(f"Calculating {n_top_genes} HVGs for reduce_data.")
             sc.pp.highly_variable_genes(
-                adata,
-                n_top_genes=n_top_genes,
-                n_bins=n_bins,
-                flavor=flavor
+                adata, n_top_genes=n_top_genes, n_bins=n_bins, flavor=flavor
             )
 
         n_hvg = np.sum(adata.var["highly_variable"])
-        print(f'Computed {n_hvg} highly variable genes')
+        print(f"Computed {n_hvg} highly variable genes")
 
     if pca:
         print("PCA")
         use_hvgs = not overwrite_hvg and "highly_variable" in adata.var
-        sc.tl.pca(adata,
-                  n_comps=pca_comps,
-                  use_highly_variable=use_hvgs,
-                  svd_solver='arpack',
-                  return_info=True)
+        sc.tl.pca(
+            adata,
+            n_comps=pca_comps,
+            use_highly_variable=use_hvgs,
+            svd_solver="arpack",
+            return_info=True,
+        )
 
     if neighbors:
         print("Nearest Neigbours")
@@ -586,8 +622,8 @@ def reduce_data(
         sc.tl.umap(adata)
 
 
-### Cell Cycle
-def score_cell_cycle(adata, organism='mouse'):
+# Cell Cycle
+def score_cell_cycle(adata, organism="mouse"):
     """Score cell cycle score given an organism
 
     Wrapper function for `scanpy.tl.score_genes_cell_cycle`_
@@ -605,22 +641,31 @@ def score_cell_cycle(adata, organism='mouse'):
     :return: tuple of ``(s_genes, g2m_genes)`` of S-phase genes and G2- and M-phase genes scores
     """
     import pathlib
-    root = pathlib.Path(__file__).parent
 
-    cc_files = {'mouse': [root / 'resources/s_genes_tirosh.txt',
-                          root / 'resources/g2m_genes_tirosh.txt'],
-                'human': [root / 'resources/s_genes_tirosh_hm.txt',
-                          root / 'resources/g2m_genes_tirosh_hm.txt']}
+    root = pathlib.Path(__file__).parent
 
-    with open(cc_files[organism][0], "r") as f:
+    cc_files = {
+        "mouse": [
+            root / "resources/s_genes_tirosh.txt",
+            root / "resources/g2m_genes_tirosh.txt",
+        ],
+        "human": [
+            root / "resources/s_genes_tirosh_hm.txt",
+            root / "resources/g2m_genes_tirosh_hm.txt",
+        ],
+    }
+
+    with open(cc_files[organism][0]) as f:
         s_genes = [x.strip() for x in f.readlines() if x.strip() in adata.var.index]
-    with open(cc_files[organism][1], "r") as f:
+    with open(cc_files[organism][1]) as f:
         g2m_genes = [x.strip() for x in f.readlines() if x.strip() in adata.var.index]
 
     if (len(s_genes) == 0) or (len(g2m_genes) == 0):
         rand_choice = np.random.randint(1, adata.n_vars, 10)
         rand_genes = adata.var_names[rand_choice].tolist()
-        raise ValueError(f"cell cycle genes not in adata\n organism: {organism}\n varnames: {rand_genes}")
+        raise ValueError(
+            f"cell cycle genes not in adata\n organism: {organism}\n varnames: {rand_genes}"
+        )
 
     sc.tl.score_genes_cell_cycle(adata, s_genes, g2m_genes)
 
@@ -641,9 +686,9 @@ def save_seurat(adata, path, batch, hvgs=None):
     import re
 
     try:
-        ro.r('library(Seurat)')
-        ro.r('library(scater)')
-    except Exception as ex:
+        ro.r("library(Seurat)")
+        ro.r("library(scater)")
+    except rpy2.rinterface_lib.embedded.RRuntimeError as ex:
         RLibraryNotFound(ex)
 
     anndata2ri.activate()
@@ -657,7 +702,7 @@ def save_seurat(adata, path, batch, hvgs=None):
             if not adata.layers[key].has_sorted_indices:
                 adata.layers[key].sort_indices()
 
-    ro.globalenv['adata'] = adata
+    ro.globalenv["adata"] = adata
 
     ro.r('sobj = as.Seurat(adata, counts="counts", data = "X")')
 
@@ -666,10 +711,10 @@ def save_seurat(adata, path, batch, hvgs=None):
     ro.r(f'Idents(sobj) = "{batch}"')
     ro.r(f'saveRDS(sobj, file="{path}")')
     if hvgs is not None:
-        hvg_out = re.sub('\.RDS$', '', path) + '_hvg.RDS'
+        hvg_out = re.sub(r"\.RDS$", "", path) + "_hvg.RDS"
         # hvg_out = path+'_hvg.rds'
-        ro.globalenv['hvgs'] = hvgs
-        ro.r('unlist(hvgs)')
+        ro.globalenv["hvgs"] = hvgs
+        ro.r("unlist(hvgs)")
         ro.r(f'saveRDS(hvgs, file="{hvg_out}")')
 
     anndata2ri.deactivate()
@@ -684,26 +729,26 @@ def read_seurat(path):
     """
 
     try:
-        ro.r('library(Seurat)')
-        ro.r('library(scater)')
-    except Exception as ex:
+        ro.r("library(Seurat)")
+        ro.r("library(scater)")
+    except rpy2.rinterface_lib.embedded.RRuntimeError as ex:
         RLibraryNotFound(ex)
 
     anndata2ri.activate()
 
     ro.r(f'sobj <- readRDS("{path}")')
-    adata = ro.r('as.SingleCellExperiment(sobj)')
+    adata = ro.r("as.SingleCellExperiment(sobj)")
 
     anndata2ri.deactivate()
 
     # Test for 'X_EMB'
-    if 'X_EMB' in adata.obsm:
-        if 'X_emb' in adata.obsm:
+    if "X_EMB" in adata.obsm:
+        if "X_emb" in adata.obsm:
             print('overwriting existing `adata.obsm["X_emb"] in the adata object')
-        adata.obsm['X_emb'] = adata.obsm['X_EMB']
-        del adata.obsm['X_EMB']
+        adata.obsm["X_emb"] = adata.obsm["X_EMB"]
+        del adata.obsm["X_EMB"]
 
-    return (adata)
+    return adata
 
 
 def read_conos(inPath, dir_path=None):
@@ -715,25 +760,27 @@ def read_conos(inPath, dir_path=None):
     :param dir_path:
     """
     from shutil import rmtree
+
     import pandas as pd
     from scipy.io import mmread
 
     if dir_path is None:
         tmpdir = tempfile.TemporaryDirectory()
-        dir_path = tmpdir.name + '/'
+        dir_path = tmpdir.name + "/"
 
     try:
-        ro.r('library(conos)')
-        ro.r('library(data.table)')
-    except Exception as ex:
+        ro.r("library(conos)")
+        ro.r("library(data.table)")
+    except rpy2.rinterface_lib.embedded.RRuntimeError as ex:
         RLibraryNotFound(ex)
 
-
     ro.r(f'con <- readRDS("{inPath}")')
-    ro.r('meta <- function(sobj) {return(sobj@meta.data)}')
-    ro.r('metalist <- lapply(con$samples, meta)')
-    ro.r('metaM <- do.call(rbind,unname(metalist))')
-    ro.r(f'saveConosForScanPy(con, output.path="{dir_path}", pseudo.pca=TRUE, pca=TRUE, metadata.df=metaM)')
+    ro.r("meta <- function(sobj) {return(sobj@meta.data)}")
+    ro.r("metalist <- lapply(con$samples, meta)")
+    ro.r("metaM <- do.call(rbind,unname(metalist))")
+    ro.r(
+        f'saveConosForScanPy(con, output.path="{dir_path}", pseudo.pca=TRUE, pca=TRUE, metadata.df=metaM)'
+    )
 
     gene_df = pd.read_csv(dir_path + "genes.csv")
 
@@ -758,13 +805,15 @@ def read_conos(inPath, dir_path=None):
 
     # Depends on which PCA you loaded
     adata.X_pca = pseudopca_df.values
-    adata.obsm['X_pca'] = pseudopca_df.values
+    adata.obsm["X_pca"] = pseudopca_df.values
 
     # Name according to embedding you saved
     adata.X_umap = embedding_df.values
-    adata.obsm['X_umap'] = embedding_df.values
+    adata.obsm["X_umap"] = embedding_df.values
 
-    adata.uns['neighbors'] = dict(connectivities=graph_conn_mtx.tocsr(), distances=graph_dist_mtx.tocsr())
+    adata.uns["neighbors"] = dict(
+        connectivities=graph_conn_mtx.tocsr(), distances=graph_dist_mtx.tocsr()
+    )
 
     # Assign raw counts to .raw slot, load in normalised counts
     # adata.raw = adata
diff --git a/scib/resources/g2m_genes_tirosh.txt b/scib/resources/g2m_genes_tirosh.txt
index a2fe82fd..d191b986 100644
--- a/scib/resources/g2m_genes_tirosh.txt
+++ b/scib/resources/g2m_genes_tirosh.txt
@@ -51,4 +51,4 @@ Nek2
 G2e3
 Gas2l3
 Cbx5
-Cenpa
\ No newline at end of file
+Cenpa
diff --git a/scib/resources/g2m_genes_tirosh_hm.txt b/scib/resources/g2m_genes_tirosh_hm.txt
index 4053db0c..d561c8f5 100644
--- a/scib/resources/g2m_genes_tirosh_hm.txt
+++ b/scib/resources/g2m_genes_tirosh_hm.txt
@@ -51,4 +51,4 @@ NEK2
 G2E3
 GAS2L3
 CBX5
-CENPA
\ No newline at end of file
+CENPA
diff --git a/scib/resources/s_genes_tirosh.txt b/scib/resources/s_genes_tirosh.txt
index 17ae40a7..bca6fad4 100644
--- a/scib/resources/s_genes_tirosh.txt
+++ b/scib/resources/s_genes_tirosh.txt
@@ -40,4 +40,4 @@ Clspn
 Pola1
 Chaf1b
 Brip1
-E2f8
\ No newline at end of file
+E2f8
diff --git a/scib/resources/s_genes_tirosh_hm.txt b/scib/resources/s_genes_tirosh_hm.txt
index 35b4aa7e..92ed9288 100644
--- a/scib/resources/s_genes_tirosh_hm.txt
+++ b/scib/resources/s_genes_tirosh_hm.txt
@@ -40,4 +40,4 @@ CLSPN
 POLA1
 CHAF1B
 BRIP1
-E2F8
\ No newline at end of file
+E2F8
diff --git a/scib/trajectory_inference.py b/scib/trajectory_inference.py
index 11ca7cb2..45d05c52 100644
--- a/scib/trajectory_inference.py
+++ b/scib/trajectory_inference.py
@@ -5,9 +5,8 @@
 from . import utils
 
 
-def paga(adata, groups='louvain'):
-    """
-    """
+def paga(adata, groups="louvain"):
+    """ """
     utils.check_adata(adata)
 
     sc.pp.neighbors(adata)
@@ -16,26 +15,32 @@ def paga(adata, groups='louvain'):
 
     fig1, ax1 = plt.subplots()
     sc.pl.umap(adata, size=40, ax=ax1, show=False)
-    sc.pl.paga(adata, pos=adata.uns['paga']['pos'],
-               show=False, node_size_scale=10,
-               node_size_power=1, ax=ax1, text_kwds={'alpha': 0})
+    sc.pl.paga(
+        adata,
+        pos=adata.uns["paga"]["pos"],
+        show=False,
+        node_size_scale=10,
+        node_size_power=1,
+        ax=ax1,
+        text_kwds={"alpha": 0},
+    )
     plt.show()
 
 
-def dpt(adata, group, root, opt='min', comp=0):
+def dpt(adata, group, root, opt="min", comp=0):
     utils.check_adata()
 
     # TODO compute diffmap before
 
     # get root
     stem_mask = np.isin(adata.obs[group], root)
-    if opt == 'min':
-        opt_stem_id = np.argmin(adata.obsm['X_diffmap'][stem_mask, comp])
-    elif opt == 'max':
-        opt_stem_id = np.argmax(adata.obsm['X_diffmap'][stem_mask, comp])
+    if opt == "min":
+        opt_stem_id = np.argmin(adata.obsm["X_diffmap"][stem_mask, comp])
+    elif opt == "max":
+        opt_stem_id = np.argmax(adata.obsm["X_diffmap"][stem_mask, comp])
     else:
         raise ("invalid optimum", opt)
     root_id = np.arange(len(stem_mask))[stem_mask][opt_stem_id]
-    adata.uns['iroot'] = root_id
+    adata.uns["iroot"] = root_id
     # compute pseudotime
     sc.tl.dpt(adata)
diff --git a/scib/utils.py b/scib/utils.py
index 69ccba31..dbff0318 100644
--- a/scib/utils.py
+++ b/scib/utils.py
@@ -4,22 +4,22 @@
 # checker functions for data sanity
 def check_adata(adata):
     if type(adata) is not anndata.AnnData:
-        raise TypeError('Input is not a valid AnnData object')
+        raise TypeError("Input is not a valid AnnData object")
 
 
 def check_batch(batch, obs, verbose=False):
     if batch not in obs:
-        raise ValueError(f'column {batch} is not in obs')
+        raise ValueError(f"column {batch} is not in obs")
     elif verbose:
-        print(f'Object contains {obs[batch].nunique()} batches.')
+        print(f"Object contains {obs[batch].nunique()} batches.")
 
 
 def check_hvg(hvg, adata_var):
     if type(hvg) is not list:
-        raise TypeError('HVG list is not a list')
+        raise TypeError("HVG list is not a list")
     else:
         if not all(i in adata_var.index for i in hvg):
-            raise ValueError('Not all HVGs are in the adata object')
+            raise ValueError("Not all HVGs are in the adata object")
 
 
 def check_sanity(adata, batch, hvg):
@@ -41,7 +41,7 @@ def split_batches(adata, batch, hvg=None, return_categories=False):
     return split
 
 
-def merge_adata(adata_list, sep='-'):
+def merge_adata(adata_list, sep="-"):
     """
     merge adatas from list and remove duplicated obs and var columns
     """
@@ -49,25 +49,32 @@ def merge_adata(adata_list, sep='-'):
     if len(adata_list) == 1:
         return adata_list[0]
 
-    adata = adata_list[0].concatenate(*adata_list[1:], index_unique=None, batch_key='tmp')
-    del adata.obs['tmp']
+    adata = adata_list[0].concatenate(
+        *adata_list[1:], index_unique=None, batch_key="tmp"
+    )
+    del adata.obs["tmp"]
 
     if len(adata.obs.columns) > 0:
         # if there is a column with separator
         if sum(adata.obs.columns.str.contains(sep)) > 0:
-            columns_to_keep = [name.split(sep)[1] == '0' for name in adata.var.columns.values]
+            columns_to_keep = [
+                name.split(sep)[1] == "0" for name in adata.var.columns.values
+            ]
             clean_var = adata.var.loc[:, columns_to_keep]
         else:
             clean_var = adata.var
 
     if len(adata.var.columns) > 0:
         if sum(adata.var.columns.str.contains(sep)) > 0:
-            adata.var = clean_var.rename(columns={name: name.split('-')[0] for name in clean_var.columns.values})
+            adata.var = clean_var.rename(
+                columns={name: name.split("-")[0] for name in clean_var.columns.values}
+            )
 
     return adata
 
 
 def todense(adata):
     import scipy
+
     if isinstance(adata.X, scipy.sparse.csr_matrix):
         adata.X = adata.X.todense()
diff --git a/setup.cfg b/setup.cfg
index d969df02..5232a002 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -17,15 +17,15 @@ author = Malte D. Luecken, Maren Buettner, Daniel C. Strobl, Michaela F. Mueller
 author_email = malte.luecken@helmholtz-muenchen.de, michaela.mueller@helmholtz-muenchen.de
 license = MIT
 url = https://github.com/theislab/scib
-project_urls = 
+project_urls =
 	Pipeline = https://github.com/theislab/scib-pipeline
 	Reproducibility = https://theislab.github.io/scib-reproducibility
 	Bug Tracker = https://github.com/theislab/scib/issues
-keywords = 
+keywords =
 	benchmarking
 	single cell
 	data integration
-classifiers = 
+classifiers =
 	Development Status :: 3 - Alpha
 	Intended Audience :: Developers
 	Intended Audience :: Science/Research
@@ -39,11 +39,11 @@ classifiers =
 build_number = 1
 
 [options]
-packages = 
+packages =
 	scib
 	scib.metrics
 python_requires = >=3.7
-install_requires = 
+install_requires =
 	numpy
 	pandas
 	seaborn
@@ -66,13 +66,13 @@ install_requires =
 zip_safe = False
 
 [options.package_data]
-scib = 
+scib =
 	resources/*.txt
 	knn_graph/*
 
 [options.extras_require]
-test = pytest; pytest-runner; pytest-icdiff
-dev = build; twine; isort; bump2version
+test = pytest; pytest-icdiff; pytest-cov; codecov
+dev = build; twine; isort; bump2version; pre-commit
 docs = sphinx; sphinx_rtd_theme; myst_parser; sphinx-automodapi
 bbknn = bbknn ==1.3.9
 scanorama = scanorama ==1.7.0
@@ -82,3 +82,69 @@ scvi = scvi  ==0.6.7
 trvae = trvae  ==1.1.2
 trvaep = trvaep  ==0.1.0
 desc = desc  ==2.0.3
+
+[tool.isort]
+include_trailing_comma = true
+multi_line_output = 3
+profile = black
+skip_glob = docs/*
+
+[tool.black]
+line-length = 120
+target-version = py38
+include = \.pyi?$
+exclude =
+  .eggs
+  .git
+  .venv
+  build
+  dist
+
+[flake8]
+max-line-length = 88
+ignore =
+    # line break before a binary operator -> black does not adhere to PEP8
+    W503
+    # line break occurred after a binary operator -> black does not adhere to PEP8
+    W504
+    # line too long -> we accept long comment lines; black gets rid of long code lines
+    E501
+    # whitespace before : -> black does not adhere to PEP8
+    E203
+    # missing whitespace after ,', ';', or ':' -> black does not adhere to PEP8
+    E231
+    # continuation line over-indented for hanging indent -> black does not adhere to PEP8
+    E126
+    # E266 too many leading '#' for block comment -> this is fine for indicating sections
+    E262
+    # Do not assign a lambda expression, use a def -> lambda expression assignments are convenient
+    E731
+    # allow I, O, l as variable names -> I is the identity matrix, i, j, k, l is reasonable indexing notation
+    E741
+    # Missing docstring in public package
+    D104
+    # ... imported but unused
+    F401
+    # Missing docstring in public module
+    D100
+    # Missing docstring in __init__
+    D107
+    # Do not perform function calls in argument defaults.
+    B008
+    # line break before binary operator
+    W503
+    # Missing docstring in magic method
+    D105
+    # whitespace before ':'
+    E203
+    # format string does contain non-indexed parameters
+    P101
+    # Unknown interpreted text role "func"
+    RST304
+    # Unnecessary dict call - rewrite as a literal
+    C408
+exclude = .git,__pycache__,build,docs/_build,dist
+per-file-ignores =
+    scib/*: D
+    tests/*: D
+    */__init__.py: F401
diff --git a/setup.py b/setup.py
index ac02561a..45677ce0 100644
--- a/setup.py
+++ b/setup.py
@@ -4,12 +4,18 @@
 from setuptools import setup
 
 try:
-    cmd = 'g++ -std=c++11 -O3 scib/knn_graph/knn_graph.cpp -o scib/knn_graph/knn_graph.o'
-    sys.stdout.write('Compile knn_graph C++ code for LISI metric...\n')
+    cmd = (
+        "g++ -std=c++11 -O3 scib/knn_graph/knn_graph.cpp -o scib/knn_graph/knn_graph.o"
+    )
+    sys.stdout.write("Compile knn_graph C++ code for LISI metric...\n")
     sys.stdout.flush()
-    subprocess.check_output(cmd, stderr=subprocess.STDOUT, shell=True, universal_newlines=True)
+    subprocess.check_output(
+        cmd, stderr=subprocess.STDOUT, shell=True, universal_newlines=True
+    )
 except subprocess.CalledProcessError as exc:
-    sys.stdout.write(f'Failed to compile knn_graph for LISI - skipping...\n{exc.returncode}\n{exc.output}')
+    sys.stdout.write(
+        f"Failed to compile knn_graph for LISI - skipping...\n{exc.returncode}\n{exc.output}"
+    )
     sys.stdout.flush()
 
 setup()


Test patch:
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
new file mode 100644
index 00000000..6c9784a2
--- /dev/null
+++ b/.github/workflows/test.yml
@@ -0,0 +1,66 @@
+name: Test
+
+on:
+    push:
+        branches: [main]
+    pull_request:
+        branches: [main]
+
+jobs:
+    lint:
+        runs-on: ubuntu-latest
+        steps:
+        -   uses: actions/checkout@v2
+        -   name: Set up Python
+            uses: actions/setup-python@v2
+            with:
+                python-version: 3.9
+
+        -   uses: actions/cache@v2
+            with:
+                path: ~/.cache/pre-commit
+                key: precommit-${{ env.pythonLocation }}-${{ hashFiles('**/.pre-commit-config.yaml') }}
+        -   name: Install dependencies
+            run: |
+                python -m pip install --upgrade pre-commit
+
+        -   name: Run pre-commit
+            run: |
+                pre-commit run --color=always --show-diff-on-failure --all-files
+
+    test:
+        runs-on: ${{ matrix.os }}
+        strategy:
+            matrix:
+                python: [3.7, 3.9]
+                os: [ubuntu-latest, macos-latest]
+        steps:
+        -   uses: actions/checkout@v2
+        -   name: Set up Python ${{ matrix.python }}
+            uses: actions/setup-python@v2
+            with:
+                python-version: ${{ matrix.python }}
+
+        -   name: Install build dependencies
+            run: |
+                python -m pip install --upgrade pip
+                pip install build twine
+
+        -   name: Build a binary wheel and a source tarball
+            run: |
+                python -m build --sdist --wheel --outdir dist/
+                twine check dist/*
+
+        -   name: Install package with test dependencies
+            run: |
+                pip install dist/*.whl
+                pip install '.[test]'
+
+        -   name: Test with pytest
+            run: |
+                python -m pytest --cov --cov-append --cov-report=term-missing -vv
+        -   name: Upload coverage
+            env:
+                CODECOV_NAME: ${{ matrix.os }}-${{ matrix.python }}
+            run: |
+                codecov --no-color --required --flags unittest
diff --git a/tests/common.py b/tests/common.py
index 8b71f35b..f45834ca 100644
--- a/tests/common.py
+++ b/tests/common.py
@@ -1,14 +1,11 @@
-import pytest
-import scib
-import scanpy as sc
-import numpy as np
-import pandas as pd
+import logging
 import os
+import subprocess
 import warnings
 
-warnings.filterwarnings('ignore')
-import subprocess
-import logging
+import scanpy as sc
+
+warnings.filterwarnings("ignore")
 
 CORES = "1"
 LOGGER = logging.getLogger(__name__)
@@ -23,14 +20,24 @@ def create_if_missing(dir):
         os.mkdir(dir)
 
 
-def run(cmd, dir_path, report_stdout=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, **kwargs):
+def run(
+    cmd,
+    dir_path,
+    report_stdout=True,
+    stdout=subprocess.PIPE,
+    stderr=subprocess.PIPE,
+    **kwargs,
+):
     LOGGER.info(cmd)
     shell = isinstance(cmd, str)
     response = subprocess.run(
-        cmd, cwd=dir_path, shell=shell,
-        stdout=stdout, stderr=stderr,
+        cmd,
+        cwd=dir_path,
+        shell=shell,
+        stdout=stdout,
+        stderr=stderr,
         universal_newlines=True,
-        **kwargs
+        **kwargs,
     )
     try:
         response.check_returncode()
@@ -47,14 +54,14 @@ def runR(r_cmd, dir_path, report_stdout=False):
 
 
 def add_embed(adata, type_):
-    if type_ == 'pca':
-        if 'X_pca' in adata.obsm:
-            mtx = adata.obsm['X_pca']
+    if type_ == "pca":
+        if "X_pca" in adata.obsm:
+            mtx = adata.obsm["X_pca"]
         else:
-            mtx = sc.tl.pca(adata, copy=True).obsm['X_pca']
-    elif type_ == 'full':
+            mtx = sc.tl.pca(adata, copy=True).obsm["X_pca"]
+    elif type_ == "full":
         mtx = adata.X
     else:
         raise ValueError(f"'{type_}' not a valid embedding type")
-    adata.obsm['X_emb'] = mtx
+    adata.obsm["X_emb"] = mtx
     return adata
diff --git a/tests/conftest.py b/tests/conftest.py
index 13e11b4b..fc75b606 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,21 +1,21 @@
-from .common import *
+import numpy as np
+import pytest
+import scanpy as sc
+
+import scib
 
 
 @pytest.fixture(scope="session")
 def adata_paul15_template():
     adata = sc.datasets.paul15()
-    adata.obs['celltype'] = adata.obs['paul15_clusters']
+    adata.obs["celltype"] = adata.obs["paul15_clusters"]
     np.random.seed(42)
-    adata.obs['batch'] = np.random.randint(1, 5, adata.n_obs)
-    adata.obs['batch'] = adata.obs['batch'].astype(str)
-    adata.obs['batch'] = adata.obs['batch'].astype("category")
-    adata.layers['counts'] = adata.X
+    adata.obs["batch"] = np.random.randint(1, 5, adata.n_obs)
+    adata.obs["batch"] = adata.obs["batch"].astype(str)
+    adata.obs["batch"] = adata.obs["batch"].astype("category")
+    adata.layers["counts"] = adata.X
     scib.preprocessing.reduce_data(
-        adata,
-        pca=False,
-        n_top_genes=None,
-        neighbors=False,
-        umap=False
+        adata, pca=False, n_top_genes=None, neighbors=False, umap=False
     )
     yield adata
     del adata
@@ -27,7 +27,7 @@ def adata_pbmc_template():
     # quick fix for broken dataset paths, should be removed with scanpy>=1.6.0
     adata_ref = sc.read(
         "pbmc3k_processed.h5ad",
-        backup_url="https://raw.githubusercontent.com/chanzuckerberg/cellxgene/main/example-dataset/pbmc3k.h5ad"
+        backup_url="https://raw.githubusercontent.com/chanzuckerberg/cellxgene/main/example-dataset/pbmc3k.h5ad",
     )
     adata = sc.datasets.pbmc68k_reduced()
 
@@ -40,14 +40,16 @@ def adata_pbmc_template():
     sc.tl.umap(adata_ref)
 
     # merge cell type labels
-    sc.tl.ingest(adata, adata_ref, obs='louvain')
-    adata_concat = adata_ref.concatenate(adata, batch_categories=['ref', 'new'])
-    adata_concat.obs.louvain = adata_concat.obs.louvain.astype('category')
+    sc.tl.ingest(adata, adata_ref, obs="louvain")
+    adata_concat = adata_ref.concatenate(adata, batch_categories=["ref", "new"])
+    adata_concat.obs.louvain = adata_concat.obs.louvain.astype("category")
     # fix category ordering
-    adata_concat.obs.louvain.cat.reorder_categories(adata_ref.obs.louvain.cat.categories, inplace=True)
-    adata_concat.obs['celltype'] = adata_concat.obs['louvain']
+    adata_concat.obs.louvain.cat.reorder_categories(
+        adata_ref.obs.louvain.cat.categories, inplace=True
+    )
+    adata_concat.obs["celltype"] = adata_concat.obs["louvain"]
 
-    del adata_concat.obs['louvain']
+    del adata_concat.obs["louvain"]
     del adata_concat.uns
     del adata_concat.obsm
     del adata_concat.varm
@@ -74,11 +76,7 @@ def adata(adata_pbmc_template):
 def adata_pca(adata):
     adata_obj = adata
     scib.pp.reduce_data(
-        adata_obj,
-        pca=True,
-        n_top_genes=200,
-        neighbors=False,
-        umap=False
+        adata_obj, pca=True, n_top_genes=200, neighbors=False, umap=False
     )
     yield adata_obj
 
@@ -87,11 +85,7 @@ def adata_pca(adata):
 def adata_neighbors(adata):
     adata_obj = adata
     scib.pp.reduce_data(
-        adata_obj,
-        pca=True,
-        n_top_genes=200,
-        neighbors=True,
-        umap=False
+        adata_obj, pca=True, n_top_genes=200, neighbors=True, umap=False
     )
     yield adata_obj
 
@@ -100,9 +94,6 @@ def adata_neighbors(adata):
 def adata_clustered(adata_neighbors):
     adata_obj = adata_neighbors
     scib.cl.opt_louvain(
-        adata_obj,
-        cluster_key='cluster',
-        label_key='celltype',
-        verbose=True
+        adata_obj, cluster_key="cluster", label_key="celltype", verbose=True
     )
     yield adata_obj
diff --git a/tests/metrics/test_all.py b/tests/metrics/test_all.py
index 6ee84d29..4f5b7b63 100644
--- a/tests/metrics/test_all.py
+++ b/tests/metrics/test_all.py
@@ -1,13 +1,16 @@
-from tests.common import *
+import scanpy as sc
+
+import scib
+from tests.common import LOGGER
 
 
 def test_fast(adata_neighbors):
     metrics_df = scib.me.metrics_fast(
         adata_neighbors,
         adata_neighbors,
-        batch_key='batch',
-        label_key='celltype',
-        embed='X_pca',
+        batch_key="batch",
+        label_key="celltype",
+        embed="X_pca",
     )
 
     for score in metrics_df:
@@ -23,9 +26,9 @@ def test_slim(adata_paul15):
     metrics_df = scib.me.metrics_slim(
         adata_paul15,
         adata_paul15,
-        batch_key='batch',
-        label_key='celltype',
-        embed='X_pca',
+        batch_key="batch",
+        label_key="celltype",
+        embed="X_pca",
     )
 
     for score in metrics_df:
diff --git a/tests/metrics/test_beyond_label_metrics.py b/tests/metrics/test_beyond_label_metrics.py
index 65b56156..77818dcf 100644
--- a/tests/metrics/test_beyond_label_metrics.py
+++ b/tests/metrics/test_beyond_label_metrics.py
@@ -1,4 +1,7 @@
-from tests.common import *
+import pandas as pd
+
+import scib
+from tests.common import LOGGER, assert_near_exact
 
 
 def test_cell_cycle(adata_paul15):
@@ -9,10 +12,10 @@ def test_cell_cycle(adata_paul15):
     score = scib.me.cell_cycle(
         adata,
         adata_int,
-        batch_key='batch',
-        organism='mouse',
+        batch_key="batch",
+        organism="mouse",
         # recompute_cc=True,
-        verbose=True
+        verbose=True,
     )
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 1, diff=1e-12)
@@ -26,25 +29,20 @@ def test_cell_cycle_all(adata_paul15):
     scores_df = scib.me.cell_cycle(
         adata,
         adata_int,
-        batch_key='batch',
-        organism='mouse',
+        batch_key="batch",
+        organism="mouse",
         # recompute_cc=True,
         agg_func=None,
-        verbose=True
+        verbose=True,
     )
     LOGGER.info(f"\nscore: {scores_df}")
     assert isinstance(scores_df, pd.DataFrame)
-    for i in scores_df['score']:
+    for i in scores_df["score"]:
         assert_near_exact(i, 1, diff=1e-12)
 
 
 def test_hvg_overlap(adata):
     adata_int = adata.copy()
-    score = scib.me.hvg_overlap(
-        adata_int,
-        adata,
-        batch='batch',
-        n_hvg=500
-    )
+    score = scib.me.hvg_overlap(adata_int, adata, batch="batch", n_hvg=500)
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 1, diff=1e-12)
diff --git a/tests/metrics/test_clisi.py b/tests/metrics/test_clisi.py
index 9fa42a8c..ba5da653 100644
--- a/tests/metrics/test_clisi.py
+++ b/tests/metrics/test_clisi.py
@@ -1,14 +1,15 @@
-from tests.common import *
+import scib
+from tests.common import LOGGER, assert_near_exact
 
 
 def test_clisi_full(adata):
     score = scib.me.clisi_graph(
         adata,
-        batch_key='batch',
-        label_key='celltype',
+        batch_key="batch",
+        label_key="celltype",
         scale=True,
-        type_='full',
-        verbose=True
+        type_="full",
+        verbose=True,
     )
 
     LOGGER.info(f"score: {score}")
@@ -16,14 +17,14 @@ def test_clisi_full(adata):
 
 
 def test_clisi_embed(adata_neighbors):
-    adata_neighbors.obsm['X_emb'] = adata_neighbors.obsm['X_pca']
+    adata_neighbors.obsm["X_emb"] = adata_neighbors.obsm["X_pca"]
     score = scib.me.clisi_graph(
         adata_neighbors,
-        batch_key='batch',
-        label_key='celltype',
+        batch_key="batch",
+        label_key="celltype",
         scale=True,
-        type_='embed',
-        verbose=True
+        type_="embed",
+        verbose=True,
     )
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.982, diff=1e-2)
@@ -32,11 +33,11 @@ def test_clisi_embed(adata_neighbors):
 def test_clisi_knn(adata_neighbors):
     score = scib.me.clisi_graph(
         adata_neighbors,
-        batch_key='batch',
-        label_key='celltype',
+        batch_key="batch",
+        label_key="celltype",
         scale=True,
-        type_='graph',
-        verbose=True
+        type_="graph",
+        verbose=True,
     )
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.982, diff=1e-2)
@@ -46,12 +47,12 @@ def test_clisi_full_parallel(adata):
     # test parallel
     score = scib.me.clisi_graph(
         adata,
-        batch_key='batch',
-        label_key='celltype',
+        batch_key="batch",
+        label_key="celltype",
         scale=True,
-        type_='full',
+        type_="full",
         n_cores=2,
-        verbose=True
+        verbose=True,
     )
 
     LOGGER.info(f"score: {score}")
@@ -59,15 +60,15 @@ def test_clisi_full_parallel(adata):
 
 
 def test_clisi_embed_parallel(adata_neighbors):
-    adata_neighbors.obsm['X_emb'] = adata_neighbors.obsm['X_pca']
+    adata_neighbors.obsm["X_emb"] = adata_neighbors.obsm["X_pca"]
     score = scib.me.clisi_graph(
         adata_neighbors,
-        batch_key='batch',
-        label_key='celltype',
+        batch_key="batch",
+        label_key="celltype",
         scale=True,
-        type_='embed',
+        type_="embed",
         n_cores=2,
-        verbose=True
+        verbose=True,
     )
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.982, diff=1e-2)
@@ -76,12 +77,12 @@ def test_clisi_embed_parallel(adata_neighbors):
 def test_clisi_knn_parallel(adata_neighbors):
     score = scib.me.clisi_graph(
         adata_neighbors,
-        batch_key='batch',
-        label_key='celltype',
+        batch_key="batch",
+        label_key="celltype",
         scale=True,
-        type_='graph',
+        type_="graph",
         n_cores=2,
-        verbose=True
+        verbose=True,
     )
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.982, diff=1e-2)
diff --git a/tests/metrics/test_cluster_metrics.py b/tests/metrics/test_cluster_metrics.py
index 819e9b0c..e868f07c 100644
--- a/tests/metrics/test_cluster_metrics.py
+++ b/tests/metrics/test_cluster_metrics.py
@@ -1,34 +1,35 @@
-from tests.common import *
+import scib
+from tests.common import LOGGER, assert_near_exact
 
 
 def test_nmi_trivial(adata):
-    score = scib.me.nmi(adata, 'celltype', 'celltype')
+    score = scib.me.nmi(adata, "celltype", "celltype")
     assert_near_exact(score, 1, diff=1e-12)
 
 
 def test_ari_trivial(adata):
-    score = scib.me.ari(adata, 'celltype', 'celltype')
+    score = scib.me.ari(adata, "celltype", "celltype")
     assert_near_exact(score, 1, diff=1e-12)
 
 
 def test_nmi(adata_neighbors):
     _, _, nmi_all = scib.cl.opt_louvain(
         adata_neighbors,
-        label_key='celltype',
-        cluster_key='cluster',
+        label_key="celltype",
+        cluster_key="cluster",
         function=scib.me.nmi,
         plot=False,
         inplace=True,
         force=True,
-        verbose=True
+        verbose=True,
     )
 
-    for score in nmi_all['score']:
+    for score in nmi_all["score"]:
         assert 0 <= score <= 1
 
 
 def test_ari(adata_clustered):
-    score = scib.me.ari(adata_clustered, group1='cluster', group2='celltype')
+    score = scib.me.ari(adata_clustered, group1="cluster", group2="celltype")
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.7614422905830917, diff=1e-2)
 
@@ -36,11 +37,11 @@ def test_ari(adata_clustered):
 def test_isolated_labels_F1(adata_neighbors):
     score = scib.me.isolated_labels(
         adata_neighbors,
-        label_key='celltype',
-        batch_key='batch',
-        embed='X_pca',
+        label_key="celltype",
+        batch_key="batch",
+        embed="X_pca",
         cluster=True,
-        verbose=True
+        verbose=True,
     )
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.5581395348837209, diff=1e-12)
diff --git a/tests/metrics/test_graph_connectivity.py b/tests/metrics/test_graph_connectivity.py
index 13144d25..23dc8989 100644
--- a/tests/metrics/test_graph_connectivity.py
+++ b/tests/metrics/test_graph_connectivity.py
@@ -1,7 +1,8 @@
-from tests.common import *
+import scib
+from tests.common import LOGGER, assert_near_exact
 
 
 def test_graph_connectivity(adata_neighbors):
-    score = scib.me.graph_connectivity(adata_neighbors, label_key='celltype')
+    score = scib.me.graph_connectivity(adata_neighbors, label_key="celltype")
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.9670013350457753, diff=1e-3)
diff --git a/tests/metrics/test_ilisi.py b/tests/metrics/test_ilisi.py
index 9721a04b..d59ca8e8 100644
--- a/tests/metrics/test_ilisi.py
+++ b/tests/metrics/test_ilisi.py
@@ -1,13 +1,10 @@
-from tests.common import *
+import scib
+from tests.common import LOGGER, assert_near_exact
 
 
 def test_ilisi_full(adata):
     score = scib.me.ilisi_graph(
-        adata,
-        batch_key='batch',
-        scale=True,
-        type_='full',
-        verbose=True
+        adata, batch_key="batch", scale=True, type_="full", verbose=True
     )
 
     LOGGER.info(f"score: {score}")
@@ -15,13 +12,9 @@ def test_ilisi_full(adata):
 
 
 def test_ilisi_embed(adata_neighbors):
-    adata_neighbors.obsm['X_emb'] = adata_neighbors.obsm['X_pca']
+    adata_neighbors.obsm["X_emb"] = adata_neighbors.obsm["X_pca"]
     score = scib.me.ilisi_graph(
-        adata_neighbors,
-        batch_key='batch',
-        scale=True,
-        type_='embed',
-        verbose=True
+        adata_neighbors, batch_key="batch", scale=True, type_="embed", verbose=True
     )
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.238, diff=1e-2)
@@ -29,11 +22,7 @@ def test_ilisi_embed(adata_neighbors):
 
 def test_ilisi_knn(adata_neighbors):
     score = scib.me.ilisi_graph(
-        adata_neighbors,
-        batch_key='batch',
-        scale=True,
-        type_='graph',
-        verbose=True
+        adata_neighbors, batch_key="batch", scale=True, type_="graph", verbose=True
     )
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.238, diff=1e-2)
@@ -41,11 +30,7 @@ def test_ilisi_knn(adata_neighbors):
 
 def test_ilisi_full_parallel(adata):
     score = scib.me.ilisi_graph(
-        adata,
-        batch_key='batch',
-        scale=True,
-        type_='full',
-        verbose=True
+        adata, batch_key="batch", scale=True, type_="full", verbose=True
     )
 
     LOGGER.info(f"score: {score}")
@@ -53,13 +38,9 @@ def test_ilisi_full_parallel(adata):
 
 
 def test_ilisi_embed_parallel(adata_neighbors):
-    adata_neighbors.obsm['X_emb'] = adata_neighbors.obsm['X_pca']
+    adata_neighbors.obsm["X_emb"] = adata_neighbors.obsm["X_pca"]
     score = scib.me.ilisi_graph(
-        adata_neighbors,
-        batch_key='batch',
-        scale=True,
-        type_='embed',
-        verbose=True
+        adata_neighbors, batch_key="batch", scale=True, type_="embed", verbose=True
     )
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.238, diff=1e-2)
@@ -67,11 +48,7 @@ def test_ilisi_embed_parallel(adata_neighbors):
 
 def test_ilisi_knn_parallel(adata_neighbors):
     score = scib.me.ilisi_graph(
-        adata_neighbors,
-        batch_key='batch',
-        scale=True,
-        type_='graph',
-        verbose=True
+        adata_neighbors, batch_key="batch", scale=True, type_="graph", verbose=True
     )
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.238, diff=1e-2)
diff --git a/tests/metrics/test_kbet.py b/tests/metrics/test_kbet.py
index 1371ae77..5f5520a5 100644
--- a/tests/metrics/test_kbet.py
+++ b/tests/metrics/test_kbet.py
@@ -1,14 +1,12 @@
 import numpy as np
 
-from tests.common import *
+import scib
+from tests.common import LOGGER
 
 
 def test_kbet(adata_pca):
     score = scib.me.kBET(
-        adata_pca,
-        batch_key='batch',
-        label_key='celltype',
-        embed='X_pca'
+        adata_pca, batch_key="batch", label_key="celltype", embed="X_pca"
     )
     LOGGER.info(f"score: {score}")
     assert np.isnan(score)
diff --git a/tests/metrics/test_pcr_metrics.py b/tests/metrics/test_pcr_metrics.py
index 85a5cea0..167c1798 100644
--- a/tests/metrics/test_pcr_metrics.py
+++ b/tests/metrics/test_pcr_metrics.py
@@ -1,4 +1,5 @@
-from tests.common import *
+import scib
+from tests.common import LOGGER, add_embed, assert_near_exact
 
 
 def test_pc_regression(adata):
@@ -8,17 +9,14 @@ def test_pc_regression(adata):
 def test_pcr_batch(adata):
     # no PCA precomputed
     score = scib.me.pcr_comparison(
-        adata, adata,
-        covariate='batch',
-        n_comps=50,
-        scale=True
+        adata, adata, covariate="batch", n_comps=50, scale=True
     )
     LOGGER.info(f"no PCA precomputed: {score}")
     assert_near_exact(score, 0, diff=1e-6)
 
 
 def test_pcr_batch_precomputed(adata_pca):
-    score = scib.me.pcr_comparison(adata_pca, adata_pca, covariate='batch', scale=True)
+    score = scib.me.pcr_comparison(adata_pca, adata_pca, covariate="batch", scale=True)
     LOGGER.info(f"precomputed PCA: {score}")
     assert_near_exact(score, 0, diff=1e-6)
 
@@ -27,11 +25,11 @@ def test_pcr_batch_embedding(adata):
     # use different embedding
     score = scib.me.pcr_comparison(
         adata_pre=adata,
-        adata_post=add_embed(adata, type_='full'),
-        covariate='batch',
-        embed='X_emb',
+        adata_post=add_embed(adata, type_="full"),
+        covariate="batch",
+        embed="X_emb",
         n_comps=50,
-        scale=True
+        scale=True,
     )
     LOGGER.info(f"using embedding: {score}")
     assert_near_exact(score, 0, diff=1e-6)
diff --git a/tests/metrics/test_silhouette_metrics.py b/tests/metrics/test_silhouette_metrics.py
index 33a0df3a..abae67a3 100644
--- a/tests/metrics/test_silhouette_metrics.py
+++ b/tests/metrics/test_silhouette_metrics.py
@@ -1,12 +1,10 @@
-from tests.common import *
+import scib
+from tests.common import LOGGER, assert_near_exact
 
 
 def test_silhouette(adata_pca):
     score = scib.me.silhouette(
-        adata_pca,
-        group_key='celltype',
-        embed='X_pca',
-        scale=True
+        adata_pca, group_key="celltype", embed="X_pca", scale=True
     )
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.5626532882452011, diff=1e-3)
@@ -15,11 +13,11 @@ def test_silhouette(adata_pca):
 def test_silhouette_batch(adata_pca):
     score = scib.me.silhouette_batch(
         adata_pca,
-        batch_key='batch',
-        group_key='celltype',
-        embed='X_pca',
+        batch_key="batch",
+        group_key="celltype",
+        embed="X_pca",
         scale=True,
-        verbose=False
+        verbose=False,
     )
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.9014384369842835, diff=1e-3)
@@ -28,11 +26,11 @@ def test_silhouette_batch(adata_pca):
 def test_isolated_labels_silhouette(adata_pca):
     score = scib.me.isolated_labels(
         adata_pca,
-        label_key='celltype',
-        batch_key='batch',
-        embed='X_pca',
+        label_key="celltype",
+        batch_key="batch",
+        embed="X_pca",
         cluster=False,
-        verbose=True
+        verbose=True,
     )
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.6101431176066399, diff=1e-3)
diff --git a/tests/metrics/test_trajectory.py b/tests/metrics/test_trajectory.py
index 6a947d64..40f50cb0 100644
--- a/tests/metrics/test_trajectory.py
+++ b/tests/metrics/test_trajectory.py
@@ -1,18 +1,21 @@
-from tests.common import *
+import scanpy as sc
+
+import scib
+from tests.common import LOGGER, assert_near_exact
 
 
 def test_trajectory(adata_neighbors):
     adata_int = adata_neighbors.copy()
 
-    adata_neighbors.uns['iroot'] = 0
+    adata_neighbors.uns["iroot"] = 0
     sc.tl.diffmap(adata_neighbors)
     sc.tl.dpt(adata_neighbors)
 
     score = scib.me.trajectory_conservation(
         adata_pre=adata_neighbors,
         adata_post=adata_int,
-        label_key='celltype',
-        pseudotime_key='dpt_pseudotime'
+        label_key="celltype",
+        pseudotime_key="dpt_pseudotime",
     )
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.9561, diff=1e-5)
@@ -21,16 +24,16 @@ def test_trajectory(adata_neighbors):
 def test_trajectory_batch(adata_neighbors):
     adata_int = adata_neighbors.copy()
 
-    adata_neighbors.uns['iroot'] = 0
+    adata_neighbors.uns["iroot"] = 0
     sc.tl.diffmap(adata_neighbors)
     sc.tl.dpt(adata_neighbors)
 
     score = scib.me.trajectory_conservation(
         adata_pre=adata_neighbors,
         adata_post=adata_int,
-        label_key='celltype',
-        batch_key='batch',
-        pseudotime_key='dpt_pseudotime'
+        label_key="celltype",
+        batch_key="batch",
+        pseudotime_key="dpt_pseudotime",
     )
     LOGGER.info(f"score: {score}")
     assert_near_exact(score, 0.96317, diff=1e-5)
diff --git a/tests/preprocessing/test_clustering.py b/tests/preprocessing/test_clustering.py
index 5e841751..4c65904b 100644
--- a/tests/preprocessing/test_clustering.py
+++ b/tests/preprocessing/test_clustering.py
@@ -1,13 +1,15 @@
-from tests.common import *
+import pandas as pd
+
+import scib
 
 
 def test_cluster(adata_neighbors):
     _, _, score_all, clustering = scib.cl.opt_louvain(
         adata_neighbors,
-        label_key='celltype',
-        cluster_key='cluster',
+        label_key="celltype",
+        cluster_key="cluster",
         plot=False,
-        inplace=False
+        inplace=False,
     )
     assert isinstance(score_all, pd.DataFrame)
     assert isinstance(clustering, pd.Series)
diff --git a/tests/preprocessing/test_preprocessing.py b/tests/preprocessing/test_preprocessing.py
index 0184a781..e29b4287 100644
--- a/tests/preprocessing/test_preprocessing.py
+++ b/tests/preprocessing/test_preprocessing.py
@@ -1,11 +1,12 @@
-import scanpy as sc
 import numpy as np
+import scanpy as sc
+
 import scib
 
 
 def test_scale():
     adata = sc.datasets.blobs()
-    scib.pp.scale_batch(adata, 'blobs')
-    split = scib.utils.split_batches(adata, 'blobs')
+    scib.pp.scale_batch(adata, "blobs")
+    split = scib.utils.split_batches(adata, "blobs")
     for i in split:
-        assert np.allclose(i.X.mean(0), np.zeros((0,adata.n_vars)))
+        assert np.allclose(i.X.mean(0), np.zeros((0, adata.n_vars)))