diff --git a/.github/workflows/python-plain-run-test.yml b/.github/workflows/python-plain-run-test.yml
index c48b72e97..d11c20875 100644
--- a/.github/workflows/python-plain-run-test.yml
+++ b/.github/workflows/python-plain-run-test.yml
@@ -9,6 +9,21 @@ on:
     branches: [ master ]
 
 jobs:
+  check-code:
+    name: Check code quality
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout branch
+        uses: actions/checkout@v2
+      - name: Setup python
+        uses: actions/setup-python@v2
+        with:
+          python-version: '3.9'
+          architecture: x64
+      - name: Install dependencies
+        run: pip install isort black
+      - name: Check code quality
+        run: isort --profile black --check dynamo tests && black --check dynamo tests && echo OK
   build:
     runs-on: ubuntu-latest
     strategy:
@@ -26,20 +41,11 @@ jobs:
       run: |
         python -m pip install --upgrade pip
         pip install git+https://github.com/pypa/pip.git@611e9253ff74ca5b8ce438c2423373ab0d44c961
-        pip install flake8 pytest
+        pip install pytest
         if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
     - name: Install package
       run: |
         pip install -e .[test]
-    - 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: Test with pytest
-#      run: |
-#        pytest
     - name: Run tests and collect coverage
       run: |
         pip install pytest-cov
diff --git a/dynamo/external/hodge.py b/dynamo/external/hodge.py
index 8864434bf..ba0d2a827 100644
--- a/dynamo/external/hodge.py
+++ b/dynamo/external/hodge.py
@@ -10,7 +10,6 @@
 from scipy.sparse import csr_matrix, issparse
 
 from ..dynamo_logger import main_finish_progress, main_info, main_log_time
-
 # from ..vectorfield.scVectorField import graphize_vecfld
 from ..tools.graph_calculus import divergence, graphize_velocity, potential
 from ..tools.sampling import sample_by_velocity, trn
@@ -21,7 +20,9 @@
     div,
     potential,
 )"""
-from ..tools.connectivity import check_and_recompute_neighbors, generate_neighbor_keys
+
+from ..tools.connectivity import (generate_neighbor_keys,
+                                  check_and_recompute_neighbors)
 
 
 def ddhodge(
@@ -218,7 +219,11 @@ def func(x):
             W.dot(ddhodge_div),
             W.dot(potential_),
         )
-        (adata.obs[prefix + "ddhodge_sampled"], adata.obs[prefix + "ddhodge_div"], adata.obs[prefix + "potential"],) = (
+        (
+            adata.obs[prefix + "ddhodge_sampled"],
+            adata.obs[prefix + "ddhodge_div"],
+            adata.obs[prefix + "potential"],
+        ) = (
             False,
             0,
             0,
diff --git a/dynamo/prediction/fate.py b/dynamo/prediction/fate.py
index bf82b8384..14325a5af 100755
--- a/dynamo/prediction/fate.py
+++ b/dynamo/prediction/fate.py
@@ -511,7 +511,6 @@ def fate_bias(
         while True:
             is_dist_larger_than_threshold = distances.flatten() < dist_threshold * median_dist
             if any(is_dist_larger_than_threshold):
-
                 # let us diffuse one step further to identify cells from terminal cell types in case
                 # cells with indices are all close to some random progenitor cells.
                 if hasattr(nbrs, "query"):
diff --git a/dynamo/preprocessing/Preprocessor.py b/dynamo/preprocessing/Preprocessor.py
index 87abfffe5..fc84f9020 100644
--- a/dynamo/preprocessing/Preprocessor.py
+++ b/dynamo/preprocessing/Preprocessor.py
@@ -9,21 +9,13 @@
 from anndata import AnnData
 
 from ..configuration import DKM
-from ..dynamo_logger import (
-    LoggerManager,
-    main_debug,
-    main_info,
-    main_info_insert_adata,
-    main_warning,
-)
+from ..dynamo_logger import (LoggerManager, main_debug, main_info,
+                             main_info_insert_adata, main_warning)
 from ..tools.connectivity import neighbors as default_neighbors
 from ..tools.utils import update_dict
 from .cell_cycle import cell_cycle_scores
-from .external import (
-    normalize_layers_pearson_residuals,
-    sctransform,
-    select_genes_by_pearson_residuals,
-)
+from .external import (normalize_layers_pearson_residuals, sctransform,
+                       select_genes_by_pearson_residuals)
 from .gene_selection import select_genes_by_seurat_recipe, select_genes_monocle
 from .normalization import calc_sz_factor, normalize
 from .pca import pca
@@ -32,15 +24,9 @@
 from .QC import filter_genes_by_outliers as monocle_filter_genes_by_outliers
 from .QC import regress_out_parallel
 from .transform import Freeman_Tukey, log, log1p, log2
-from .utils import (
-    _infer_labeling_experiment_type,
-    calc_new_to_total_ratio,
-    collapse_species_adata,
-    convert2symbol,
-    convert_layers2csr,
-    detect_experiment_datatype,
-    unique_var_obs_adata,
-)
+from .utils import (_infer_labeling_experiment_type, calc_new_to_total_ratio,
+                    collapse_species_adata, convert2symbol, convert_layers2csr,
+                    detect_experiment_datatype, unique_var_obs_adata)
 
 
 class Preprocessor:
diff --git a/dynamo/preprocessing/QC.py b/dynamo/preprocessing/QC.py
index 4d370c169..ea0cb1d63 100644
--- a/dynamo/preprocessing/QC.py
+++ b/dynamo/preprocessing/QC.py
@@ -13,16 +13,10 @@
 from scipy.sparse import issparse, spmatrix
 
 from ..configuration import DKM
-from ..dynamo_logger import (
-    LoggerManager,
-    main_debug,
-    main_exception,
-    main_finish_progress,
-    main_info,
-    main_info_insert_adata_obs,
-    main_log_time,
-    main_warning,
-)
+from ..dynamo_logger import (LoggerManager, main_debug, main_exception,
+                             main_finish_progress, main_info,
+                             main_info_insert_adata_obs, main_log_time,
+                             main_warning)
 from .utils import get_inrange_shared_counts_mask
 
 
diff --git a/dynamo/preprocessing/__init__.py b/dynamo/preprocessing/__init__.py
index a309d76da..e5663206f 100755
--- a/dynamo/preprocessing/__init__.py
+++ b/dynamo/preprocessing/__init__.py
@@ -2,22 +2,13 @@
 """
 
 from .cell_cycle import cell_cycle_scores
-from .deprecated import (
-    calc_sz_factor_legacy,
-    cook_dist,
-    filter_cells_legacy,
-    normalize_cell_expr_by_size_factors,
-    recipe_monocle,
-    recipe_velocyto,
-)
+from .deprecated import (calc_sz_factor_legacy, cook_dist, filter_cells_legacy,
+                         normalize_cell_expr_by_size_factors, recipe_monocle,
+                         recipe_velocyto)
 from .dynast import lambda_correction
-from .external import (
-    harmony_debatch,
-    integrate,
-    normalize_layers_pearson_residuals,
-    sctransform,
-    select_genes_by_pearson_residuals,
-)
+from .external import (harmony_debatch, integrate,
+                       normalize_layers_pearson_residuals, sctransform,
+                       select_genes_by_pearson_residuals)
 from .normalization import calc_sz_factor, normalize
 from .pca import pca, top_pca_genes
 from .QC import (
@@ -45,12 +36,8 @@
 normalize_cells = normalize
 
 from .CnmfPreprocessor import CnmfPreprocessor
-from .gene_selection import (
-    calc_dispersion_by_svr,
-    calc_Gini,
-    highest_frac_genes,
-    select_genes_monocle,
-)
+from .gene_selection import (calc_dispersion_by_svr, calc_Gini,
+                             highest_frac_genes, select_genes_monocle)
 from .Preprocessor import Preprocessor
 
 __all__ = [
diff --git a/dynamo/preprocessing/deprecated.py b/dynamo/preprocessing/deprecated.py
index 522a5777a..025231e35 100644
--- a/dynamo/preprocessing/deprecated.py
+++ b/dynamo/preprocessing/deprecated.py
@@ -19,35 +19,24 @@
 from sklearn.decomposition import FastICA
 
 from ..configuration import DKM, DynamoAdataConfig, DynamoAdataKeyManager
-from ..dynamo_logger import (
-    LoggerManager,
-    main_debug,
-    main_info,
-    main_info_insert_adata_obsm,
-    main_warning,
-)
+from ..dynamo_logger import (LoggerManager, main_debug, main_info,
+                             main_info_insert_adata_obsm, main_warning)
+from ..tools.utils import update_dict
 from ..utils import copy_adata
 from .cell_cycle import cell_cycle_scores
 from .gene_selection import calc_dispersion_by_svr
-from .normalization import calc_sz_factor, get_sz_exprs, normalize_mat_monocle, sz_util
+from .normalization import (calc_sz_factor, get_sz_exprs,
+                            normalize_mat_monocle, sz_util)
 from .pca import pca
 from .QC import basic_stats, filter_genes_by_clusters, filter_genes_by_outliers
 from .transform import _Freeman_Tukey
-from .utils import (
-    _infer_labeling_experiment_type,
-    add_noise_to_duplicates,
-    calc_new_to_total_ratio,
-    collapse_species_adata,
-    compute_gene_exp_fraction,
-    convert2symbol,
-    convert_layers2csr,
-    detect_experiment_datatype,
-    get_inrange_shared_counts_mask,
-    get_nan_or_inf_data_bool_mask,
-    get_svr_filter,
-    merge_adata_attrs,
-    unique_var_obs_adata,
-)
+from .utils import (_infer_labeling_experiment_type, add_noise_to_duplicates,
+                    calc_new_to_total_ratio, collapse_species_adata,
+                    compute_gene_exp_fraction, convert2symbol,
+                    convert_layers2csr, detect_experiment_datatype,
+                    get_inrange_shared_counts_mask,
+                    get_nan_or_inf_data_bool_mask, get_svr_filter,
+                    merge_adata_attrs, unique_var_obs_adata)
 
 
 def deprecated(func):
diff --git a/dynamo/preprocessing/external/__init__.py b/dynamo/preprocessing/external/__init__.py
index bb73ed945..0d21bd4aa 100644
--- a/dynamo/preprocessing/external/__init__.py
+++ b/dynamo/preprocessing/external/__init__.py
@@ -1,8 +1,6 @@
 from .integration import harmony_debatch, integrate
-from .pearson_residual_recipe import (
-    normalize_layers_pearson_residuals,
-    select_genes_by_pearson_residuals,
-)
+from .pearson_residual_recipe import (normalize_layers_pearson_residuals,
+                                      select_genes_by_pearson_residuals)
 from .sctransform import sctransform
 
 __all__ = [
diff --git a/dynamo/preprocessing/external/pearson_residual_recipe.py b/dynamo/preprocessing/external/pearson_residual_recipe.py
index b4361306c..8dac00483 100644
--- a/dynamo/preprocessing/external/pearson_residual_recipe.py
+++ b/dynamo/preprocessing/external/pearson_residual_recipe.py
@@ -17,13 +17,10 @@
 from scipy.sparse import issparse
 
 from ...configuration import DKM
-from ...dynamo_logger import (
-    LoggerManager,
-    main_info,
-    main_info_insert_adata_layer,
-    main_warning,
-)
-from ...preprocessing.utils import is_nonnegative_integer_arr, seurat_get_mean_var
+from ...dynamo_logger import (LoggerManager, main_info,
+                              main_info_insert_adata_layer, main_warning)
+from ...preprocessing.utils import (is_nonnegative_integer_arr,
+                                    seurat_get_mean_var)
 from ..QC import filter_genes_by_outliers
 
 main_logger = LoggerManager.main_logger
@@ -109,7 +106,6 @@ def _highly_variable_pearson_residuals(
     # Get pearson residuals for each batch separately
     residual_gene_vars_by_batch = []
     for batch in np.unique(batch_info):
-
         adata_subset = adata[batch_info == batch]
 
         # Filter out zero genes
diff --git a/dynamo/preprocessing/external/sctransform.py b/dynamo/preprocessing/external/sctransform.py
index 4d85ba03c..9f9280b8d 100644
--- a/dynamo/preprocessing/external/sctransform.py
+++ b/dynamo/preprocessing/external/sctransform.py
@@ -20,7 +20,8 @@
 from scipy import stats
 
 from ...configuration import DKM
-from ...dynamo_logger import main_info, main_info_insert_adata_layer, main_warning
+from ...dynamo_logger import (main_info, main_info_insert_adata_layer,
+                              main_warning)
 from ..utils import get_gene_selection_filter
 
 _EPS = np.finfo(float).eps
diff --git a/dynamo/preprocessing/gene_selection.py b/dynamo/preprocessing/gene_selection.py
index 305b3208d..6253cb7b6 100644
--- a/dynamo/preprocessing/gene_selection.py
+++ b/dynamo/preprocessing/gene_selection.py
@@ -16,24 +16,13 @@
 from scipy.sparse import csr_matrix, issparse
 
 from ..configuration import DKM
-from ..dynamo_logger import (
-    LoggerManager,
-    main_critical,
-    main_debug,
-    main_info,
-    main_info_insert_adata_uns,
-    main_info_insert_adata_var,
-    main_warning,
-)
+from ..dynamo_logger import (LoggerManager, main_critical, main_debug,
+                             main_info, main_info_insert_adata_uns,
+                             main_info_insert_adata_var, main_warning)
 from .pca import pca
-from .utils import (
-    compute_gene_exp_fraction,
-    get_gene_selection_filter,
-    get_nan_or_inf_data_bool_mask,
-    get_svr_filter,
-    merge_adata_attrs,
-    seurat_get_mean_var,
-)
+from .utils import (compute_gene_exp_fraction, get_gene_selection_filter,
+                    get_nan_or_inf_data_bool_mask, get_svr_filter,
+                    merge_adata_attrs, seurat_get_mean_var)
 
 
 def calc_Gini(adata: AnnData, layers: Union[Literal["all"], List[str]] = "all") -> AnnData:
@@ -247,7 +236,11 @@ def calc_dispersion_by_svr(
         adata_ori.uns[key] = {"mean": mean, "cv": cv, "svr_gamma": svr_gamma}
 
         prefix = "" if layer == "X" else layer + "_"
-        (adata.var[prefix + "log_m"], adata.var[prefix + "log_cv"], adata.var[prefix + "score"],) = (
+        (
+            adata.var[prefix + "log_m"],
+            adata.var[prefix + "log_cv"],
+            adata.var[prefix + "score"],
+        ) = (
             np.nan,
             np.nan,
             -np.inf,
diff --git a/dynamo/preprocessing/normalization.py b/dynamo/preprocessing/normalization.py
index 245cda265..5054c3d67 100644
--- a/dynamo/preprocessing/normalization.py
+++ b/dynamo/preprocessing/normalization.py
@@ -14,12 +14,8 @@
 from scipy.sparse import issparse
 
 from ..configuration import DKM
-from ..dynamo_logger import (
-    main_debug,
-    main_info_insert_adata_layer,
-    main_info_insert_adata_obsm,
-    main_warning,
-)
+from ..dynamo_logger import (main_debug, main_info_insert_adata_layer,
+                             main_info_insert_adata_obsm, main_warning)
 from .utils import merge_adata_attrs
 
 
diff --git a/dynamo/preprocessing/pca.py b/dynamo/preprocessing/pca.py
index c3437d784..245b46240 100644
--- a/dynamo/preprocessing/pca.py
+++ b/dynamo/preprocessing/pca.py
@@ -16,7 +16,8 @@
 from sklearn.utils.sparsefuncs import mean_variance_axis
 
 from ..configuration import DKM
-from ..dynamo_logger import main_info_insert_adata_obsm, main_info_insert_adata_var
+from ..dynamo_logger import (main_info_insert_adata_obsm,
+                             main_info_insert_adata_var)
 
 
 def _truncatedSVD_with_center(
diff --git a/dynamo/simulation/__init__.py b/dynamo/simulation/__init__.py
index 23e1f1ff2..2175d2f10 100755
--- a/dynamo/simulation/__init__.py
+++ b/dynamo/simulation/__init__.py
@@ -1,25 +1,11 @@
 from .evaluation import evaluate
 from .Gillespie import Gillespie
-from .ODE import (
-    Simulator,
-    Ying_model,
-    jacobian_bifur2genes,
-    neurongenesis,
-    ode_bifur2genes,
-    state_space_sampler,
-    toggle,
-)
-from .simulate_anndata import (
-    AnnDataSimulator,
-    BifurcationTwoGenes,
-    CellularModelSimulator,
-    KinLabelingSimulator,
-    Neurongenesis,
-    OscillationTwoGenes,
-    bifur2genes_params,
-    bifur2genes_splicing_params,
-    neurongenesis_params,
-    osc2genes_params,
-    osc2genes_splicing_params,
-)
+from .ODE import (Simulator, Ying_model, jacobian_bifur2genes, neurongenesis,
+                  ode_bifur2genes, state_space_sampler, toggle)
+from .simulate_anndata import (AnnDataSimulator, BifurcationTwoGenes,
+                               CellularModelSimulator, KinLabelingSimulator,
+                               Neurongenesis, OscillationTwoGenes,
+                               bifur2genes_params, bifur2genes_splicing_params,
+                               neurongenesis_params, osc2genes_params,
+                               osc2genes_splicing_params)
 from .utils import CellularSpecies, directMethod
diff --git a/dynamo/simulation/simulate_anndata.py b/dynamo/simulation/simulate_anndata.py
index 0cd79af21..6d1a82fd5 100644
--- a/dynamo/simulation/simulate_anndata.py
+++ b/dynamo/simulation/simulate_anndata.py
@@ -5,24 +5,12 @@
 import pandas as pd
 
 # dynamo logger related
-from ..dynamo_logger import (
-    LoggerManager,
-    main_critical,
-    main_exception,
-    main_info,
-    main_tqdm,
-    main_warning,
-)
+from ..dynamo_logger import (LoggerManager, main_critical, main_exception,
+                             main_info, main_tqdm, main_warning)
 from ..tools.sampling import sample
 from ..tools.utils import flatten, isarray
-from .ODE import (
-    hill_act_func,
-    hill_inh_func,
-    neurongenesis,
-    ode_bifur2genes,
-    ode_neurongenesis,
-    ode_osc2genes,
-)
+from .ODE import (hill_act_func, hill_inh_func, neurongenesis, ode_bifur2genes,
+                  ode_neurongenesis, ode_osc2genes)
 from .utils import CellularSpecies, GillespieReactions, Reaction
 
 bifur2genes_params = {
@@ -248,7 +236,6 @@ def generate_anndata(self, remove_empty_cells: bool = False) -> anndata.AnnData:
             The AnnData object.
         """
         if self.T is not None and self.C is not None:
-
             obs = pd.DataFrame(
                 {
                     "cell_name": np.arange(self.get_n_cells()),
diff --git a/dynamo/simulation/utils.py b/dynamo/simulation/utils.py
index 328885c4a..6b0a567a6 100644
--- a/dynamo/simulation/utils.py
+++ b/dynamo/simulation/utils.py
@@ -3,14 +3,8 @@
 import numpy as np
 
 # dynamo logger related
-from ..dynamo_logger import (
-    LoggerManager,
-    main_critical,
-    main_exception,
-    main_info,
-    main_tqdm,
-    main_warning,
-)
+from ..dynamo_logger import (LoggerManager, main_critical, main_exception,
+                             main_info, main_tqdm, main_warning)
 
 
 def directMethod(