Repository: bioinf-mcb/Metagenomic-DeepFRI (language: python)

Issue:
Build fails 
Build fails with the following error: 
  ```
Building wheel for mDeepFRI (pyproject.toml) ... error
  error: subprocess-exited-with-error
  
  × Building wheel for mDeepFRI (pyproject.toml) did not run successfully.
  │ exit code: 1
  ╰─> [51 lines of output]
      /tmp/pip-build-env-svcl1p9q/overlay/lib/python3.10/site-packages/setuptools/dist.py:745: SetuptoolsDeprecationWarning: Invalid dash-separated options
      !!
      
              ********************************************************************************
              Usage of dash-separated 'description-file' will not be supported in future
              versions. Please use the underscore name 'description_file' instead.
      
              By 2023-Sep-26, you need to update your project and remove deprecated calls
              or your builds will no longer be supported.
      
              See https://setuptools.pypa.io/en/latest/userguide/declarative_config.html for details.
              ********************************************************************************
      
      !!
        opt = self.warn_dash_deprecation(opt, section)
      running bdist_wheel
      running build
      running build_py
      running egg_info
      writing mDeepFRI.egg-info/PKG-INFO
      writing dependency_links to mDeepFRI.egg-info/dependency_links.txt
      writing entry points to mDeepFRI.egg-info/entry_points.txt
      writing top-level names to mDeepFRI.egg-info/top_level.txt
      reading manifest file 'mDeepFRI.egg-info/SOURCES.txt'
      adding license file 'LICENSE'
      writing manifest file 'mDeepFRI.egg-info/SOURCES.txt'
      running build_ext
      cmake /home/ubuntu/Metagenomic-DeepFRI -DCMAKE_LIBRARY_OUTPUT_DIRECTORY=/home/ubuntu/Metagenomic-DeepFRI/build/lib.linux-x86_64-cpython-310/mDeepFRI/CPP_lib -DPY_INCLUDE_PATH=/home/ubuntu/miniconda3/envs/deepfri/include/python3.10/ -DCMAKE_BUILD_TYPE=Release
      -- Configuring done (0.0s)
      -- Generating done (0.0s)
      -- Build files have been written to: /home/ubuntu/Metagenomic-DeepFRI/build/temp.linux-x86_64-cpython-310
      cmake --build . --config Release -- -j4
      [ 25%] Building CXX object CMakeFiles/AtomDistanceIO.dir/mDeepFRI/CPP_lib/load_contact_maps.cpp.o
      /home/ubuntu/Metagenomic-DeepFRI/mDeepFRI/CPP_lib/load_contact_maps.cpp: In function ‘std::pair<bool*, int> LoadDenseContactMap(const string&, float)’:
      /home/ubuntu/Metagenomic-DeepFRI/mDeepFRI/CPP_lib/load_contact_maps.cpp:58:24: error: invalid initialization of reference of type ‘std::unique_ptr<float []>&’ from expression of type ‘float*’
         58 |           if (Distance(atoms_positions, atom_a, atom_b) <=
            |                        ^~~~~~~~~~~~~~~
      /home/ubuntu/Metagenomic-DeepFRI/mDeepFRI/CPP_lib/load_contact_maps.cpp:11:49: note: in passing argument 1 of ‘float Distance(std::unique_ptr<float []>&, int, int)’
         11 | static float Distance(std::unique_ptr<float[]> &array, int i, int j) {
            |                       ~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~
      /home/ubuntu/Metagenomic-DeepFRI/mDeepFRI/CPP_lib/load_contact_maps.cpp: In function ‘std::vector<std::pair<int, int> >* LoadSparseContactMap(const string&, float)’:
      /home/ubuntu/Metagenomic-DeepFRI/mDeepFRI/CPP_lib/load_contact_maps.cpp:106:24: error: invalid initialization of reference of type ‘std::unique_ptr<float []>&’ from expression of type ‘float*’
        106 |           if (Distance(atoms_positions, atom_a, atom_b) <=
            |                        ^~~~~~~~~~~~~~~
      /home/ubuntu/Metagenomic-DeepFRI/mDeepFRI/CPP_lib/load_contact_maps.cpp:11:49: note: in passing argument 1 of ‘float Distance(std::unique_ptr<float []>&, int, int)’
         11 | static float Distance(std::unique_ptr<float[]> &array, int i, int j) {
            |                       ~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~
      make[2]: *** [CMakeFiles/AtomDistanceIO.dir/build.make:104: CMakeFiles/AtomDistanceIO.dir/mDeepFRI/CPP_lib/load_contact_maps.cpp.o] Error 1
      make[1]: *** [CMakeFiles/Makefile2:83: CMakeFiles/AtomDistanceIO.dir/all] Error 2
      make: *** [Makefile:91: all] Error 2
      error: command '/home/ubuntu/miniconda3/envs/deepfri/bin/cmake' failed with exit code 2
      [end of output]
  
  note: This error originates from a subprocess, and is likely not a problem with pip.
  ERROR: Failed building wheel for mDeepFRI
Failed to build mDeepFRI
ERROR: Could not build wheels for mDeepFRI, which is required to install pyproject.toml-based projects
```


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 8fe3c9c..bf2a4f1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,11 +1,9 @@
 # project specific folder
 test_database/
-trained_models
-trained_models.tar.gz
-newest_trained_models.tar.gz
-biopython_test_output
-parasail_test_output
 notebooks
+mDeepFRI/*.cpp
+cython_tests.ipynb
+foldcomp
 
 docs/build/*
 build/
diff --git a/CMakeLists.txt b/CMakeLists.txt
deleted file mode 100644
index 09587b4..0000000
--- a/CMakeLists.txt
+++ /dev/null
@@ -1,20 +0,0 @@
-cmake_minimum_required(VERSION 3.19)
-project(GenomicLibrary)
-
-set(CMAKE_CXX_STANDARD 17)
-
-add_library(AtomDistanceIO SHARED
-        mDeepFRI/CPP_lib/library_definition.cpp
-        mDeepFRI/CPP_lib/atoms_file_io.h
-        mDeepFRI/CPP_lib/atoms_file_io.cpp
-        mDeepFRI/CPP_lib/python_utils.h
-        mDeepFRI/CPP_lib/load_contact_maps.h
-        mDeepFRI/CPP_lib/load_contact_maps.cpp)
-
-
-target_include_directories(AtomDistanceIO PUBLIC ${PY_INCLUDE_PATH})
-
-FIND_PACKAGE( Boost COMPONENTS python numpy REQUIRED )
-INCLUDE_DIRECTORIES( ${Boost_INCLUDE_DIR} )
-
-TARGET_LINK_LIBRARIES( AtomDistanceIO LINK_PUBLIC ${Boost_LIBRARIES} )
diff --git a/README.md b/README.md
index a19de08..8b1d804 100644
--- a/README.md
+++ b/README.md
@@ -13,8 +13,9 @@ This is the right project for this task! Pipeline in a nutshell:
 
 * [DeepFRI](https://github.com/SoliareofAstora/DeepFRI)
 * [MMseqs2](https://github.com/soedinglab/MMseqs2)
-* [Boost.Python](https://www.boost.org/doc/libs/1_75_0/libs/python/doc/html/index.html)
-* [Parasail](https://github.com/jeffdaily/parasail)
+* [FoldComp](https://github.com/steineggerlab/foldcomp)
+* [pyOpal](https://github.com/steineggerlab/foldcomp)
+* [ONNX](https://github.com/onnx/onnx)
 
 # Installation
 
@@ -41,53 +42,27 @@ mDeepFRI --help
 ```
 
 # Usage
+## Prepare structural database
+Download the database from the [website](https://foldcomp.steineggerlab.workers.dev/). The app was tested with `afdb_swissprot_v4`. You can use different databases, but be mindful that computation time might increase exponentially with the size of the database and the format of protein names might differ and app will crash.
 ## 1. Download models
 Run command:
 ```
 mDeepFRI get-models --output path/to/weights/folder
 ```
 
-## 2. Prepare database
-
-1. Upload structure (`.pdb` or `.mmcif`) files to a folder in your system.
-2. Run command:
+## 2. Predict protein function & capture log
 ```
-mDeepFRI build-db --input path/to/folder/with/strucures --output path/to/database -t threads
+mDeepFRI predict-function -i /path/to/protein/sequences -d /path/to/foldcomp/database/ -w /path/to/deepfri/weights/folder -o /output_path 2> log.txt
 ```
-**Tip:** building a database from AF2Swissprot (~550k predicted structures) on 32 CPU cores took ~30 min.
-
-Use parameter `-max_len` to define maximal length of the protein. Due to initial DeepFRI training set default value is set to `1000`.
-
-Main feature of this project is its ability to generate query contact map on the fly
-using results from mmseqs2 target database search for similar protein sequences with known structures.
-Later in the `metagenomic_deepfri.py` contact map alignment is performed to use it as input to DeepFRI GCN.
-(implemented in CPP_lib/load_contact_maps.h)
-
-The command will search for structure files,
-process them and store protein sequence and atoms positions inside `database/seq_atom_db`.
-It will also create a mmseqs2 database within `database/`.
-
-You can also use `--input DIR_1 FILE_2 ...` argument list to parse structures from multiple sources.
-Accepted formats are: `.pdb`, `.cif`, `.ent` both raw and compressed with `.gz`
-
-Protein ID is used as a filename. A new protein whose ID already exists in the database will be skipped.
-Use `--overwrite` flag to overwrite existing sequences and atoms positions.
-
-## 3. Predict protein function
-```
-mDeepFRI predict-function -i /path/to/protein/sequences -d /path/to/database/folder/from/previous/step -w /path/to/deepfri/weights/folder -o /output_path
-```
-**Attention:** Single instance of DeepFRI on GPU requires 10GB VRAM.
 
+Logging module writes output into `stderr`, so use `2>` to redirect it to file.
 Other available parameters can be found upon command `mDeepFRI --help`.
-
 ## Results
 Finished folder will contain:
-1. `query_files/*` - directory containing all input query files.
-2. `mmseqs2_search_results.m8`
-3. `alignments.json` - results of alignment search implemented in `utils.search_alignments.py`
-4. `metadata*` - files with some useful info
-5. `results*` - multiple files from DeepFRI. Organized by model type (`GCN` or `CNN`) and its mode (`mf`, `bp`, `cc`, `ec`) for the total of 8 files.
+1. `mmseqs2_search_results.m8`
+2. `metadata_skipped_ids_due_to_length` - too long or too short queries
+3. `queryDB` + index from MMSeqs2 search.
+4. `results*.tsv` - multiple files from DeepFRI. Organized by model type (`GCN` or `CNN`) and its mode (`mf`, `bp`, `cc`, `ec`).
 Sometimes results from one model can be missing which means that all query proteins sequences were aligned correctly or none of them were aligned.
    ```
    mf = molecular_function
@@ -96,15 +71,21 @@ Sometimes results from one model can be missing which means that all query prote
    ec = enzyme_commission
    ```
 
-## GPU / CPU utilization
-If CUDA is installed on your machine, `metaDeepFRI` will automatically use it for prediction, no additional installations are needed. If not, the model will use CPUs. If argument `threads` is provided, the prediction will run on multiple CPU cores.
+## Temporary files
+The first run of `mDeepFRI` with the database will create temporary files, needed for the pipeline. If you don't want to keep them for the next run use
+flag `--no-keep-temporary`.
+
+## CPU / GPU utilization
+If argument `threads` is provided, the app will parallelize certain steps (alignment, contact map alignment, inference).
+If CUDA is installed on your machine, `mDeepFRI` will automatically use it for prediction. If not, the model will use CPUs.
+**Attention:** Single instance of DeepFRI on GPU requires 2GB VRAM.
 
 ## Citations
 If you use this software please cite:
 - Gligorijević et al. "Structure-based protein function prediction using graph convolutional networks" Nat. Comms. (2021). https://doi.org/10.1038/s41467-021-23303-9
 - Steinegger & Söding "MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets" Nat. Biotechnol. (2017) https://doi.org/10.1038/nbt.3988
+- Kim, Midrita & Steinegger "Foldcomp: a library and format for compressing and indexing large protein structure sets" Bioinformatics (2023) https://doi.org/10.1093/bioinformatics/btad153
 - Maranga et al. "Comprehensive Functional Annotation of Metagenomes and Microbial Genomes Using a Deep Learning-Based Method" mSystems (2023) https://doi.org/10.1128/msystems.01178-22
-- Daily "Parasail: SIMD C library for global, semi-global, and local pairwise sequence alignments" BMC Bioinform. (2016) https://doi.org/10.1186/s12859-016-0930-z
 
 ## Contributing
 
@@ -112,5 +93,5 @@ If you have a suggestion that would make this project better, please send an e-m
 
 ### Contact
 
-Piotr Kucharski - soliareofastorauj@gmail.com \
-Valentyn Bezshapkin - valentyn.bezshapkin@micro.biol.ethz.ch
+Valentyn Bezshapkin - valentyn.bezshapkin@micro.biol.ethz.ch \
+Piotr Kucharski - soliareofastorauj@gmail.com
diff --git a/environment.yml b/environment.yml
index c093434..6e5445b 100644
--- a/environment.yml
+++ b/environment.yml
@@ -3,16 +3,19 @@ channels:
   - conda-forge
   - defaults
 dependencies:
-  - parasail-python=1.3.3
   - click
-  - boost>=1.80.0
-  - boost-cpp>=1.80.0
   - numpy>=1.19
-  - pandas>=1.2
+  - scipy
+  - biopython
   - mmseqs2=14.7e284
   - pip
-  - pysam>=0.16.0
   - pytest
-  - python
-  - cmake>=3.19
-  - onnxruntime>=1.10
+  - python>=3.11
+  - libdeflate
+  - pip:
+    - pyopal
+    - pysam>=0.21.0
+    - foldcomp
+    - onnxruntime-gpu
+    - torch
+    - pre-commit
diff --git a/mDeepFRI/CPP_lib/README.md b/mDeepFRI/CPP_lib/README.md
deleted file mode 100644
index a0cb37e..0000000
--- a/mDeepFRI/CPP_lib/README.md
+++ /dev/null
@@ -1,30 +0,0 @@
-# CPP_lib
-Here you can find a few performance critical functions used in metagenomic_deepfri_pipeline
-
-* `library_definition` contains functions definitions that are accessible in python using BOOST
-* In `atoms_file_io` you can find how protein structures are stored in binary format
-* `python_utils` implements quite interesting logic of [handling ownership of memory to python](https://stackoverflow.com/questions/57068443/setting-owner-in-boostpythonndarray-so-that-data-is-owned-and-managed-by-pyt)
-* `load_contact_maps` contains the most interesting functions
-
-### Build from source
-For more information please refer to the [cmake documentation](https://cmake.org/runningcmake/).
-
-First we need to install boost
-```bash
-wget https://boostorg.jfrog.io/artifactory/main/release/1.81.0/source/boost_1_81_0.tar.gz
-tar -xf boost_1_81_0.tar.gz
-cd boost_1_81_0
-./bootstrap.sh
-./b2
-```
-
-After building boost you will have to copy and paste include path.
-Additional libraries will be installed in `boost_1_81_0/stage/lib` directory.
-You will need to find and paste path to `libboost_python310.so.1.81.0`
-and `libboost_numpy310.so.1.81.0` files.
-
-```bash
-cd ../meta_deepFRI/CPP_lib
-ccmake .
-make
-```
diff --git a/mDeepFRI/CPP_lib/__init__.py b/mDeepFRI/CPP_lib/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/mDeepFRI/CPP_lib/atoms_file_io.cpp b/mDeepFRI/CPP_lib/atoms_file_io.cpp
deleted file mode 100644
index 1750a67..0000000
--- a/mDeepFRI/CPP_lib/atoms_file_io.cpp
+++ /dev/null
@@ -1,45 +0,0 @@
-//
-// Created by soliareofastora on 24.12.2021.
-//
-#include <boost/python.hpp>
-#include <boost/python/numpy.hpp>
-#include <fstream>
-
-namespace py = boost::python;
-namespace np = py::numpy;
-
-static void SaveAtomsFile(const np::ndarray &position_array,
-                          const np::ndarray &groups_array,
-                          const std::string &save_path) {
-  int chain_length = (int)groups_array.shape(0);
-  int *group_indexes = reinterpret_cast<int *>(groups_array.get_data());
-  int atom_count = group_indexes[chain_length - 1];
-  float *atom_positions = reinterpret_cast<float *>(position_array.get_data());
-
-  std::ofstream writer(save_path, std::ios::out | std::ios::binary);
-  writer.write(reinterpret_cast<const char *>(&chain_length), 4);
-  writer.write(reinterpret_cast<const char *>(group_indexes), 4 * chain_length);
-  writer.write(reinterpret_cast<const char *>(atom_positions),
-               4 * atom_count * 3);
-  writer.close();
-};
-
-static std::tuple<int, int *, float *>
-LoadAtomsFile(const std::string &file_path) {
-  std::ifstream reader(file_path, std::ios::in | std::ios::binary);
-
-  int chain_length;
-  reader.read(reinterpret_cast<char *>(&chain_length), 4);
-
-  int *group_indexes = new int[chain_length];
-  reader.read(reinterpret_cast<char *>(group_indexes), 4 * chain_length);
-
-  chain_length--;
-  int atom_count = group_indexes[chain_length];
-
-  float *atoms_positions = new float[atom_count * 3];
-  reader.read(reinterpret_cast<char *>(atoms_positions), 4 * atom_count * 3);
-  reader.close();
-
-  return std::make_tuple(chain_length, group_indexes, atoms_positions);
-};
diff --git a/mDeepFRI/CPP_lib/atoms_file_io.h b/mDeepFRI/CPP_lib/atoms_file_io.h
deleted file mode 100644
index fdbfb519..0000000
--- a/mDeepFRI/CPP_lib/atoms_file_io.h
+++ /dev/null
@@ -1,18 +0,0 @@
-#ifndef ATOMS_FILE_IO
-#define ATOMS_FILE_IO
-
-#include <boost/python.hpp>
-#include <boost/python/numpy.hpp>
-#include <fstream>
-
-namespace py = boost::python;
-namespace np = py::numpy;
-
-static void SaveAtomsFile(const np::ndarray &position_array,
-                          const np::ndarray &groups_array,
-                          const std::string &save_path);
-
-static std::tuple<int, int *, float *>
-LoadAtomsFile(const std::string &file_path);
-
-#endif
diff --git a/mDeepFRI/CPP_lib/library_definition.cpp b/mDeepFRI/CPP_lib/library_definition.cpp
deleted file mode 100644
index 536a2cf..0000000
--- a/mDeepFRI/CPP_lib/library_definition.cpp
+++ /dev/null
@@ -1,19 +0,0 @@
-//
-// Created by soliareofastora on 27.05.2021.
-//
-
-#include <boost/python.hpp>
-
-#include "atoms_file_io.h"
-#include "load_contact_maps.h"
-#include "python_utils.h"
-
-namespace py = boost::python;
-namespace np = py::numpy;
-
-BOOST_PYTHON_MODULE(libAtomDistanceIO) {
-  py::def("initialize", Initialize);
-  py::def("save_atoms", SaveAtomsFile);
-  py::def("load_contact_map", LoadContactMap);
-  py::def("load_aligned_contact_map", LoadAlignedContactMap);
-}
diff --git a/mDeepFRI/CPP_lib/load_contact_maps.cpp b/mDeepFRI/CPP_lib/load_contact_maps.cpp
deleted file mode 100644
index f33e178..0000000
--- a/mDeepFRI/CPP_lib/load_contact_maps.cpp
+++ /dev/null
@@ -1,199 +0,0 @@
-#include <cmath>
-#include <iostream>
-#include <queue>
-#include <string>
-#include <sys/stat.h>
-
-#include "atoms_file_io.h"
-#include "load_contact_maps.h"
-#include "python_utils.h"
-
-static float Distance(std::unique_ptr<float[]> &array, int i, int j) {
-  int xIndex = i * 3;
-  int yIndex = xIndex + 1;
-  int zIndex = xIndex + 2;
-
-  float xDist = array.get()[xIndex] - array.get()[xIndex + 3 * j];
-  float yDist = array.get()[yIndex] - array.get()[yIndex + 3 * j];
-  float zDist = array.get()[zIndex] - array.get()[zIndex + 3 * j];
-
-  float xSquared = xDist * xDist;
-  float ySquared = yDist * yDist;
-  float zSquared = zDist * zDist;
-
-  float distance = std::hypot(xSquared, ySquared, zSquared);
-
-  return distance;
-}
-
-inline bool exists(const std::string &name) {
-  struct stat buffer;
-  return (stat(name.c_str(), &buffer) == 0);
-}
-
-static std::pair<bool *, int>
-LoadDenseContactMap(const std::string &file_path,
-                    const float angstrom_contact_threshold) {
-  int chain_length;
-  int *group_indexes;
-  float *atoms_positions;
-  std::tie(chain_length, group_indexes, atoms_positions) =
-      LoadAtomsFile(file_path);
-
-  // allocate output array
-  bool *const output_data = new bool[(int)pow(chain_length, 2)];
-  std::memset(output_data, 0, (int)pow(chain_length, 2));
-
-  // fill up output array with atom contacts
-  for (int group_a = 0; group_a < chain_length; ++group_a) {
-    output_data[group_a * chain_length + group_a] = true;
-
-    for (int group_b = group_a + 1; group_b < chain_length; ++group_b) {
-      bool group_connected = false;
-
-      for (int atom_a = group_indexes[group_a];
-           atom_a < group_indexes[group_a + 1]; ++atom_a) {
-        for (int atom_b = group_indexes[group_b];
-             atom_b < group_indexes[group_b + 1]; ++atom_b) {
-          if (Distance(atoms_positions, atom_a, atom_b) <=
-              angstrom_contact_threshold) {
-            group_connected = true;
-            output_data[group_a * chain_length + group_b] = true;
-            output_data[group_a + group_b * chain_length] = true;
-            break;
-          }
-        }
-
-        if (group_connected)
-          break;
-      }
-    }
-  }
-
-  delete[] group_indexes;
-  delete[] atoms_positions;
-  return std::make_pair(output_data, chain_length);
-}
-
-static std::vector<std::pair<int, int>> *
-LoadSparseContactMap(const std::string &file_path,
-                     const float angstrom_contact_threshold) {
-  int chain_length;
-  int *group_indexes;
-  float *atoms_positions;
-
-  // verbose error if file_path doesn't exist
-  if (!exists(file_path)) {
-    std::cout << "Error: file " << file_path << " doesn't exist" << std::endl;
-    return nullptr;
-  }
-
-  std::tie(chain_length, group_indexes, atoms_positions) =
-      LoadAtomsFile(file_path);
-
-  // fill up vector with sparse atom contacts
-  auto *sparse_contacts = new std::vector<std::pair<int, int>>();
-  sparse_contacts->reserve(chain_length * 10);
-
-  for (int group_a = 0; group_a < chain_length; ++group_a) {
-    for (int group_b = group_a + 1; group_b < chain_length; ++group_b) {
-      bool group_connected = false;
-
-      for (int atom_a = group_indexes[group_a];
-           atom_a < group_indexes[group_a + 1]; ++atom_a) {
-        for (int atom_b = group_indexes[group_b];
-             atom_b < group_indexes[group_b + 1]; ++atom_b) {
-          if (Distance(atoms_positions, atom_a, atom_b) <=
-              angstrom_contact_threshold) {
-            group_connected = true;
-            sparse_contacts->emplace_back(group_a, group_b);
-            break;
-          }
-        }
-
-        if (group_connected)
-          break;
-      }
-    }
-  }
-
-  delete[] group_indexes;
-  delete[] atoms_positions;
-  return sparse_contacts;
-}
-
-static np::ndarray LoadContactMap(const std::string &file_path,
-                                  const float angstrom_contact_threshold) {
-  bool *contact_map;
-  int chain_length;
-  std::tie(contact_map, chain_length) =
-      LoadDenseContactMap(file_path, angstrom_contact_threshold);
-  return CreateNumpyArray(contact_map, chain_length);
-}
-
-static np::ndarray LoadAlignedContactMap(const std::string &file_path,
-                                         float angstrom_contact_threshold,
-                                         const std::string &query_alignment,
-                                         const std::string &target_alignment,
-                                         const int generated_contacts) {
-  std::vector<std::pair<int, int>> *sparse_target_contacts =
-      LoadSparseContactMap(file_path, angstrom_contact_threshold);
-  std::vector<std::pair<int, int>> sparse_query_contacts;
-  sparse_query_contacts.reserve(sparse_target_contacts->size());
-
-  int target_index = 0;
-  int query_index = 0;
-  std::vector<int> target_to_query_indexes;
-  target_to_query_indexes.reserve(query_alignment.size());
-
-  for (int i = 0; i < query_alignment.size(); ++i) {
-    if (query_alignment[i] == '-') {
-      target_to_query_indexes[target_index] = -1;
-      ++target_index;
-    } else if (target_alignment[i] == '-') {
-      for (int j = 1; j < generated_contacts + 1; ++j) {
-        sparse_query_contacts.emplace_back(query_index - j, query_index);
-        sparse_query_contacts.emplace_back(query_index + j, query_index);
-      }
-      ++query_index;
-    } else {
-      target_to_query_indexes[target_index] = query_index;
-      ++target_index;
-      ++query_index;
-    }
-  }
-
-  for (int i = 0; i < sparse_target_contacts->size(); ++i) {
-    int contact_x =
-        target_to_query_indexes[sparse_target_contacts->operator[](i).first];
-    if (contact_x < 0) {
-      continue;
-    }
-    int contact_y =
-        target_to_query_indexes[sparse_target_contacts->operator[](i).second];
-    if (contact_y < 0) {
-      continue;
-    }
-    sparse_query_contacts.emplace_back(contact_x, contact_y);
-  }
-
-  bool *const output_data = new bool[(int)pow(query_index, 2)];
-  std::memset(output_data, 0, (int)pow(query_index, 2));
-
-  for (int i = 0; i < query_index; ++i) {
-    output_data[i * query_index + i] = true;
-  }
-  for (std::pair<int, int> pair : sparse_query_contacts) {
-    if (pair.first < 0) {
-      continue;
-    }
-    if (pair.first >= query_index) {
-      continue;
-    }
-    output_data[pair.first * query_index + pair.second] = true;
-    output_data[pair.second * query_index + pair.first] = true;
-  }
-
-  delete sparse_target_contacts;
-  return CreateNumpyArray(output_data, query_index);
-}
diff --git a/mDeepFRI/CPP_lib/load_contact_maps.h b/mDeepFRI/CPP_lib/load_contact_maps.h
deleted file mode 100644
index b8e8c5c..0000000
--- a/mDeepFRI/CPP_lib/load_contact_maps.h
+++ /dev/null
@@ -1,34 +0,0 @@
-//
-// Created by soliareofastora on 04.01.2022.
-//
-
-#ifndef LOAD_CONTACT_MAPS
-#define LOAD_CONTACT_MAPS
-
-#include <cmath>
-#include <iostream>
-#include <queue>
-#include <string>
-#include <sys/stat.h>
-
-#include "atoms_file_io.h"
-#include "python_utils.h"
-
-static std::pair<bool *, int>
-LoadDenseContactMap(const std::string &file_path,
-                    const float angstrom_contact_threshold);
-
-static std::vector<std::pair<int, int>> *
-LoadSparseContactMap(const std::string &file_path,
-                     const float angstrom_contact_threshold);
-
-static np::ndarray LoadContactMap(const std::string &file_path,
-                                  const float angstrom_contact_threshold);
-
-static np::ndarray LoadAlignedContactMap(const std::string &file_path,
-                                         float angstrom_contact_threshold,
-                                         const std::string &query_alignment,
-                                         const std::string &target_alignment,
-                                         const int generated_contacts);
-
-#endif
diff --git a/mDeepFRI/CPP_lib/python_utils.h b/mDeepFRI/CPP_lib/python_utils.h
deleted file mode 100644
index 0379c46..0000000
--- a/mDeepFRI/CPP_lib/python_utils.h
+++ /dev/null
@@ -1,40 +0,0 @@
-//
-// Created by soliareofastora on 24.12.2021.
-//
-
-#ifndef PYTHON_UTILS
-#define PYTHON_UTILS
-
-#include <boost/python.hpp>
-#include <boost/python/numpy.hpp>
-
-namespace py = boost::python;
-namespace np = py::numpy;
-
-// Call this function before using other functions!!
-static void Initialize() {
-  Py_Initialize();
-  boost::python::numpy::initialize();
-}
-
-static void DestroyCapsule(PyObject* self) {
-  auto* b = reinterpret_cast<bool*>(PyCapsule_GetPointer(self, nullptr));
-  delete[] b;
-}
-
-static np::ndarray CreateNumpyArray(bool* array, int size){
-  // handling over memory management to python
-  // https://stackoverflow.com/questions/57068443/
-  //setting-owner-in-boostpythonndarray-so-that-data-is-owned-and-managed-by-python
-  PyObject* capsule = PyCapsule_New((void*) array, nullptr, (PyCapsule_Destructor) &DestroyCapsule);
-  py::handle<> capsule_handle{capsule};
-  py::object capsule_owner{capsule_handle};
-
-  return np::from_data(array,
-                       np::dtype::get_builtin<bool>(),
-                       py::make_tuple(size, size),
-                       py::make_tuple(sizeof(bool) * size, sizeof(bool)),
-                       capsule_owner);
-}
-
-#endif
diff --git a/mDeepFRI/__init__.py b/mDeepFRI/__init__.py
index 12c7e7c..56bb19d 100644
--- a/mDeepFRI/__init__.py
+++ b/mDeepFRI/__init__.py
@@ -1,45 +1,39 @@
-import os
-
 __version__ = "2.0.0"
-
-repo_link = "https://huggingface.co/crusher083/mDeepFRI/resolve/main"
-
-model_names = [
-    "DeepCNN-MERGED_bp.onnx",
-    "DeepCNN-MERGED_cc.onnx",
-    "DeepCNN-MERGED_mf.onnx",
-    "DeepCNN-MERGED_ec.onnx",
-    "DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_bp.onnx",
-    "DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_cc.onnx",
-    "DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_ec.onnx",
-    "DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_mf.onnx",
-]
-
-model_links = [
-    os.path.join(repo_link, model_name) for model_name in model_names
-]
-
-config_names = [
-    "DeepCNN-MERGED_bp_model_params.json",
-    "DeepCNN-MERGED_cc_model_params.json",
-    "DeepCNN-MERGED_mf_model_params.json",
-    "DeepCNN-MERGED_ec_model_params.json",
-    "DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_bp_model_params.json",
-    "DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_cc_model_params.json",
-    "DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_ec_model_params.json",
-    "DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_mf_model_params.json",
-    "model_config.json"
-]
+cnn_model_links = {
+    "bp":
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepCNN-MERGED_bp.onnx",
+    "cc":
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepCNN-MERGED_cc.onnx",
+    "mf":
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepCNN-MERGED_mf.onnx",
+    "ec":
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepCNN-MERGED_ec.onnx",
+}
+
+gcn_model_links = {
+    "bp":
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_bp.onnx",
+    "cc":
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_cc.onnx",
+    "mf":
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_mf.onnx",
+    "ec":
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_ec.onnx",
+}
 
 config_links = [
-    os.path.join(repo_link, config_name) for config_name in config_names
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepCNN-MERGED_bp_model_params.json",
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepCNN-MERGED_cc_model_params.json",
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepCNN-MERGED_mf_model_params.json",
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepCNN-MERGED_ec_model_params.json",
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_bp_model_params.json",
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_cc_model_params.json",
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_ec_model_params.json",
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/DeepFRI-MERGED_GraphConv_gcd_512-512-512_fcd_1024_ca_10.0_mf_model_params.json",
+    "https://huggingface.co/crusher083/mDeepFRI/resolve/main/model_config.json"
 ]
 
 TARGET_MMSEQS_DB_NAME = "targetDB"
-SEQUENCES = "seq"
-ATOMS = "atom"
-SEQ_ATOMS_DATASET_PATH = "seq_atom_db"
 
-ALIGNMENTS = "alignments.json"
-MERGED_SEQUENCES = 'merged_sequences.faa'
-MMSEQS_SEARCH_RESULTS = 'mmseqs2_search_results.m8'
+MERGED_SEQUENCES = "merged_sequences.faa"
+MMSEQS_SEARCH_RESULTS = "mmseqs2_search_results.m8"
diff --git a/mDeepFRI/alignment.py b/mDeepFRI/alignment.py
new file mode 100644
index 0000000..0a9678b
--- /dev/null
+++ b/mDeepFRI/alignment.py
@@ -0,0 +1,129 @@
+import logging
+from functools import partial
+from multiprocessing.pool import ThreadPool
+
+import pyopal
+
+from mDeepFRI.alignment_utils import alignment_identity
+
+
+def insert_gaps(sequence, reference, alignment_string):
+    sequence = list(sequence)
+    reference = list(reference)
+    alignment_string = list(alignment_string)
+
+    for i, a in enumerate(alignment_string):
+        if a == "I":
+            sequence.insert(i, "-")
+        elif a == "D":
+            reference.insert(i, "-")
+    return "".join(sequence), "".join(reference)
+
+
+class AlignmentResult:
+    def __init__(self,
+                 query_name,
+                 query_sequence,
+                 target_name,
+                 target_sequence,
+                 alignment,
+                 gapped_sequence="",
+                 gapped_target="",
+                 identity=None):
+        self.query_name = query_name
+        self.query_sequence = query_sequence
+        self.target_name = target_name
+        self.target_sequence = target_sequence
+        self.alignment = alignment
+        self.gapped_sequence = gapped_sequence
+        self.gapped_target = gapped_target
+        self.identity = identity
+
+    def __str__(self):
+        return f"AlignmentResult(query_name={self.query_name}, target_name={self.target_name}, " \
+               f"identity={self.identity})"
+
+    def __repr__(self):
+        return f"AlignmentResult(query_name={self.query_name}, target_name={self.target_name}, " \
+               f"identity={self.identity})"
+
+    def insert_gaps(self):
+        seq = list(self.query_sequence)
+        ref = list(self.target_sequence)
+        aln = list(self.alignment)
+
+        for i, a in enumerate(aln):
+            if a == "I":
+                seq.insert(i, "-")
+            elif a == "D":
+                ref.insert(i, "-")
+        self.gapped_sequence = "".join(seq)
+        self.gapped_target = "".join(ref)
+
+        return self
+
+    def calculate_identity(self):
+        self.identity = alignment_identity(self.gapped_sequence,
+                                           self.gapped_target)
+
+        return self
+
+
+def align_best_score(query_item: str, database: pyopal.Database,
+                     target_names: list, gap_open: int, gap_extend: int,
+                     identity_threshold: float):
+
+    query_name, query_sequence = query_item
+    results = database.search(query_sequence,
+                              mode="score",
+                              algorithm="nw",
+                              gap_open=gap_open,
+                              gap_extend=gap_extend)
+
+    best = max(results, key=lambda x: x.score)
+    best_sequence = database[best.target_index]
+    best_name = target_names[best.target_index]
+    single_seq = pyopal.Database([best_sequence])
+    alignment = single_seq.search(query_sequence,
+                                  mode="full",
+                                  algorithm="nw",
+                                  gap_open=gap_open,
+                                  gap_extend=gap_extend)[0].alignment
+
+    result = AlignmentResult(query_name, query_sequence, best_name,
+                             best_sequence, alignment)
+    result.insert_gaps().calculate_identity()
+
+    if result.identity >= identity_threshold:
+        return result
+
+    return None
+
+
+def align_query(query_seqs: dict,
+                target_seqs: dict,
+                alignment_gap_open: float,
+                alignment_gap_extend: float,
+                identity_threshold: str = 0.3,
+                threads: int = 1):
+
+    logging.info("Pairwise alignment started.")
+    logging.info("Aligning %i queries against %i database sequences.",
+                 len(query_seqs), len(target_seqs))
+    # Default substitution matrix BLOSUM50
+    database = pyopal.Database(list(target_seqs.values()))
+    align_against_db = partial(align_best_score,
+                               database=database,
+                               target_names=list(target_seqs.keys()),
+                               gap_open=alignment_gap_open,
+                               gap_extend=alignment_gap_extend,
+                               identity_threshold=identity_threshold)
+
+    with ThreadPool(threads) as pool:
+        all_alignments = pool.map(align_against_db, query_seqs.items())
+
+    filtered_alignments = list(filter(None, all_alignments))
+    logging.info("Found %i alignments.", len(filtered_alignments))
+    logging.info("Pairwise alignment finished.")
+
+    return filtered_alignments
diff --git a/mDeepFRI/alignment_utils.pyx b/mDeepFRI/alignment_utils.pyx
new file mode 100644
index 0000000..92bf488
--- /dev/null
+++ b/mDeepFRI/alignment_utils.pyx
@@ -0,0 +1,37 @@
+cimport cython
+from libc.stdlib cimport free, malloc
+from libc.string cimport strlen
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef float alignment_sequences_identity(char* query, char* target) nogil:
+    cdef int matches = 0
+    cdef int i
+    cdef float seq_identity
+    cdef int query_len = strlen(query)
+
+    # count matches excluding gaps
+    for i in range(query_len):
+        if query[i] == target[i] and query[i] != b'-' and target[i] != b'-':
+            matches += 1
+
+    # calculate identity
+    seq_identity = matches / query_len
+    return seq_identity
+
+
+def alignment_identity(str query, str target):
+    """
+    Calculate sequence identity between two proteins in an alignment.
+    Gaps are excluded. The length of sequence is irrelevant, because
+    lengths should be identical.
+
+    Args:
+        query (str): Query sequence.
+        target (str): Target sequence.
+
+    Returns:
+        float: Sequence identity between two sequences in an alignment.
+    """
+    return alignment_sequences_identity(query.encode(), target.encode())
diff --git a/mDeepFRI/bio_utils.py b/mDeepFRI/bio_utils.py
new file mode 100644
index 0000000..3045239
--- /dev/null
+++ b/mDeepFRI/bio_utils.py
@@ -0,0 +1,242 @@
+from io import StringIO
+from typing import Dict, List, Tuple
+
+import foldcomp
+import numpy as np
+from Bio.PDB import PDBParser
+from pysam import FastaFile, FastxFile
+from scipy.spatial.distance import pdist, squareform
+
+
+def load_fasta_as_dict(fasta_file):
+    """Load FASTA file as dict"""
+    with FastxFile(fasta_file) as fasta:
+        fasta_dict = {entry.name: entry.sequence for entry in fasta}
+    return fasta_dict
+
+
+def retrieve_fasta_entries_as_dict(fasta_file, entries):
+    """Retrieve selected FASTA entries as dict"""
+    fasta_dict = dict()
+    with FastaFile(fasta_file) as fasta:
+        for name in entries:
+            fasta_dict[name] = fasta.fetch(name)
+    return fasta_dict
+
+
+# copied from Biopython to remove dependency
+protein_letters_1to3 = {
+    "A": "Ala",
+    "C": "Cys",
+    "D": "Asp",
+    "E": "Glu",
+    "F": "Phe",
+    "G": "Gly",
+    "H": "His",
+    "I": "Ile",
+    "K": "Lys",
+    "L": "Leu",
+    "M": "Met",
+    "N": "Asn",
+    "P": "Pro",
+    "Q": "Gln",
+    "R": "Arg",
+    "S": "Ser",
+    "T": "Thr",
+    "V": "Val",
+    "W": "Trp",
+    "Y": "Tyr",
+}
+protein_letters_1to3_extended = {
+    **protein_letters_1to3,
+    **{
+        "B": "Asx",
+        "X": "Xaa",
+        "Z": "Glx",
+        "J": "Xle",
+        "U": "Sec",
+        "O": "Pyl"
+    },
+}
+
+protein_letters_3to1_extended = {
+    value: key
+    for key, value in protein_letters_1to3_extended.items()
+}
+
+PROTEIN_LETTERS = dict()
+for k, v in protein_letters_3to1_extended.items():
+    PROTEIN_LETTERS[str.upper(k)] = v
+PROTEIN_LETTERS["UNK"] = "X"
+
+
+def seq2onehot(seq):
+    """Create 26-dim embedding"""
+    chars = [
+        '-', 'D', 'G', 'U', 'L', 'N', 'T', 'K', 'H', 'Y', 'W', 'C', 'P', 'V',
+        'S', 'O', 'I', 'E', 'F', 'X', 'Q', 'A', 'B', 'Z', 'R', 'M'
+    ]
+    vocab_size = len(chars)
+    vocab_embed = dict(zip(chars, range(vocab_size)))
+
+    # Convert vocab to one-hot
+    vocab_one_hot = np.zeros((vocab_size, vocab_size), int)
+    for _, val in vocab_embed.items():
+        vocab_one_hot[val, val] = 1
+
+    embed_x = [vocab_embed[v] for v in seq]
+    seqs_x = np.array([vocab_one_hot[j, :] for j in embed_x])
+
+    return seqs_x
+
+
+def calculate_cmap(pdb_string, max_seq_len=1000, threshold=6.0, mode="matrix"):
+
+    parser = PDBParser()
+    structure = parser.get_structure("", StringIO(pdb_string))
+
+    residues = [x for x in structure.get_residues()][:max_seq_len]
+
+    coords = np.array([residue["CA"].get_coord() for residue in residues])
+    distances = squareform(pdist(coords))
+
+    cmap = (distances < threshold).astype(np.int32)
+
+    if mode == "sparse":
+        cmap = np.argwhere(cmap == 1).astype(np.uint32)
+    else:
+        pass
+
+    return cmap
+
+
+def retrieve_structure(idx, db):
+    with foldcomp.open(db, ids=[idx]) as db:
+        for _, pdb in db:
+            structure = pdb
+
+        return structure
+
+
+def align_contact_map(query_alignment: str,
+                      target_alignment: str,
+                      sparse_target_contact_map: List[Tuple[int, int]],
+                      generated_contacts: int = 2) -> np.ndarray:
+    """
+    Aligns a contact map based on the alignments of query and target sequences.
+
+    Args:
+        query_alignment: The alignment of the query sequence.
+        target_alignment: The alignment of the target sequence.
+        sparse_target_contact_map: The sparse contact map of the target
+                                   sequence represented as a list of tuples (i, j)
+                                   indicating contacts between residues iand j.
+        generated_contacts: The number of generated contacts to add for gapped
+                            regions in the query alignment. Defaults to 2.
+
+    Returns:
+        The aligned contact map as a numpy array.
+
+    Algorithm:
+    1. Initialize an empty list `sparse_query_contact_map` to store the contacts in the aligned contact map.
+    2. Initialize variables `target_index` and `query_index` to track the indices of residues in the target
+    and query proteins, respectively.
+    3. Initialize an empty dictionary `target_to_query_indices` to map target residues to query residues
+    using shift resulting from the alignments.
+    4. Iterate over each position in the query alignment:
+        - If the query residue is '-', increment the `target_index` and do not add a contact
+        to the aligned contact map.
+        - If the query residue is not '-', check the target residue:
+            - If the target residue is '-', add contacts for the generated region in the query alignment:
+                - For each generated contact, add the contact (query_index + j, query_index)
+                and (query_index - j, query_index) to the `sparse_query_contact_map` according to generated_contacts.
+                - Increment the `query_index`.
+            - If the target residue is not '-', map the target residue to the query residue by adding
+            an entry in the `target_to_query_indices` dictionary.
+                - Increment both the `query_index` and `target_index`.
+    5. Translate the target residue indices to query residue indices
+    in the `sparse_target_contact_map` by using the `target_to_query_indices` dictionary.
+    6. Filter out the contacts that are not present in the query alignment by removing contacts
+    with '-1' indices from the `sparse_target_contact_map`.
+    7. Add the filtered contacts from the filtered `sparse_target_contact_map` to the `sparse_query_contact_map`.
+    8. Build the output contact map with dimensions (query_index, query_index) initialized as all zeros.
+    Query index is the number of residues in the query sequence.
+    9. Set the diagonal elements of the output contact map to 1.
+    10. For each contact (i, j) in the `sparse_query_contact_map`:
+        - If i is less than 0 or greater than or equal to `query_index`, skip the contact.
+        - Otherwise, set the corresponding elements in the output contact map to 1 symmetrically.
+    11. Return the aligned contact map as a numpy array.
+    """
+    # The sparse contact map of the query sequence will contain all contacts. Will be used to create dense contact map
+    sparse_query_contact_map: List[Tuple[int, int]] = []
+
+    # The index of the residues in sequences
+    target_index: int = 0
+    query_index: int = 0
+
+    # Map target residues to query residues based on the alignments
+    target_to_query_indices: Dict[int, int] = {}
+
+    # Map target residues to query residues based on the alignments
+    for i in range(len(query_alignment)):
+        # If the query residue is a gap, skip target residue
+        if query_alignment[i] == "-":
+            target_to_query_indices[target_index] = -1
+            target_index += 1
+        else:
+            # If the target residue is a gap, add contacts to the query residue
+            # connected to generated_contacts nearest residues
+            if target_alignment[i] == "-":
+                for j in range(1, generated_contacts + 1):
+                    sparse_query_contact_map.append(
+                        (query_index + j, query_index))
+                    sparse_query_contact_map.append(
+                        (query_index - j, query_index))
+                query_index += 1
+            else:
+                # If there is an alignment match, map target residue to query residue
+                target_to_query_indices[target_index] = query_index
+                query_index += 1
+                target_index += 1
+
+    # Translate the target residues index to query residues index
+    sparse_map = list(
+        map(
+            lambda x:
+            (target_to_query_indices[x[0]], target_to_query_indices[x[1]]),
+            sparse_target_contact_map))
+    # Filter out the contacts that are not in the query alignment by removing columns and rows from sparse contact map
+    sparse_map = list(filter(lambda x: x[0] != -1 and x[1] != -1, sparse_map))
+    # Add the contacts to the output contact map
+    sparse_query_contact_map.extend(sparse_map)
+
+    # Build the output contact map
+    output_contact_map = np.zeros((query_index, query_index))
+    # Fill the diagonal
+    for i in range(query_index):
+        output_contact_map[i, i] = 1
+    # Fill the contacts from the sparse query contact map
+    for i, j in sparse_query_contact_map:
+        if i < 0:
+            continue
+        if i >= query_index:
+            continue
+        # Apply symmetry
+        output_contact_map[i, j] = 1
+        output_contact_map[j, i] = 1
+
+    return output_contact_map
+
+
+def retrieve_align_contact_map(alignment, database: str, max_seq_len: int,
+                               threshold: float, generated_contacts: int):
+    idx = alignment.target_name.rsplit(".", 1)[0]
+    pdb_string = retrieve_structure(idx, database)
+    cmap = calculate_cmap(pdb_string,
+                          max_seq_len=max_seq_len,
+                          threshold=threshold,
+                          mode="sparse")
+    aligned_cmap = align_contact_map(alignment.gapped_sequence,
+                                     alignment.gapped_target, cmap,
+                                     generated_contacts)
+    return (alignment, aligned_cmap)
diff --git a/mDeepFRI/cli.py b/mDeepFRI/cli.py
index 1d409e3..a8eca98 100644
--- a/mDeepFRI/cli.py
+++ b/mDeepFRI/cli.py
@@ -4,9 +4,8 @@
 import click
 
 from mDeepFRI import __version__
-from mDeepFRI.database import build_database
-from mDeepFRI.pipeline import metagenomic_deepfri
-from mDeepFRI.utils.utils import download_model_weights
+from mDeepFRI.pipeline import predict_protein_function
+from mDeepFRI.utils import download_model_weights
 
 logger = logging.getLogger(__name__)
 logger.setLevel(logging.INFO)
@@ -16,12 +15,11 @@
 @click.option("--debug/--no-debug", default=False)
 @click.version_option(version=__version__)
 @click.pass_context
-def cli(ctx, debug):
+def main(ctx, debug):
     """mDeepFRI"""
 
     ctx.ensure_object(dict)
     ctx.obj["debug"] = debug
-
     if debug:
         logger.setLevel(logging.DEBUG)
         logger.debug("Debug mode is on.")
@@ -32,7 +30,7 @@ def cli(ctx, debug):
         logging.getLogger("requests").setLevel(logging.INFO)
 
 
-@cli.command
+@main.command
 @click.option(
     "-o",
     "--output",
@@ -43,7 +41,6 @@ def cli(ctx, debug):
 @click.pass_context
 def get_models(ctx, output):
     """Download model weights for mDeepFRI."""
-
     if ctx.obj["debug"] is True:
         logger.setLevel(logging.DEBUG)
         logging.getLogger("requests").setLevel(logging.DEBUG)
@@ -53,60 +50,13 @@ def get_models(ctx, output):
         logging.getLogger("requests").setLevel(logging.INFO)
         logging.getLogger("urllib3").setLevel(logging.INFO)
 
+    logger.info("Downloading DeepFRI models.")
     output_path = Path(output)
     output_path.mkdir(parents=True, exist_ok=True)
     download_model_weights(output_path)
 
 
-# add subcommand for building database
-@cli.command()
-@click.option(
-    "-i",
-    "--input",
-    required=True,
-    multiple=True,
-    type=click.Path(exists=True),
-    help="Path to structure files or folders containing structure files.",
-)
-@click.option(
-    "-o",
-    "--output",
-    required=True,
-    type=click.Path(exists=False),
-    help="Path to folder where the database will be created.",
-)
-@click.option(
-    "-t",
-    "--threads",
-    default=1,
-    type=int,
-    help="Number of threads to use. Default is 1.",
-)
-@click.option(
-    "-m",
-    "--max-protein-length",
-    default=1_000,
-    type=int,
-    help=
-    "If protein is longer than this value, it will be truncated. Default is 1000 aa.",
-)
-@click.pass_context
-def build_db(ctx, input, output, threads, max_protein_length):
-    """Build a database for meta_deepFRI."""
-
-    if ctx.obj["debug"] is True:
-        logger.setLevel(logging.DEBUG)
-
-    input_seqs = [Path(seqs) for seqs in input]
-    output_path = Path(output)
-    build_database(input_paths=input_seqs,
-                   output_path=output_path,
-                   overwrite=True,
-                   threads=threads,
-                   max_protein_length=max_protein_length)
-
-
-@cli.command()
+@main.command()
 @click.option(
     "-i",
     "--input",
@@ -119,7 +69,7 @@ def build_db(ctx, input, output, threads, max_protein_length):
     "--db-path",
     required=True,
     type=click.Path(exists=True),
-    help="Path to a structures database folder.",
+    help="Path to a structures database compessed with FoldComp.",
 )
 @click.option(
     "-w",
@@ -135,20 +85,13 @@ def build_db(ctx, input, output, threads, max_protein_length):
     type=click.Path(exists=False),
     help="Path to output file.",
 )
-@click.option(
-    "-f",
-    "--output-format",
-    default="tsv",
-    type=click.Choice(["tsv", "csv", "json"]),
-    help="Output format. Default is TSV.",
-)
 @click.option(
     "-p",
     "--processing-modes",
     default=["bp", "cc", "ec", "mf"],
     type=click.Choice(["bp", "cc", "ec", "mf"]),
     multiple=True,
-    help="Processing modes. Default is all "
+    help="Processing modes. Default is all"
     "(biological process, cellular component, enzyme comission, molecular function).",
 )
 @click.option(
@@ -183,11 +126,11 @@ def build_db(ctx, input, output, threads, max_protein_length):
     help="Minimum identity for MMseqs2 alignment.",
 )
 @click.option(
-    "--alignment-matrix",
-    default="blosum62",
-    type=click.Choice(
-        ["blosum62", "pam250", "blosum80", "blosum45", "pam120", "pam160"]),
-    help="Alignment matrix for contact map alignment.",
+    "--top-k",
+    default=30,
+    type=int,
+    help="Number of top MMSeqs2 alignment for"
+    "precise pairwise alignment check. Default is 30.",
 )
 @click.option(
     "--alignment-gap-open",
@@ -207,6 +150,11 @@ def build_db(ctx, input, output, threads, max_protein_length):
     type=float,
     help="Minimum identity for contact map alignment.",
 )
+@click.option(
+    "--keep-intermediate/--no-keep-intermediate",
+    default=True,
+    help="Keep intermediate files. Default is True.",
+)
 @click.option(
     "-t",
     "--threads",
@@ -215,27 +163,38 @@ def build_db(ctx, input, output, threads, max_protein_length):
     help="Number of threads to use. Default is 1.",
 )
 @click.pass_context
-def predict_function(ctx, input, db_path, weights, output, output_format,
-                     processing_modes, angstrom_contact_thresh,
-                     generate_contacts, mmseqs_min_bit_score,
-                     mmseqs_max_evalue, mmseqs_min_identity, alignment_matrix,
-                     alignment_gap_open, alignment_gap_extend,
-                     alignment_min_identity, threads):
+def predict_function(ctx, input, db_path, weights, output, processing_modes,
+                     angstrom_contact_thresh, generate_contacts,
+                     mmseqs_min_bit_score, mmseqs_max_evalue,
+                     mmseqs_min_identity, top_k, alignment_gap_open,
+                     alignment_gap_extend, alignment_min_identity,
+                     keep_intermediate, threads):
     """Predict protein function from sequence."""
+    logger.info("Starting Metagenomic-DeepFRI.")
     if ctx.obj["debug"] is True:
         logger.setLevel(logging.DEBUG)
 
     output_path = Path(output)
     output_path.mkdir(parents=True, exist_ok=True)
 
-    metagenomic_deepfri(Path(input), Path(db_path), Path(weights), output_path,
-                        output_format, processing_modes,
-                        angstrom_contact_thresh, generate_contacts,
-                        mmseqs_min_bit_score, mmseqs_max_evalue,
-                        mmseqs_min_identity, alignment_matrix,
-                        alignment_gap_open, alignment_gap_extend,
-                        alignment_min_identity, threads)
+    predict_protein_function(
+        query_file=Path(input),
+        database=Path(db_path),
+        weights=Path(weights),
+        output_path=output_path,
+        deepfri_processing_modes=processing_modes,
+        angstrom_contact_threshold=angstrom_contact_thresh,
+        generate_contacts=generate_contacts,
+        mmseqs_min_bit_score=mmseqs_min_bit_score,
+        mmseqs_max_eval=mmseqs_max_evalue,
+        mmseqs_min_identity=mmseqs_min_identity,
+        top_k=top_k,
+        alignment_gap_open=alignment_gap_open,
+        alignment_gap_continuation=alignment_gap_extend,
+        identity_threshold=alignment_min_identity,
+        keep_intermediate=keep_intermediate,
+        threads=threads)
 
 
 if __name__ == "__main__":
-    cli(prog_name="mDeepFRI")
+    main()
diff --git a/mDeepFRI/database.py b/mDeepFRI/database.py
index bb22b8b..c961503 100644
--- a/mDeepFRI/database.py
+++ b/mDeepFRI/database.py
@@ -1,19 +1,9 @@
-import json
 # Create logger
 import logging
-import multiprocessing
-import pathlib
-from itertools import repeat
-from typing import List
+from pathlib import Path
 
-import numpy as np
-
-from mDeepFRI import SEQ_ATOMS_DATASET_PATH
-from mDeepFRI.CPP_lib import libAtomDistanceIO  # type: ignore[attr-defined]
-from mDeepFRI.structure_files.parse_structure_file import (
-    process_structure_file, search_structure_files)
-from mDeepFRI.utils.mmseqs import create_target_database
-from mDeepFRI.utils.utils import shutdown
+from mDeepFRI import MERGED_SEQUENCES, TARGET_MMSEQS_DB_NAME
+from mDeepFRI.mmseqs import create_target_database, extract_fasta_foldcomp
 
 logging.basicConfig(
     level=logging.DEBUG,
@@ -24,102 +14,44 @@
 
 
 def build_database(
-    input_paths: List[str],
+    input_path: str,
     output_path: str,
-    overwrite: bool,
     threads: int,
-    max_protein_length: int = 1_000,
 ) -> None:
     """
-    Searches for structure files in input (both files & folders are accepted, folders are globbed),
-    parses them, creates a database of available structures & sequences.
-    From sequences creates and indexes MMSeqs 2 database.
-    Additionally, creates a .json of added structures.
+    Extracts FASTA file from FoldComp database. Creates MMSeqs2 database and index.
 
     Args:
-        input_paths (list): list of paths to structure files or folders containing structure files.
-        output_path (str): path to folder where the database will be created.
-        overwrite (bool): flag to override existing sequences and atom positions.
+        input_path (str): path to a FoldComp database with compressed structures.
+        output_path (str): path to folder where the database for DeepFRI will be created.
         threads (int): number of threads to use.
-        max_protein_length (int): if protein is longer than this value, it will be truncated. Default is 1_000.
 
     Returns:
         None
     """
-    output_path = pathlib.Path(output_path)
-
-    seq_atoms_path = output_path / SEQ_ATOMS_DATASET_PATH
-    if not seq_atoms_path.exists():
-        seq_atoms_path.mkdir(parents=True)
-
-    # report max_protein_length
-    logger.info("MAX_PROTEIN_LENGTH: %s", max_protein_length)
-
-    logger.info("Searching for structure files in following paths:")
-    for input_path in input_paths:
-        logger.info("%s", str(input_path))
-    # dict [ structure_file_id ] = pathlib.Path( structure_file_path )
-    structure_files_paths = search_structure_files(input_paths)
 
-    n_structures = len(structure_files_paths)
-    if n_structures == 0:
-        shutdown("No structure files found")
+    logging.info("Building MMSeqs2 database from %s", input_path)
+    input_path = Path(input_path)
+    output_path = Path(output_path)
+    output_path.mkdir(parents=True, exist_ok=True)
+    output_sequences = output_path / MERGED_SEQUENCES
+    # check if files exist in output directory
+    if output_sequences.exists():
+        logging.info("Found %s in %s", MERGED_SEQUENCES, output_path)
+        logging.info(
+            "Skipping extraction of FASTA file from FoldComp database.")
     else:
-        logger.info("Found %s structure files to extract", n_structures)
-
-    # search for already processed protein_ids to skip them
-    if not overwrite:
-        duplicated_ids_counter = 0
-        # remove duplicates from structure_files_paths
-        for structure_id in structure_files_paths:
-            if (seq_atoms_path / structure_id).exists():
-                logger.info("Skipping %s - already in database", structure_id)
-                del structure_files_paths[structure_id]
-                duplicated_ids_counter += 1
-        logger.info("Found %s duplicated IDs", duplicated_ids_counter)
-
-    # actual processing of structure files takes place here
-    logger.info("Processing %s files", len(structure_files_paths))
-    libAtomDistanceIO.initialize()
-    with multiprocessing.Pool(processes=threads) as p:
-        processing_status = p.starmap(
-            process_structure_file,
-            zip(structure_files_paths.values(),
-                repeat(seq_atoms_path.absolute()), repeat(max_protein_length)))
-
-    # General stats over database
-    status, status_count = np.unique(processing_status, return_counts=True)
+        logging.info("Extracting FASTA file from FoldComp database.")
+        extract_fasta_foldcomp(input_path, output_sequences, threads)
+        logging.info("FASTA file extracted to %s", output_sequences)
 
-    # reversing the list, so that FAIL is the last one printed
-    for stat, count in list(zip(status, status_count))[::-1]:
-        logger.info("%s - %i", stat, count)
-
-    # Print out failed structures
-    for struct, status in zip(structure_files_paths.values(),
-                              processing_status):
-        if status.startswith("FAIL"):
-            logger.error("%s - %s", struct, status)
-
-    freshly_added_ids = []
-    structure_file_ids = list(structure_files_paths.keys())
-    for status, structure_id in zip(processing_status, structure_file_ids):
-        if status.startswith("SUCCESS"):
-            freshly_added_ids.append(structure_id)
-
-    if len(freshly_added_ids) == 0:
-        message = "\nNo new protein structures added.\n No new target database will be created."
-        shutdown(message)
-
-    # save db params
-    param_dict = {}
-    param_dict["sequences"] = sorted(freshly_added_ids)
-    param_dict["MAX_PROTEIN_LENGTH"] = max_protein_length
-    param_dict["input_structures_path"] = [
-        str(path) for path in sorted(input_paths)
-    ]
-
-    json.dump(param_dict,
-              open(output_path / "db_params.json", "w", encoding="utf-8"),
-              indent=4)
+    if (output_path / TARGET_MMSEQS_DB_NAME).exists():
+        logging.info("Found %s in %s", TARGET_MMSEQS_DB_NAME, output_path)
+        logging.info("Skipping creation of MMSeqs2 database.")
+    else:
+        logging.info("Creating MMSeqs2 database.")
+        create_target_database(output_sequences,
+                               output_path / TARGET_MMSEQS_DB_NAME)
+        logging.info("Database created at %s", output_path)
 
-    create_target_database(seq_atoms_path, output_path)
+    return (output_sequences, output_path / TARGET_MMSEQS_DB_NAME)
diff --git a/mDeepFRI/mmseqs.py b/mDeepFRI/mmseqs.py
new file mode 100644
index 0000000..3b33ae5
--- /dev/null
+++ b/mDeepFRI/mmseqs.py
@@ -0,0 +1,170 @@
+import logging
+import tempfile
+from functools import partial
+from multiprocessing.pool import ThreadPool
+from pathlib import Path
+
+import numpy as np
+
+import mDeepFRI
+from mDeepFRI import MMSEQS_SEARCH_RESULTS
+from mDeepFRI.utils import run_command
+
+MMSEQS_COLUMN_NAMES = [
+    "query", "target", "identity", "alignment_length", "mismatches",
+    "gap_openings", "query_start", "query_end", "target_start", "target_end",
+    "e_value", "bit_score"
+]
+
+logging.basicConfig(
+    level=logging.DEBUG,
+    format='[%(asctime)s] %(module)s.%(funcName)s %(levelname)s: %(message)s',
+    datefmt='%Y-%m-%d %H:%M:%S')
+
+logger = logging.getLogger(__name__)
+
+
+def createdb(sequences_file, db_path):
+    """
+    Converts FASTA file to a DB format needed for MMseqs2.
+    This should generate five files,
+    e.g. queryDB, queryDB_h and its corresponding index file queryDB.index,
+    queryDB_h.index and queryDB.lookup from the FASTA QUERY.fasta input sequences.
+
+    sequence_file (str): path to FASTA file.
+    db_path (str): path to output db file.
+
+    Returns:
+        None
+    """
+    run_command(f"mmseqs createdb {sequences_file} {db_path} --dbtype 1")
+
+
+def createindex(db_path):
+    with tempfile.TemporaryDirectory() as tmp_path:
+        run_command(f"mmseqs createindex {db_path} {tmp_path}")
+
+
+def search(query_db, target_db, result_db):
+    with tempfile.TemporaryDirectory() as tmp_path:
+        run_command(
+            f"mmseqs search {query_db} {target_db} {result_db} {tmp_path}")
+
+
+def convertalis(query_db, target_db, result_db, output_file):
+    run_command(
+        f"mmseqs convertalis {query_db} {target_db} {result_db} {output_file}")
+
+
+def extract_fasta_foldcomp(foldcomp_db: str,
+                           output_file: str,
+                           threads: int = 1):
+    """
+    Extracts FASTA from database
+    """
+    foldcomp_bin = Path(mDeepFRI.__path__[0]).parent / "foldcomp_bin"
+
+    # run command
+    run_command(
+        f"{foldcomp_bin} extract --fasta -t {threads} {foldcomp_db} {output_file}"
+    )
+
+
+def create_target_database(foldcomp_fasta_path: Path,
+                           mmseqs_db_path: Path) -> None:
+    """
+    Extracts sequences from compressed FoldComp database.
+
+    Args:
+        foldcomp_db_path (pathlib.Path): Path to FoldComp database.
+        new_db_path (pathlib.Path): Path to new MMSeqs database.
+
+    Returns:
+        None
+    """
+    createdb(foldcomp_fasta_path, mmseqs_db_path)
+    logging.info("Indexing new target mmseqs2 database %s",
+                 str(mmseqs_db_path))
+    createindex(mmseqs_db_path)
+
+
+def run_mmseqs_search(query_file: Path, target_db: Path,
+                      output_path: Path) -> Path:
+    """Creates a database from query sequences and runs mmseqs2 search against database.
+
+    Args:
+        query_file (pathlib.Path): Path to query FASTA file.
+        target_db (pathlib.Path): Path to target MMSeqs2 database.
+        output_path (pathlib.Path): Path to output folder.
+
+    Returns:
+        output_file (pathlib.Path): Path to MMSeqs2 search results.
+    """
+    output_path.mkdir(parents=True, exist_ok=True)
+    output_file = output_path / MMSEQS_SEARCH_RESULTS
+    query_db = output_path / 'queryDB'
+    createdb(query_file, query_db)
+
+    with tempfile.TemporaryDirectory() as tmp_path:
+
+        result_db = Path(tmp_path) / 'search_resultDB'
+        search(query_db, target_db, result_db)
+
+        # Convert results to tabular format
+        convertalis(query_db, target_db, result_db, output_file)
+
+    return output_file
+
+
+def filter_mmseqs_results(results_file: Path,
+                          min_bit_score: float = None,
+                          max_evalue: float = None,
+                          min_identity: float = None,
+                          k_best_hits: int = 30,
+                          threads: int = 1) -> np.recarray:
+    """
+    Filters MMSeqs results retrieving only k best hits based on identity
+    above specified thresholds. Allows number of paiwise alignments
+    in the next step of pipeline.
+
+    Args:
+        results_file (pathlib.Path): Path to MMSeqs2 search results.
+        min_bit_score (float): Minimum bit score.
+        max_evalue (float): Maximum e-value.
+        min_identity (float): Minimum identity.
+        k_best_hits (int): Number of best hits to keep.
+        threads (int): Number of threads to use.
+
+    Returns:
+        output (numpy.recarray): Filtered results.
+    """
+    def select_top_k(query, db, k=30):
+        return db[db["query"] == query][-k:]
+
+    output = np.recfromcsv(results_file,
+                           delimiter="\t",
+                           encoding="utf-8",
+                           names=MMSEQS_COLUMN_NAMES)
+
+    logging.info("%i MMSeqs2 hits in the database.", output.shape[0])
+
+    # MMSeqs2 alginment filters
+    if min_identity:
+        filtered = output[output['identity'] >= min_identity]
+    if min_bit_score:
+        filtered = output[output['bit_score'] >= min_bit_score]
+    if max_evalue:
+        filtered = output[output['e_value'] <= max_evalue]
+
+    # Get k best hits
+    filtered.sort(order=["query", "identity", "e_value"], kind="quciksort")
+    top_k_db = partial(select_top_k, db=filtered, k=k_best_hits)
+    with ThreadPool(threads) as pool:
+        top_k_chunks = pool.map(top_k_db, np.unique(filtered["query"]))
+
+    final_database = np.concatenate(top_k_chunks)
+
+    logging.info("%i pairs after filtering with k=%i best hits.",
+                 final_database.shape[0], k_best_hits)
+
+    return final_database
diff --git a/mDeepFRI/pipeline.py b/mDeepFRI/pipeline.py
index a48ff28..55b5f55 100644
--- a/mDeepFRI/pipeline.py
+++ b/mDeepFRI/pipeline.py
@@ -2,17 +2,20 @@
 import logging
 import os.path
 import pathlib
-from typing import List, Tuple
-
-from pysam.libcfaidx import FastxFile
-
-from mDeepFRI import ATOMS, SEQ_ATOMS_DATASET_PATH, TARGET_MMSEQS_DB_NAME
-from mDeepFRI.CPP_lib import libAtomDistanceIO  # type: ignore[attr-defined]
+from functools import partial
+from multiprocessing import Pool
+from typing import Dict, List
+
+import numpy as np
+
+from mDeepFRI import MERGED_SEQUENCES, TARGET_MMSEQS_DB_NAME
+from mDeepFRI.alignment import align_query
+from mDeepFRI.bio_utils import (load_fasta_as_dict, retrieve_align_contact_map,
+                                retrieve_fasta_entries_as_dict)
+from mDeepFRI.database import build_database
+from mDeepFRI.mmseqs import filter_mmseqs_results, run_mmseqs_search
 from mDeepFRI.predict import Predictor
-from mDeepFRI.utils.fasta_file_io import SeqFileLoader
-from mDeepFRI.utils.mmseqs import run_mmseqs_search
-from mDeepFRI.utils.search_alignments import search_alignments
-from mDeepFRI.utils.utils import load_deepfri_config
+from mDeepFRI.utils import remove_temporary
 
 logging.basicConfig(
     level=logging.DEBUG,
@@ -21,49 +24,25 @@
 
 logger = logging.getLogger(__name__)
 
-###########################################################################
-# in a nutshell:
-#
-#   load_and_verify_job_data
-#   1.  select first .faa file inside task_path
-#   2.  filter out proteins that are too long
-#   3.  find target database
-#
-#   metagenomic_deepfri
-#   4.  run mmseqs2 search on query and target database
-#   5.  find the best alignment for pairs found by mmseqs2 search.
-#   6.  If alignment for query exists:
-#           DeepFRI GCN for query sequence with aligned target contact map
-#       else:
-#           DeepFRI CNN for query sequence alone
-###########################################################################
-
-
-def check_inputs(
-    query_file: pathlib.Path, database: pathlib.Path, output_path: pathlib.Path
-) -> Tuple[pathlib.Path, dict, pathlib.Path, SeqFileLoader]:
+
+def load_query_sequences(query_file, output_path) -> Dict[str, str]:
     """
-    Check if input files and directories exist and are valid. Filters out proteins that are too long
-    or too short.
+    Loads query protein sequences from FASTA file. Filters
+    out sequences that are too short or too long.
 
     Args:
-        query_file (pathlib.Path): Path to a query file with protein sequences.
-        database (pathlib.Path): Path to a directory with a pre-built database.
-        output_path (pathlib.Path): Path to a directory where results will be saved.
-
-    Raises:
-        ValueError: Query file does not contain parsable protein sequences.
-        FileNotFoundError: MMSeqs2 database appears to be corrupted.
+        query_file (str): Path to FASTA file with query protein sequences.
+        output_path (str): Path to output folder.
 
     Returns:
-        Tuple[pathlib.Path, dict, pathlib.Path, SeqFileLoader]: Tuple of query file path, query sequences,
-            path to target MMSeqs2 database and target sequences.
+        query_seqs (dict): Dictionary with query protein headers to sequences.
     """
 
+    # By DeepFRI design
     MIN_PROTEIN_LENGTH = 60
+    MAX_PROTEIN_LENGTH = 1_000
 
-    with FastxFile(query_file) as fasta:
-        query_seqs = {record.name: record.sequence for record in fasta}
+    query_seqs = load_fasta_as_dict(query_file)
 
     if len(query_seqs) == 0:
         raise ValueError(
@@ -72,9 +51,6 @@ def check_inputs(
     logging.info("Found total of %i protein sequences in %s", len(query_seqs),
                  query_file)
 
-    with open(database / "db_params.json", "r", encoding="utf-8") as f:
-        MAX_PROTEIN_LENGTH = json.load(f)["MAX_PROTEIN_LENGTH"]
-
     prot_len_outliers = {}
     for prot_id, sequence in query_seqs.items():
         prot_len = len(sequence)
@@ -89,17 +65,35 @@ def check_inputs(
             "Skipping %i proteins due to sequence length outside range %i-%i aa.",
             len(prot_len_outliers), MIN_PROTEIN_LENGTH, MAX_PROTEIN_LENGTH)
         logging.info("Skipped protein ids will be saved in " \
-                     "metadata_skipped_ids_length.json")
+                     "metadata_skipped_ids_length.json.")
         json.dump(prot_len_outliers,
                   open(output_path / 'metadata_skipped_ids_due_to_length.json',
                        "w",
                        encoding="utf-8"),
                   indent=4,
                   sort_keys=True)
-        if len(query_seqs) == 0:
-            logging.info(
-                "All sequences in %s were too long. No sequences will be processed.",
-                query_file)
+
+    assert len(query_seqs
+               ) > 0, "All proteins were filtered out due to sequence length."
+
+    return query_seqs
+
+
+def check_mmseqs_database(database: pathlib.Path):
+    """
+    Check if MMSeqs2 database is intact.
+
+    Args:
+        query_file (pathlib.Path): Path to a query file with protein sequences.
+        database (pathlib.Path): Path to a directory with a pre-built database.
+        output_path (pathlib.Path): Path to a directory where results will be saved.
+
+    Raises:
+        FileNotFoundError: MMSeqs2 database appears to be corrupted.
+
+    Returns:
+        target_db (pathlib.Path): Path to MMSeqs2 database.
+    """
 
     # Verify all the files for MMSeqs2 database
     mmseqs2_ext = [
@@ -115,14 +109,12 @@ def check_inputs(
         raise FileNotFoundError(
             "MMSeqs2 database appears to be corrupted. Please, rebuild it.")
 
-    target_seqs = SeqFileLoader(database / SEQ_ATOMS_DATASET_PATH)
-
-    return query_file, query_seqs, target_db, target_seqs
+    return target_db
 
 
-def check_deepfri_weights(weights: pathlib.Path) -> pathlib.Path:
+def load_deepfri_config(weights: pathlib.Path) -> pathlib.Path:
     """
-    Check if DeepFRI weights are valid.
+    Check if DeepFRI weights are valid and load config.
     Args:
         weights :
 
@@ -150,72 +142,73 @@ def check_deepfri_weights(weights: pathlib.Path) -> pathlib.Path:
             ), f"DeepFRI weights are missing {model_type} model at {model_name}"
             assert config_name.exists(
             ), f"DeepFRI weights are missing {model_type} model config at {config_name}"
-
-    return config_path
+            # correct config
+            models_config[net]["models"][model_type] = str(
+                model_name.absolute())
+    return models_config
 
 
 ## TODO: structure output folder
-def metagenomic_deepfri(query_file: pathlib.Path, database: pathlib.Path,
-                        weights: pathlib.Path, output_path: pathlib.Path,
-                        output_format: List[str],
-                        deepfri_processing_modes: List[str],
-                        angstrom_contact_threshold: float,
-                        generate_contacts: int, mmseqs_min_bit_score: float,
-                        mmseqs_max_eval: float, mmseqs_min_identity: float,
-                        alignment_matrix: str, alignment_gap_open: float,
-                        alignment_gap_continuation: float,
-                        alignment_min_identity: float, threads: int):
+def predict_protein_function(
+        query_file: str,
+        database: str,
+        weights: str,
+        output_path: str,
+        deepfri_processing_modes: List[str] = ["ec", "bp", "mf", "cc"],
+        angstrom_contact_threshold: float = 6,
+        generate_contacts: int = 2,
+        mmseqs_min_bit_score: float = None,
+        mmseqs_max_eval: float = 10e-5,
+        mmseqs_min_identity: float = 0.3,
+        top_k: int = 30,
+        alignment_gap_open: float = 10,
+        alignment_gap_continuation: float = 1,
+        identity_threshold: float = 0.3,
+        keep_intermediate=True,
+        threads: int = 1):
     """
-    Run metagenomic-DeepFRI.
-    Args:
-        query_file:
-        database:
-        weights:
-        output_path:
-        output_format:
-        deepfri_processing_modes:
-        angstrom_contact_threshold:
-        generate_contacts:
-        mmseqs_min_bit_score:
-        mmseqs_max_eval:
-        mmseqs_min_identity:
-        alignment_match:
-        alignment_missmatch:
-        alignment_gap_open:
-        alignment_gap_continuation:
-        alignment_min_identity:
-        threads:
 
-    Returns:
 
     """
-    logging.info("Starting metagenomic-DeepFRI.")
-    model_config_json = check_deepfri_weights(weights)
-
-    query_file, query_seqs, target_db, target_seqs = check_inputs(
-        query_file, database, output_path)
-
-    logging.info("Running metagenomic-DeepFRI for %i sequences",
-                 len(query_seqs))
-    logging.info("Running MMSeqs2 search for the query against database")
-    mmseqs_search_output = run_mmseqs_search(query_file, target_db,
-                                             output_path, mmseqs_min_bit_score,
-                                             mmseqs_max_eval,
-                                             mmseqs_min_identity)
-
-    mmseqs2_found_proteins = mmseqs_search_output.shape[0]
-    logging.info("Found %i proteins in the database", mmseqs2_found_proteins)
-
-    alignments = search_alignments(query_seqs, mmseqs_search_output,
-                                   target_seqs, output_path, alignment_matrix,
-                                   alignment_gap_open,
-                                   alignment_gap_continuation,
-                                   alignment_min_identity, threads)
-
-    unaligned_queries = query_seqs.keys() - alignments.keys()
-
-    libAtomDistanceIO.initialize()
-    deepfri_models_config = load_deepfri_config(model_config_json)
+    MAX_SEQ_LEN = 1000
+    query_file = pathlib.Path(query_file)
+    database = pathlib.Path(database)
+    intermediate = build_database(
+        input_path=database,
+        output_path=database.parent,
+        threads=threads,
+    )
+
+    weights = pathlib.Path(weights)
+    output_path = pathlib.Path(output_path)
+    output_path.mkdir(parents=True, exist_ok=True)
+
+    query_seqs = load_query_sequences(query_file, output_path)
+    logging.info("Aligning %i sequences with MMSeqs2.", len(query_seqs))
+    target_db = check_mmseqs_database(database.parent)
+
+    mmseqs_results = run_mmseqs_search(query_file, target_db, output_path)
+    filtered_mmseqs_results = filter_mmseqs_results(mmseqs_results,
+                                                    mmseqs_min_bit_score,
+                                                    mmseqs_max_eval,
+                                                    mmseqs_min_identity,
+                                                    k_best_hits=top_k,
+                                                    threads=threads)
+
+    target_ids = np.unique(filtered_mmseqs_results["target"]).tolist()
+    target_seqs = retrieve_fasta_entries_as_dict(
+        database.parent / MERGED_SEQUENCES, target_ids)
+
+    alignments = align_query(query_seqs, target_seqs, alignment_gap_open,
+                             alignment_gap_continuation, identity_threshold,
+                             threads)
+
+    aligned_queries = [aln.query_name for aln in alignments]
+    unaligned_queries = {
+        k: v
+        for k, v in query_seqs.items() if k not in aligned_queries
+    }
+    deepfri_models_config = load_deepfri_config(weights)
 
     # deepfri_processing_modes = ['mf', 'bp', 'cc', 'ec']
     # mf = molecular_function
@@ -223,44 +216,39 @@ def metagenomic_deepfri(query_file: pathlib.Path, database: pathlib.Path,
     # cc = cellular_component
     # ec = enzyme_commission
 
+    gcn_prots, cnn_prots = len(aligned_queries), len(unaligned_queries)
+
+    if gcn_prots > 0:
+        partial_align = partial(retrieve_align_contact_map,
+                                database=database,
+                                max_seq_len=MAX_SEQ_LEN,
+                                threshold=angstrom_contact_threshold,
+                                generated_contacts=generate_contacts)
+        logging.info("Aligning contact maps for %i proteins", gcn_prots)
+        with Pool(threads) as p:
+            aligned_cmaps = p.map(partial_align, alignments)
+
+        logging.info("Aligned %i contact maps", len(aligned_cmaps))
+
     for mode in deepfri_processing_modes:
         logging.info("Processing mode: %s", mode)
         # GCN for queries with aligned contact map
-        gcn_prots, cnn_prots = len(alignments), len(unaligned_queries)
 
-        if gcn_prots > 0:
-            logging.info("Predicting with GCN: %i proteins", gcn_prots)
-            output_file_name = output_path / f"results_gcn_{mode}"
+        logging.info("Predicting with GCN: %i proteins", gcn_prots)
+        output_file_name = output_path / f"results_gcn_{mode}"
+        gcn_params = deepfri_models_config["gcn"]["models"][mode]
 
-            gcn_params = deepfri_models_config["gcn"]["models"][mode]
-            gcn = Predictor(gcn_params, gcn=True, threads=threads)
+        gcn = Predictor(gcn_params, threads=threads)
 
-            for query_id, alignment in alignments.items():
-                logging.info("Predicting %s", query_id)
-                query_seq = query_seqs[query_id]
-                target_id = alignment["target_id"]
-
-                generated_query_contact_map = libAtomDistanceIO.load_aligned_contact_map(
-                    str(database / SEQ_ATOMS_DATASET_PATH / ATOMS /
-                        (target_id + ".bin")),
-                    angstrom_contact_threshold,
-                    alignment["query_sequence"],  # query alignment
-                    alignment["target_sequence"],  # target alignment
-                    generate_contacts)
-
-                # running the actual prediction
-                gcn.predict_function(seqres=query_seq,
-                                     cmap=generated_query_contact_map,
-                                     chain=query_id)
-
-                if "tsv" in output_format:
-                    gcn.export_tsv(output_file_name.with_suffix('.tsv'))
-                if "csv" in output_format:
-                    gcn.export_csv(output_file_name.with_suffix('.csv'))
-                if "json" in output_format:
-                    gcn.export_json(output_file_name.with_suffix('.json'))
+        for aln, aligned_cmap in aligned_cmaps:
+            logging.info("Predicting %s", aln.query_name)
+            # running the actual prediction
+            gcn.predict_function(seqres=aln.query_sequence,
+                                 cmap=aligned_cmap,
+                                 chain=aln.query_name)
 
-            del gcn
+        gcn.export_tsv(str(output_file_name.with_suffix('.tsv')))
+        del gcn
 
         # CNN for queries without satisfying alignments
         if cnn_prots > 0:
@@ -268,19 +256,16 @@ def metagenomic_deepfri(query_file: pathlib.Path, database: pathlib.Path,
             output_file_name = output_path / f"results_cnn_{mode}"
 
             cnn_params = deepfri_models_config["cnn"]["models"][mode]
-            cnn = Predictor(cnn_params, gcn=False, threads=threads)
+            cnn = Predictor(cnn_params, threads=threads)
             for query_id in unaligned_queries:
                 logging.info("Predicting %s", query_id)
                 cnn.predict_function(seqres=query_seqs[query_id],
                                      chain=query_id)
 
-            if "tsv" in output_format:
-                cnn.export_tsv(output_file_name.with_suffix('.tsv'))
-            if "csv" in output_format:
-                cnn.export_csv(output_file_name.with_suffix('.csv'))
-            if "json" in output_format:
-                cnn.export_json(output_file_name.with_suffix('.json'))
-
+            cnn.export_tsv(str(output_file_name.with_suffix('.tsv')))
             del cnn
 
+    if not keep_intermediate:
+        remove_temporary(intermediate)
+
     logging.info("meta-DeepFRI finished successfully")
diff --git a/mDeepFRI/predict.py b/mDeepFRI/predict.pyx
similarity index 70%
rename from mDeepFRI/predict.py
rename to mDeepFRI/predict.pyx
index b324873..516fe87 100644
--- a/mDeepFRI/predict.py
+++ b/mDeepFRI/predict.pyx
@@ -2,25 +2,41 @@
 import json
 
 import numpy as np
+
+cimport numpy as np
+
+np.import_array()
+
 import onnxruntime as rt
+import torch
 
-from mDeepFRI.utils.bio_utils import seq2onehot
+from mDeepFRI.bio_utils import seq2onehot
 
 
 ## TODO: do not accumulate predictions, write them out to csv asap to avoid unneccessary RAM usage
-class Predictor(object):
+cdef class Predictor(object):
     """
     Class for loading trained models and computing GO/EC predictions and class activation maps (CAMs).
     """
-    def __init__(self, model_prefix: str, gcn: bool = True, threads: int = 0):
-        self.model_prefix = model_prefix
-        self.gcn = gcn
+
+    cdef public bint gcn
+    cdef public str model_path
+    cdef public int threads
+    cdef public dict prot2goterms
+    cdef public dict goidx2chains
+    cdef public np.ndarray gonames
+    cdef public np.ndarray goterms
+    cdef public np.ndarray thresh
+    cdef public session
+    cdef public np.ndarray Y_hat
+    cdef public dict data
+    cdef public list test_prot_list
+
+    def __init__(self, model_path: str, threads: int = 0):
+        self.model_path = model_path
+        self.threads = threads
 
         # Not clear how parameter influences GPU exec
-        if rt.get_device() == 'CPU':
-            self.threads = threads
-        elif rt.get_device() == 'GPU':
-            self.threads = 0
 
         self._load_model()
         self.prot2goterms = {}
@@ -28,15 +44,21 @@ def __init__(self, model_prefix: str, gcn: bool = True, threads: int = 0):
 
     def _load_model(self):
         session_options = rt.SessionOptions()
-        session_options.intra_op_num_threads = self.threads
+
+        if rt.get_device() == 'CPU':
+            session_options.intra_op_num_threads = self.threads
+        elif rt.get_device() == 'GPU':
+            session_options.intra_op_num_threads = 1
+
         self.session = rt.InferenceSession(
-            self.model_prefix + '.onnx',
-            providers=['CUDAExecutionProvider', 'CPUExecutionProvider'],
+            self.model_path,
+            providers=[('CUDAExecutionProvider', {"cudnn_conv_algo_search": "DEFAULT"}),
+                        'CPUExecutionProvider'],
             sess_options=session_options,
         )
 
         # load parameters
-        with open(self.model_prefix + "_model_params.json") as json_file:
+        with open(self.model_path.rsplit(".", 1)[0] + "_model_params.json") as json_file:
             metadata = json.load(json_file)
 
         self.gonames = np.asarray(metadata['gonames'])
@@ -45,8 +67,8 @@ def _load_model(self):
 
     def predict_function(
         self,
-        seqres: np.ndarray,
-        cmap: np.ndarray = None,
+        seqres: str,
+        cmap = None,
         chain: str = "",
     ):
         """
@@ -60,13 +82,20 @@ def predict_function(
         Returns:
             None
         """
+
+        cdef np.ndarray A
+        cdef np.ndarray prediction
+        cdef np.ndarray y
+
+
         self.Y_hat = np.zeros((1, len(self.goterms)), dtype=float)
         self.data = {}
         self.test_prot_list = [chain]
         S = seq2onehot(seqres)
         S = S.reshape(1, *S.shape)
         inputDetails = self.session.get_inputs()
-        if self.gcn:
+        # if cmap present use GCN with 2 inputs - sequence + cmap
+        if cmap is not None:
             A = cmap.reshape(1, *cmap.shape)
             prediction = self.session.run(
                 None, {
@@ -74,6 +103,8 @@ def predict_function(
                     inputDetails[1].name: S.astype(np.float32)
                 })[0]
             self.data[chain] = [[A, S], seqres]
+
+        # if no cmap use CNN with 1 input
         else:
             prediction = self.session.run(
                 None, {inputDetails[0].name: S.astype(np.float32)})[0]
@@ -90,29 +121,6 @@ def predict_function(
             self.prot2goterms[chain].append(
                 (self.goterms[idx], self.gonames[idx], float(y[idx])))
 
-    def export_json(self, output_fn: str):
-        """
-        Exports predictions to .json format
-
-        Args:
-            output_fn (str): output file name
-
-        Returns:
-            None
-        """
-        data = []
-        for prot, goterms in self.prot2goterms.items():
-            sorted_rows = sorted(goterms, key=lambda x: x[2], reverse=True)
-            for row in sorted_rows:
-                data.append({
-                    'Protein': prot,
-                    'GO_term/EC_number': row[0],
-                    'Score': '{:.5f}'.format(row[2]),
-                    'GO_term/EC_number name': row[1]
-                })
-
-        with open(output_fn, 'w', encoding="utf-8") as f:
-            json.dump(data, f, indent=4)
 
     def __save_file(self,
                     output_fn: str,
@@ -142,13 +150,6 @@ def __save_file(self,
                     writer.writerow(
                         [prot, row[0], '{:.5f}'.format(row[2]), row[1]])
 
-    def export_csv(self, output_fn: str):
-        """
-        Exports predictions to .csv format
-        """
-
-        self.__save_file(output_fn, ",")
-
     def export_tsv(self, output_fn: str):
         """
         Exports predictions to .tsv format
diff --git a/mDeepFRI/structure_files/__init__.py b/mDeepFRI/structure_files/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/mDeepFRI/structure_files/parse_structure_file.py b/mDeepFRI/structure_files/parse_structure_file.py
deleted file mode 100644
index 1607fc6..0000000
--- a/mDeepFRI/structure_files/parse_structure_file.py
+++ /dev/null
@@ -1,190 +0,0 @@
-import dataclasses
-import gzip
-import logging
-import pathlib
-import traceback
-
-import numpy as np
-
-from mDeepFRI import ATOMS, SEQUENCES
-from mDeepFRI.CPP_lib import libAtomDistanceIO  # type: ignore[attr-defined]
-from mDeepFRI.structure_files.parsers import parse_mmcif, parse_pdb
-from mDeepFRI.utils import bio_utils
-
-# need to parse different type of files? Add a file name pattern with a parser in this dict.
-# read structure_files_parsers/README.md for more information about how to create new parser.
-PARSERS = {
-    '.pdb': parse_pdb,
-    '.pdb.gz': parse_pdb,
-    '.cif': parse_mmcif,
-    '.cif.gz': parse_mmcif,
-    '.ent': parse_pdb,
-    '.ent.gz': parse_pdb
-}
-
-
-@dataclasses.dataclass
-class SeqAtoms:
-    protein_id: str
-    atom_amino_group: np.ndarray
-    positions: np.ndarray
-    groups: np.ndarray
-
-
-def search_structure_files(input_paths: list):
-    """
-
-    :param input_paths:
-    :return:
-    """
-    structure_files_paths = dict()
-    # Iterate input_paths and check if file extension matches any key in structure_files.PARSERS
-    for input_path in input_paths:
-        print(f"{str(input_path)}")
-        if input_path.is_file():
-            for pattern in PARSERS:
-                if str(input_path).endswith(pattern):
-                    print(f"\tFile {input_path} is {pattern}")
-                    structure_file_id = input_path.name[:-len(pattern)]
-                    structure_files_paths[structure_file_id] = input_path
-                    continue
-        elif input_path.is_dir():
-            files_inside_input_path = list(input_path.glob("**/*"))
-            for pattern in PARSERS:
-                structure_file_paths = list(
-                    filter(lambda x: str(x).endswith(pattern),
-                           files_inside_input_path))
-                if len(structure_file_paths) == 0:
-                    continue
-                print(
-                    f"\tFound {len(structure_file_paths)} {pattern} files in {input_path}"
-                )
-                structure_file_ids = [
-                    x.name[:-len(pattern)] for x in structure_file_paths
-                ]
-                structure_files_paths.update(
-                    zip(structure_file_ids, structure_file_paths))
-        else:
-            print(f"\tUnable to find {str(input_path)}")
-    return structure_files_paths
-
-
-def read_structure_file(file_path: pathlib.Path) -> SeqAtoms:
-    """
-    Extract sequence and atom positions from structure file
-    :param file_path:
-    :return:
-    """
-    for pattern, structure_parsing_func in PARSERS.items():
-        if file_path.name.endswith(pattern):
-            if file_path.name.endswith('.gz'):
-                f = gzip.open(file_path, 'rt')
-            else:
-                f = open(file_path, 'r')
-            atom_amino_group, positions, groups = structure_parsing_func(f)
-            f.close()
-
-            protein_id = file_path.name.replace(pattern, '')
-
-    return SeqAtoms(protein_id, atom_amino_group, positions, groups)
-
-
-def save_sequence_and_atoms(seq_atoms: SeqAtoms, sequence_path: pathlib.Path,
-                            atoms_path: pathlib.Path,
-                            max_target_chain_length: int) -> str:
-    """
-    Reads the structure file extracting sequence and atom positions.
-    If the sequence is too short, it is skipped as it is considered to be corrupted or DNA sequence.
-    If the sequence is too long, it is truncated based on the max_target_chain_length defined in target_db_config.json.
-
-    Saves the extracted data:
-        sequence - protein_id.faa file containing single sequence f.write(f">{protein_id}\n{sequence}\n")
-            SEQ_ATOMS_DATASET_PATH / SEQUENCES / (protein_id + ".faa")
-        atom positions - binary file containing positions of all atoms and correlated amino acid chain index
-            SEQ_ATOMS_DATASET_PATH / ATOMS / (protein_id + ".bin")
-            For more information on how those binary are saved check out source code at CPP_lib/atoms_file_io.h
-
-    Args:
-        seq_atoms: An object of SeqAtoms containing the sequence and atom positions.
-        sequence_path: A pathlib.Path object containing the path to the sequence file to be saved.
-        atoms_path: A pathlib.Path object containing the path to the atom positions file to be saved.
-        max_target_chain_length: An integer value of the maximum length of the protein chain.
-
-    Returns:
-        If the protein chain is too short, it returns "Sequence is too short, probably DNA or corrupted".
-        If the protein chain is truncated, it returns the string with the new length of the chain.
-        Otherwise, it returns "SUCCEED".
-    """
-    _, groups_index = np.unique(seq_atoms.groups, return_index=True)
-    groups_index.sort()
-
-    if len(groups_index) < 9:
-        # Files containing protein chains shorter than 9 amino acids might be corrupted or contain DNA
-        return "FAIL - Sequence is too short, probably DNA or corrupted"
-
-    truncated = False
-    if len(groups_index) > max_target_chain_length:
-        # Protein chains with more than {max_target_chain_length} are truncated.
-        # Change CONFIG/RUNTIME_PARAMETERS.py :max_target_chain_length
-        # and rerun this script with --overwrite to process this file again
-        truncated = True
-        group_indexes = groups_index[:max_target_chain_length]
-        group_indexes = np.append(
-            group_indexes,
-            groups_index[max_target_chain_length]).astype(np.int32)
-    else:
-        group_indexes = np.append(
-            groups_index, seq_atoms.positions.shape[0]).astype(np.int32)
-
-    sequence = ''.join([
-        bio_utils.PROTEIN_LETTERS[seq_atoms.atom_amino_group[i]]
-        for i in group_indexes[:-1]
-    ])
-    with open(sequence_path, "w", encoding="utf-8") as f:
-        f.write(f">{seq_atoms.protein_id}\n{sequence}\n")
-
-    libAtomDistanceIO.save_atoms(seq_atoms.positions, group_indexes,
-                                 str(atoms_path))
-
-    if truncated:
-        return f"SUCCESS, but sequences and contact maps got truncated to {max_target_chain_length}"
-    else:
-        return "SUCCESS"
-
-
-def process_structure_file(structure_file, save_path, max_target_chain_length):
-    """
-
-    :param structure_file:
-    :param save_path:
-    :param max_target_chain_length:
-    :return:
-    """
-    try:
-        seq_atoms = read_structure_file(structure_file)
-    except Exception:
-        print("EXCEPTION WHILE READING FILE ", str(structure_file))
-        logging.error(traceback.format_exc())
-        return "file reading exceptions"
-
-    # create directories for saving if they don't exist
-    sequences_folder = save_path / SEQUENCES
-    atoms_path = save_path / ATOMS
-    sequences_folder.mkdir(parents=True, exist_ok=True)
-    atoms_path.mkdir(parents=True, exist_ok=True)
-
-    # process single protein
-    prot_sequence_path = save_path / SEQUENCES / (seq_atoms.protein_id +
-                                                  ".faa")
-    prot_atoms_path = save_path / ATOMS / (seq_atoms.protein_id + ".bin")
-
-    try:
-        return save_sequence_and_atoms(seq_atoms, prot_sequence_path,
-                                       prot_atoms_path,
-                                       max_target_chain_length)
-    except Exception:
-        print("EXCEPTION DURING FILE PROCESSING ", str(structure_file))
-        logging.error(traceback.format_exc())
-        prot_sequence_path.unlink(missing_ok=True)
-        prot_atoms_path.unlink(missing_ok=True)
-        return "file processing exceptions"
diff --git a/mDeepFRI/structure_files/parsers.py b/mDeepFRI/structure_files/parsers.py
deleted file mode 100644
index 3be468c..0000000
--- a/mDeepFRI/structure_files/parsers.py
+++ /dev/null
@@ -1,125 +0,0 @@
-from io import TextIOWrapper
-from typing import Tuple
-
-import numpy as np
-
-
-def parse_pdb(
-        file: TextIOWrapper) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
-    """Parse a PDB file.
-
-    Args:
-        file (TextIOWrapper): Open file instance.
-
-    Returns:
-        sequence (list): Aminoacid sequence.
-    """
-    sequence = []
-    positions = []
-    groups = []
-    line = file.readline()
-    while line != "":
-        if line.startswith("TER"):
-            break
-        if line.startswith("ATOM"):
-            if len(line) == 81:
-                if line[76] != 'H' and line[17] != ' ':
-                    sequence.append(line[17:20])
-                    positions.append([line[30:38], line[38:46], line[46:54]])
-                    groups.append(line[21:26])
-        line = file.readline()
-
-    sequence = np.array(sequence)
-    positions = np.array(positions, dtype=np.float32)
-    groups = np.array(groups)
-
-    return sequence, positions, groups
-
-
-def parse_mmcif(
-        file: TextIOWrapper) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
-    """Parse a mmCIF file.
-
-    Args:
-        file (TextIOWrapper): Open file instance.
-
-    Returns:
-        sequence (np.array): Aminoacid sequence.
-        positions (np.array): Positions of the atoms.
-        groups (np.array): Groups of the atoms.
-    """
-
-    n_atoms = -1
-
-    label_counter = -1
-    line = file.readline()
-    while line != "":
-        if line.startswith('_refine_hist.pdbx_number_atoms_protein'):
-            n = line.split()[-1]
-            if n != '?' and n != '0':
-                try:
-                    n_atoms = int(n)
-                except ValueError:
-                    line = file.readline()
-                    continue
-                sequence = np.empty(n_atoms, dtype=np.dtype('<U3'))
-                positions = np.empty((n_atoms, 3), dtype=np.float32)
-                groups = np.empty(n_atoms, dtype=np.dtype('<U10'))
-
-        if line.startswith('_atom_site.'):
-            label_counter += 1
-            if line.startswith('_atom_site.type_symbol '):
-                atom_symbol = label_counter
-            if line.startswith('_atom_site.label_asym_id '):
-                assembly = label_counter
-            if line.startswith('_atom_site.label_seq_id '):
-                sequence_id = label_counter
-            if line.startswith('_atom_site.label_comp_id '):
-                residue = label_counter
-            if line.startswith('_atom_site.Cartn_x '):
-                x = label_counter
-            if line.startswith('_atom_site.Cartn_y '):
-                y = label_counter
-            if line.startswith('_atom_site.Cartn_z '):
-                z = label_counter
-
-        line = file.readline()
-        if line.startswith("ATOM"):
-            if label_counter > 0:
-                break
-
-    if n_atoms > 0:
-        index = 0
-        while line != '':
-            if line.startswith('loop_'):
-                break
-            if line.startswith('ATOM'):
-                atom = line.split()
-                if len(atom[residue]) == 3 and atom[atom_symbol] != "H":
-                    sequence[index] = atom[residue]
-                    groups[index] = ''.join(
-                        [atom[assembly], atom[sequence_id]])
-                    positions[index][0] = atom[x]
-                    positions[index][1] = atom[y]
-                    positions[index][2] = atom[z]
-                    index += 1
-                    if index == n_atoms:
-                        break
-            line = file.readline()
-        sequence.resize(index)
-        positions.resize((index, 3))
-        groups.resize(index)
-    else:
-        while line != '':
-            if line.startswith('loop_'):
-                break
-            if line.startswith('ATOM'):
-                atom = line.split()
-                if len(atom[residue]) == 3 and atom[atom_symbol] != "H":
-                    sequence.append(atom[residue])
-                    groups.append(''.join([atom[assembly], atom[sequence_id]]))
-                    positions.append([atom[x], atom[y], atom[z]])
-            line = file.readline()
-        positions = np.array(positions, dtype=np.float32)
-
-    return sequence, positions, groups
diff --git a/mDeepFRI/utils/utils.py b/mDeepFRI/utils.py
similarity index 75%
rename from mDeepFRI/utils/utils.py
rename to mDeepFRI/utils.py
index fd22d47..218fe5f 100644
--- a/mDeepFRI/utils/utils.py
+++ b/mDeepFRI/utils.py
@@ -1,14 +1,16 @@
-import json
 import logging
 import pathlib
 import shlex
 import shutil
 import subprocess
 import sys
+from glob import glob
+from pathlib import Path
+from typing import Iterable
 
 import requests
 
-from mDeepFRI import config_links, model_links
+from mDeepFRI import cnn_model_links, config_links, gcn_model_links
 
 
 def run_command(command, timeout=-1):
@@ -47,6 +49,7 @@ def download_file(url, path):
 
 def download_model_weights(output_path: pathlib.Path):
 
+    model_links = list(cnn_model_links.items()) + list(gcn_model_links.items())
     total_len = len(model_links)
     for i, model_link in enumerate(model_links):
         download_file(model_link, output_path / model_link.split("/")[-1])
@@ -93,24 +96,9 @@ def shutdown(message):
     sys.exit(message)
 
 
-def load_deepfri_config(filepath_model_config_json: str) -> dict:
-    """Loads model_config.json with paths to different models.
-
-    Args:
-        filepath_model_config_json (str): path to a file within trained_models folder.
-        Distributed with original DeepFRI repo.
-
-    Returns:
-        dict: a dict of different models and paths to their weights.
-    """
-    json_filepath = pathlib.Path(filepath_model_config_json)
-    if not json_filepath.exists():
-        raise FileNotFoundError(f"Config file not found at {json_filepath}.")
-
-    # load and replace local paths to files with absolute paths
-    with open(filepath_model_config_json, "r") as json_file:
-        json_string = json_file.read().replace("./trained_models",
-                                               f"{json_filepath.parent}")
-    deepfri_config = json.loads(json_string)
-
-    return deepfri_config
+def remove_temporary(temporary_files: Iterable):
+    for file in temporary_files:
+        extensions = glob(str(file) + "*")
+        for ext in extensions:
+            logging.info(f"Removing temporary file {ext}.")
+            Path(ext).unlink()
diff --git a/mDeepFRI/utils/__init__.py b/mDeepFRI/utils/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/mDeepFRI/utils/bio_utils.py b/mDeepFRI/utils/bio_utils.py
deleted file mode 100644
index b91e5fd..0000000
--- a/mDeepFRI/utils/bio_utils.py
+++ /dev/null
@@ -1,137 +0,0 @@
-import numpy as np
-import parasail
-
-# copied from Biopython to remove dependency
-protein_letters_1to3 = {
-    "A": "Ala",
-    "C": "Cys",
-    "D": "Asp",
-    "E": "Glu",
-    "F": "Phe",
-    "G": "Gly",
-    "H": "His",
-    "I": "Ile",
-    "K": "Lys",
-    "L": "Leu",
-    "M": "Met",
-    "N": "Asn",
-    "P": "Pro",
-    "Q": "Gln",
-    "R": "Arg",
-    "S": "Ser",
-    "T": "Thr",
-    "V": "Val",
-    "W": "Trp",
-    "Y": "Tyr",
-}
-protein_letters_1to3_extended = {
-    **protein_letters_1to3,
-    **{
-        "B": "Asx",
-        "X": "Xaa",
-        "Z": "Glx",
-        "J": "Xle",
-        "U": "Sec",
-        "O": "Pyl"
-    },
-}
-
-protein_letters_3to1_extended = {
-    value: key
-    for key, value in protein_letters_1to3_extended.items()
-}
-
-PROTEIN_LETTERS = dict()
-for k, v in protein_letters_3to1_extended.items():
-    PROTEIN_LETTERS[str.upper(k)] = v
-PROTEIN_LETTERS["UNK"] = "X"
-
-
-def seq2onehot(seq):
-    """Create 26-dim embedding"""
-    chars = [
-        '-', 'D', 'G', 'U', 'L', 'N', 'T', 'K', 'H', 'Y', 'W', 'C', 'P', 'V',
-        'S', 'O', 'I', 'E', 'F', 'X', 'Q', 'A', 'B', 'Z', 'R', 'M'
-    ]
-    vocab_size = len(chars)
-    vocab_embed = dict(zip(chars, range(vocab_size)))
-
-    # Convert vocab to one-hot
-    vocab_one_hot = np.zeros((vocab_size, vocab_size), int)
-    for _, val in vocab_embed.items():
-        vocab_one_hot[val, val] = 1
-
-    embed_x = [vocab_embed[v] for v in seq]
-    seqs_x = np.array([vocab_one_hot[j, :] for j in embed_x])
-
-    return seqs_x
-
-
-# parasail martices
-substitution_matrices = {
-    "blosum30": parasail.blosum30,
-    "blosum35": parasail.blosum35,
-    "blosum40": parasail.blosum40,
-    "blosum45": parasail.blosum45,
-    "blosum50": parasail.blosum50,
-    "blosum55": parasail.blosum55,
-    "blosum60": parasail.blosum60,
-    "blosum62": parasail.blosum62,
-    "blosum65": parasail.blosum65,
-    "blosum70": parasail.blosum70,
-    "blosum75": parasail.blosum75,
-    "blosum80": parasail.blosum80,
-    "blosum85": parasail.blosum85,
-    "blosum90": parasail.blosum90,
-    "blosum100": parasail.blosum100,
-    "pam10": parasail.pam10,
-    "pam20": parasail.pam20,
-    "pam30": parasail.pam30,
-    "pam40": parasail.pam40,
-    "pam50": parasail.pam50,
-    "pam60": parasail.pam60,
-    "pam70": parasail.pam70,
-    "pam80": parasail.pam80,
-    "pam90": parasail.pam90,
-    "pam100": parasail.pam100,
-    "pam110": parasail.pam110,
-    "pam120": parasail.pam120,
-    "pam130": parasail.pam130,
-    "pam140": parasail.pam140,
-    "pam150": parasail.pam150,
-    "pam160": parasail.pam160,
-    "pam170": parasail.pam170,
-    "pam180": parasail.pam180,
-    "pam190": parasail.pam190,
-    "pam200": parasail.pam200,
-    "pam210": parasail.pam210,
-    "pam220": parasail.pam220,
-    "pam230": parasail.pam230,
-    "pam240": parasail.pam240,
-    "pam250": parasail.pam250,
-    "pam260": parasail.pam260,
-    "pam270": parasail.pam270,
-    "pam280": parasail.pam280,
-    "pam290": parasail.pam290,
-    "pam300": parasail.pam300,
-    "pam310": parasail.pam310,
-    "pam320": parasail.pam320,
-    "pam330": parasail.pam330,
-    "pam340": parasail.pam340,
-    "pam350": parasail.pam350,
-    "pam360": parasail.pam360,
-    "pam370": parasail.pam370,
-    "pam380": parasail.pam380,
-    "pam390": parasail.pam390,
-    "pam400": parasail.pam400,
-    "pam410": parasail.pam410,
-    "pam420": parasail.pam420,
-    "pam430": parasail.pam430,
-    "pam440": parasail.pam440,
-    "pam450": parasail.pam450,
-    "pam460": parasail.pam460,
-    "pam470": parasail.pam470,
-    "pam480": parasail.pam480,
-    "pam490": parasail.pam490,
-    "pam500": parasail.pam500
-}
diff --git a/mDeepFRI/utils/fasta_file_io.py b/mDeepFRI/utils/fasta_file_io.py
deleted file mode 100644
index e004fc8..0000000
--- a/mDeepFRI/utils/fasta_file_io.py
+++ /dev/null
@@ -1,25 +0,0 @@
-import pathlib
-
-from pysam.libcfaidx import FastxFile
-
-from mDeepFRI import ATOMS, SEQUENCES
-
-
-class SeqFileLoader:
-    def __init__(self, path):
-        self.path = pathlib.Path(path)
-
-    def __getitem__(self, target_id):
-        protein_file = (self.path / SEQUENCES / (target_id + ".faa"))
-        structure_file = (self.path / ATOMS / (target_id + ".bin"))
-        assert protein_file.exists(), f"Target database {self.path.name} does not contain sequence" \
-                                      f"{self.path / SEQUENCES / (target_id + '.faa')}. " \
-                                      "It should not happen. Please create new target database."
-        assert structure_file.exists(), f"Target database {self.path.name} does not contain structure" \
-                                        f"{self.path / ATOMS / (target_id + '.bin')}. " \
-                                        "It should not happen. Please create new target database."
-
-        with FastxFile(self.path / SEQUENCES / (target_id + ".faa"), "r") as f:
-            sequence = [record.sequence for record in f][0]
-
-        return sequence
diff --git a/mDeepFRI/utils/mmseqs.py b/mDeepFRI/utils/mmseqs.py
deleted file mode 100644
index cd47c70..0000000
--- a/mDeepFRI/utils/mmseqs.py
+++ /dev/null
@@ -1,107 +0,0 @@
-import tempfile
-from pathlib import Path
-
-import pandas as pd
-
-from mDeepFRI import (MERGED_SEQUENCES, MMSEQS_SEARCH_RESULTS,
-                      TARGET_MMSEQS_DB_NAME)
-from mDeepFRI.utils.utils import merge_files_binary, run_command
-
-MMSEQS_COLUMN_NAMES = [
-    "query", "target", "identity", "alignment_length", "mismatches",
-    "gap_openings", "query_start", "query_end", "target_start", "target_end",
-    "e_value", "bit_score"
-]
-
-
-def createdb(sequences_file, db_path):
-    """
-    Converts FASTA file to a DB format needed for MMseqs2.
-    This should generate five files,
-    e.g. queryDB, queryDB_h and its corresponding index file queryDB.index,
-    queryDB_h.index and queryDB.lookup from the FASTA QUERY.fasta input sequences.
-
-    sequence_file (str): path to FASTA file.
-    db_path (str): path to output db file.
-
-    Returns:
-        None
-    """
-    run_command(f"mmseqs createdb {sequences_file} {db_path} --dbtype 1")
-
-
-def createindex(db_path):
-    with tempfile.TemporaryDirectory() as tmp_path:
-        run_command(f"mmseqs createindex {db_path} {tmp_path}")
-
-
-def search(query_db, target_db, result_db):
-    with tempfile.TemporaryDirectory() as tmp_path:
-        run_command(
-            f"mmseqs search {query_db} {target_db} {result_db} {tmp_path}")
-
-
-def convertalis(query_db, target_db, result_db, output_file):
-    run_command(
-        f"mmseqs convertalis {query_db} {target_db} {result_db} {output_file}")
-
-
-def create_target_database(seq_atoms_path: Path, new_db_path: Path) -> None:
-    """
-
-    :param seq_atoms_path:
-    :param new_db_path:
-    :param freshly_added_ids:
-    :return:
-    """
-    sequence_files = list((seq_atoms_path).glob("**/*.faa"))
-    print("\nMerging " + str(len(sequence_files)) +
-          " sequence files for mmseqs2")
-    merge_files_binary(sequence_files, seq_atoms_path / MERGED_SEQUENCES)
-
-    print("Creating new target mmseqs2 database " + str(new_db_path))
-    createdb(seq_atoms_path / MERGED_SEQUENCES,
-             new_db_path / TARGET_MMSEQS_DB_NAME)
-    print("Indexing new target mmseqs2 database " + str(new_db_path))
-    createindex(new_db_path / TARGET_MMSEQS_DB_NAME)
-
-
-def run_mmseqs_search(query_file: Path, target_db: Path, output_path: Path,
-                      min_bit_score: float, max_evalue: float,
-                      min_identity: float) -> pd.DataFrame:
-    """Creates a database from query sequences and runs mmseqs2 search against database.
-
-    Args:
-        query_file (pathlib.Path): Path to query FASTA file.
-        target_db (pathlib.Path): Path to target MMSeqs2 database.
-        output_path (pathlib.Path):
-
-    Returns:
-        pd.DataFrame: Pandas DataFrame with MMSeqs2 search results.
-    """
-    output_file = output_path / MMSEQS_SEARCH_RESULTS
-
-    query_db = output_path / 'queryDB'
-    createdb(query_file, query_db)
-
-    with tempfile.TemporaryDirectory() as tmp_path:
-
-        result_db = Path(tmp_path) / 'search_resultDB'
-        search(query_db, target_db, result_db)
-
-        # Convert results to tabular format
-        convertalis(query_db, target_db, result_db, output_file)
-
-    output = pd.read_csv(output_file, sep="\t", names=MMSEQS_COLUMN_NAMES)
-
-    # MMSeqs2 alginment filters
-    if min_identity:
-        output = output.query(f"identity >= {min_identity}")
-    if min_bit_score:
-        output = output.query(f"bitscore >= {min_bit_score}")
-    if max_evalue:
-        output = output.query(f"e_value <= {max_evalue}")
-
-    output.to_csv(output_file, sep="\t", index=False)
-
-    return output
diff --git a/mDeepFRI/utils/search_alignments.py b/mDeepFRI/utils/search_alignments.py
deleted file mode 100644
index 42187ba..0000000
--- a/mDeepFRI/utils/search_alignments.py
+++ /dev/null
@@ -1,141 +0,0 @@
-import json
-import logging
-import pathlib
-from functools import partial
-from multiprocessing.pool import ThreadPool
-
-import pandas as pd
-import parasail
-
-from mDeepFRI import ALIGNMENTS
-from mDeepFRI.utils.bio_utils import substitution_matrices
-from mDeepFRI.utils.fasta_file_io import SeqFileLoader
-
-
-def alignment_sequences_identity(query: str, target: str) -> float:
-    """
-    Calculate sequence identity between two proteins in an alignment.
-    Gaps are excluded.
-
-    Args:
-        query (str): Query sequence.
-        target (str): Target sequence.
-
-    Returns:
-        float: Sequence identity between two sequences in an alignment.
-    """
-    # take the shorter sequence
-    if len(query) > len(target):
-        query, target = target, query
-
-    # count matches excluding gaps
-    matches = 0
-    for query_char, target_char in zip(query, target):
-        if query_char == "-" or target_char == "-":
-            continue
-        elif query_char == target_char:
-            matches += 1
-
-    # calculate identity
-    seq_identity = matches / len(query)
-    return seq_identity
-
-
-# skipped in tests
-def align(query: str, target: str, matrix: parasail.bindings_v2.Matrix,
-          gap_open: int, gap_extend: int) -> parasail.bindings_v2.Result:
-    """
-    Align two protein sequences using biopython pairwise2.align.globalms.
-
-    Args:
-        query_seq (str): Query protein sequence.
-        target_seq (str): Target protein sequence.
-        matrix (parasail.bindings_v2.Matrix): Substitution matrix.
-        gap_open (int): Penalty score for opening a gap.
-        gap_continuation (int): Penalty score for gap continuation.
-
-    Returns:
-        parasail.bindings_v2.Result: Outputs Result object with required info.
-    """
-
-    alignment = parasail.nw_trace(query, target, gap_open, gap_extend, matrix)
-
-    return alignment
-
-
-def search_alignments(query_seqs: dict, mmseqs_search_output: pd.DataFrame,
-                      target_seqs: SeqFileLoader, output_path: pathlib.Path,
-                      matrix: str, alignment_gap_open: float,
-                      alignment_gap_extend: float,
-                      alignment_min_identity: float, threads: int):
-
-    # format of output JSON file:
-    # alignments = dict[query_id]
-    #     "target_id": target_id,
-    #     "sequence_identity" : alignment_sequence_identity(alignment),
-    #     "query_sequence": query protein sequence,
-    #     "target_sequence": target protein sequence,
-    #     "aln_score": alignment score
-
-    alignment_output_json_path = output_path / ALIGNMENTS
-
-    query_seqs_keys = list(query_seqs.keys())
-
-    filtered = mmseqs_search_output[mmseqs_search_output['query'].isin(
-        query_seqs_keys)]
-    logging.info("Filtered %i mmseqs matches",
-                 len(mmseqs_search_output) - len(filtered))
-    logging.info("Total alignments to check %i", len(filtered))
-
-    queries = list(map(lambda x: query_seqs[x], filtered["query"]))
-    targets = list(map(lambda x: target_seqs[x], filtered["target"]))
-
-    # parametrize alignment
-    parametrized_align = partial(align,
-                                 matrix=substitution_matrices[matrix],
-                                 gap_open=alignment_gap_open,
-                                 gap_extend=alignment_gap_extend)
-
-    # align with parasail
-    with ThreadPool(threads) as pool:
-        all_alignments = pool.starmap(parametrized_align,
-                                      zip(queries, targets))
-
-    alignments_output = dict()
-    for i, alignment in enumerate(all_alignments):
-        # filter out bad alignments based on alignment sequences identity
-        query_sequence = alignment.traceback.query
-        target_sequence = alignment.traceback.ref
-
-        sequence_identity = alignment_sequences_identity(
-            query_sequence, target_sequence)
-
-        if sequence_identity > alignment_min_identity:
-            query_id = filtered["query"].iloc[i]
-            target_id = filtered["target"].iloc[i]
-
-            # if query and it's alignment is not in the list yet - add it to list
-            if query_id not in alignments_output:
-                alignments_output[query_id] = {
-                    "target_id": target_id,
-                    "query_sequence": query_sequence,
-                    "target_sequence": target_sequence,
-                    "sequence_identity": sequence_identity,
-                    "aln_score": alignment.score
-                }
-
-            # select best alignment based on pairwise2.align.globalms.score
-            elif alignment.score > alignments_output[query_id]["aln_score"]:
-                alignments_output[query_id] = {
-                    "target_id": target_id,
-                    "query_sequence": query_sequence,
-                    "target_sequence": target_sequence,
-                    "sequence_identity": sequence_identity,
-                    "aln_score": alignment.score
-                }
-
-    json.dump(alignments_output,
-              open(alignment_output_json_path, "w", encoding="utf-8"),
-              indent=4,
-              sort_keys=True)
-    return alignments_output
diff --git a/pyproject.toml b/pyproject.toml
index 5e5ece5..5d46516 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,3 +1,7 @@
+[build-system]
+build-backend = "setuptools.build_meta"
+requires = ["setuptools", "numpy", "cython", "requests"]
+
 [tool.ruff]
 line-length = 120
 
diff --git a/setup.cfg b/setup.cfg
index 76b73e1..999a811 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,3 +1,3 @@
 [metadata]
-description-file=README.md
+description_file=README.md
 license_files=LICENSE
diff --git a/setup.py b/setup.py
index 747f6e4..af7087f 100644
--- a/setup.py
+++ b/setup.py
@@ -1,66 +1,154 @@
 import os
-import pathlib
-import sys
+import re
+import shutil
+import tarfile
+from distutils.util import convert_path
+from pathlib import Path
 
-from setuptools import Extension, find_namespace_packages, setup
-from setuptools.command.build_ext import build_ext as build_ext_orig
+import numpy as np
+import requests
+from setuptools import Extension, find_packages, setup
+from setuptools.command.build_ext import build_ext as _build_ext
 
-from mDeepFRI import __version__
+try:
+    from Cython.Build import cythonize
+except ImportError as err:
+    cythonize = err
 
+# --- Utils -----------------------------------------------------------------
 
-class CMakeExtension(Extension):
-    def __init__(self, name):
-        # don't invoke the original build_ext for this special extension
-        super().__init__(name, sources=[])
 
+def _detect_target_machine(platform):
+    if platform == "win32":
+        return "x86"
+    return platform.rsplit("-", 1)[-1]
+
+
+def _detect_target_cpu(platform):
+    machine = _detect_target_machine(platform)
+    if re.match("^mips", machine):
+        return "mips"
+    elif re.match("^(aarch64|arm64)$", machine):
+        return "aarch64"
+    elif re.match("^arm", machine):
+        return "arm"
+    elif re.match("(x86_64)|(x86)|(AMD64|amd64)|(^i.86$)", machine):
+        return "x86"
+    elif re.match("^(powerpc|ppc)", machine):
+        return "ppc"
+    return None
+
+
+def _detect_target_system(platform):
+    if platform.startswith("win"):
+        return "windows"
+    elif platform.startswith("macos"):
+        return "macos"
+    elif platform.startswith("linux"):
+        return "linux_or_android"
+    elif platform.startswith("freebsd"):
+        return "freebsd"
+    return None
+
+
+FOLDCOMP_BINARIES = {
+    "linux_or_android":
+    "https://mmseqs.com/foldcomp/foldcomp-linux-x86_64.tar.gz",
+    "aarch64": "https://mmseqs.com/foldcomp/foldcomp-linux-arm64.tar.gz",
+    "macos": "https://mmseqs.com/foldcomp/foldcomp-macos-universal.tar.gz",
+    "windows": "https://mmseqs.com/foldcomp/foldcomp-windows-x64.zip"
+}
+
+main_ns = {}
+ver_path = convert_path('mDeepFRI/__init__.py')
+with open(ver_path) as ver_file:
+    exec(ver_file.read(), main_ns)
+
+
+def read(fname):
+    return open(os.path.join(os.path.dirname(__file__), fname)).read()
+
+
+def _download_file(url, path):
+    with requests.get(url, stream=True) as r:
+        with open(path, 'wb') as f:
+            shutil.copyfileobj(r.raw, f)
+
+
+def _download_foldcomp(system, cpu, output_path):
+    # select appropriate binary
+    if cpu == "aarch64" and system == "linux":
+        build = "aarch64"
+    else:
+        build = system
+    url = FOLDCOMP_BINARIES[build]
+
+    output_path = Path(output_path)
+    foldcomp_tar = output_path / "foldcomp.tar.gz"
+
+    _download_file(url, foldcomp_tar)
+    # untar file
+    with tarfile.open(foldcomp_tar, "r:gz") as archive:
+        archive.extract("foldcomp", output_path)
+    # rename binary after extraction to foldcomp_bin
+    (output_path / "foldcomp").rename(output_path / "foldcomp_bin")
+
+    # remove tar file
+    foldcomp_tar.unlink()
+
+
+class build_ext(_build_ext):
+    def initialize_options(self) -> None:
+        _build_ext.initialize_options(self)
+        self.target_machine = None
+        self.target_cpu = None
+        self.target_system = None
 
-class build_ext(build_ext_orig):
     def run(self):
-        for ext in self.extensions:
-            self.build_cmake(ext)
-        super().run()
-
-    def build_cmake(self, ext):
-        cwd = pathlib.Path().absolute()
-
-        # these dirs will be created in build_py, so if you don't have
-        # any python sources to bundle, the dirs will be missing
-        build_temp = pathlib.Path(self.build_temp)
-        build_temp.mkdir(parents=True, exist_ok=True)
-        extdir = pathlib.Path(self.get_ext_fullpath(ext.name))
-        extdir.mkdir(parents=True, exist_ok=True)
-
-        # example of cmake args
-        python_version = f"{sys.version_info.major}.{sys.version_info.minor}"
-        python_include_path = str(
-            pathlib.Path(sys.executable).parent.parent.absolute()
-        ) + f'/include/python{python_version}/'
-
-        config = 'Debug' if self.debug else 'Release'
-        cmake_args = [
-            '-DCMAKE_LIBRARY_OUTPUT_DIRECTORY=' +
-            str(extdir.parent.absolute()),
-            '-DPY_INCLUDE_PATH=' + python_include_path,
-            '-DCMAKE_BUILD_TYPE=' + config,
-        ]
-
-        # example of build args
-        build_args = ['--config', config, '--', '-j4']
-
-        os.chdir(str(build_temp))
-        self.spawn(['cmake', str(cwd)] + cmake_args)
-        if not self.dry_run:
-            self.spawn(['cmake', '--build', '.'] + build_args)
-        # Troubleshooting: if fail on line above then delete all possible
-        # temporary CMake files including "CMakeCache.txt" in top level dir.
-        os.chdir(str(cwd))
+        if isinstance(cythonize, ImportError):
+            raise RuntimeError("Failed to import Cython") from cythonize
 
+        self.extensions = cythonize(self.extensions,
+                                    compiler_directives={
+                                        "linetrace": True,
+                                        "language_level": 3
+                                    })
+        for ext in self.extensions:  # this fixes a bug with setuptools
+            ext._needs_stub = False
+        _build_ext.run(self)
+        self.target_machine = _detect_target_machine(self.plat_name)
+        self.target_cpu = _detect_target_cpu(self.plat_name)
+        self.target_system = _detect_target_system(self.plat_name)
+        _download_foldcomp(self.target_system, self.target_cpu, self.build_lib)
+
+
+SRC_DIR = "mDeepFRI"
+PACKAGES = [SRC_DIR]
+
+install_requires = ["cython", "numpy"]
+setup_requires = ["cython"]
+
+EXTENSIONS = [
+    Extension("mDeepFRI.predict",
+              sources=[SRC_DIR + "/predict.pyx"],
+              language="c++",
+              libraries=["stdc++"],
+              extra_compile_args=["-std=c++17", "-O3"]),
+    Extension("mDeepFRI.alignment_utils",
+              sources=[SRC_DIR + "/alignment_utils.pyx"],
+              language="c++",
+              libraries=["stdc++"],
+              extra_compile_args=["-std=c++17", "-O3"]),
+]
 
 setup(
     name="mDeepFRI",
-    version=__version__,
+    version=main_ns['__version__'],
     description=
     "Pipeline for searching and aligning contact maps for proteins, then running DeepFri's GCN.",
+    long_description=read("README.md"),
+    long_description_content_type='text/markdown',
+    keywords="protein function metagenomics deep neural network",
     author="Piotr Kucharski, Valentyn Bezshapkin",
     author_email=
     "soliareofastorauj@gmail.com, valentyn.bezshapkin@micro.biol.ethz.ch",
@@ -72,12 +160,12 @@ def build_cmake(self, ext):
             "mDeepFRI = mDeepFRI.cli:main",
         ],
     },
-    ext_modules=[CMakeExtension("mDeepFRI/CPP_lib/library_definition")],
-    cmdclass={
-        'build_ext': build_ext,
-    },
+    ext_modules=EXTENSIONS,
+    include_dirs=[np.get_include()],
+    install_requires=install_requires,
+    cmdclass={'build_ext': build_ext},
     license="GNU GPLv3",
-    packages=find_namespace_packages(),
+    packages=find_packages(),
     classifiers=[
         'Programming Language :: Python',
         'Programming Language :: Python :: 3',


Test patch:
diff --git a/tests/data/1S3P-A.bin.cmap.npy b/tests/data/1S3P-A.bin.cmap.npy
new file mode 100644
index 0000000..95b4538
Binary files /dev/null and b/tests/data/1S3P-A.bin.cmap.npy differ
diff --git a/tests/data/1S3P-A.bin.cmap.pkl b/tests/data/1S3P-A.bin.cmap.pkl
deleted file mode 100644
index 8bb98a7..0000000
Binary files a/tests/data/1S3P-A.bin.cmap.pkl and /dev/null differ
diff --git a/tests/data/small_query.faa b/tests/data/small_query.faa
index f6d536b..e35a897 100644
--- a/tests/data/small_query.faa
+++ b/tests/data/small_query.faa
@@ -15,4 +15,6 @@ AKHLFTSESVSEGHPDKIADQISDAVLDAILEQDPKARVACETYVKTGMVLVGGEITTSAWVDIEEITRNTVREIGYVHS
 >1YT3_1|Chain A|Ribonuclease D|Escherichia coli (562)
 MNYQMITTDDALASLCEAVRAFPAIALDTEFVRTRTYYPQLGLIQLFDGEHLALIDPLGITDWSPLKAILRDPSITKFLHAGSEDLEVFLNVFGELPQPLIDTQILAAFCGRPMSWGFASMVEEYSGVTLDKSESRTDWLARPLTERQCEYAAADVWYLLPITAKLMVETEASGWLPAALDECRLMQMRRQEVVAPEDAWRDITNAWQLRTRQLACLQLLADWRLRKARERDLAVNFVVREEHLWSVARYMPGSLGELDSLGLSGSEIRFHGKTLLALVEKAQTLPEDALPQPMLNLMDMPGYRKAFKAIKSLITDVSETHKISAELLASRRQINQLLNWHWKLKPQNNLPELISGWRGELMAEALHNLLQEYPQ
 >1NZO_1_mutant|Chain A|Penicillin-binding protein 5|Escherichia coli (562)
-DDLNIKTMIPGVPQIDAESYILIDYNSGKVLAEQNADVRRDPASLTKMMTSYVIGQAMKAGKFKETDLVTIGNDAWATGNPVFKGSSLLQSKPGMQINLQSGNVPVSQLIRGDACVAMADFAAGSQDAFVGLMNSYVNALGLKNTHFQTVHGLDADGQYSSARDMALIGQALIRDVPNEYSIYKEKEFTFNGIRQLNRNGLLWDNSLNVDGIKTGHTDKAGYNLVASATEGQMRLISAVMGGRTFKGREAESKKLLTWGFRFFETVNPLKVGKEFASEPVWFGDSDRASLGVDKDVYLTIPRGRMKDLKASYVLNSSELHAPLQKNQVVGTINFQLDGKTIEQRPLVVLQEIPEGNFGDPVID
\ No newline at end of file
+DDLNIKTMIPGVPQIDAESYILIDYNSGKVLAEQNADVRRDPASLTKMMTSYVIGQAMKAGKFKETDLVTIGNDAWATGNPVFKGSSLLQSKPGMQINLQSGNVPVSQLIRGDACVAMADFAAGSQDAFVGLMNSYVNALGLKNTHFQTVHGLDADGQYSSARDMALIGQALIRDVPNEYSIYKEKEFTFNGIRQLNRNGLLWDNSLNVDGIKTGHTDKAGYNLVASATEGQMRLISAVMGGRTFKGREAESKKLLTWGFRFFETVNPLKVGKEFASEPVWFGDSDRASLGVDKDVYLTIPRGRMKDLKASYVLNSSELHAPLQKNQVVGTINFQLDGKTIEQRPLVVLQEIPEGNFGDPVID
+>freak
+RPLVVLQEIPEGGVDKDVYLGDSDRTIPRGRMKDLKASTTTTTTTNVPVSQLIKKKWDFDFKNTDKAGYNLSSLLPVFKGMQQSKPGHFQTVHGLDADGQYSSARDDFDFTTTTTNEYSIYKDFTSYVNALGLEKWWWWWWWWWWWRNGLLWDNSLNVDGIKTGHTTIRDVPTVASATEGQMRLISAVMGGRTFKGREAELHAPLQKNQVVGTINFQLDGKTIEQYVLNSS
\ No newline at end of file
diff --git a/tests/test_alignment.py b/tests/test_alignment.py
index a194fa4..ec7697e 100644
--- a/tests/test_alignment.py
+++ b/tests/test_alignment.py
@@ -1,11 +1,26 @@
 import pytest
 
-from mDeepFRI.utils.search_alignments import alignment_sequences_identity
+from mDeepFRI.alignment import insert_gaps
+from mDeepFRI.alignment_utils import alignment_identity
 
 
 @pytest.mark.parametrize('query,target,expected', [('AASDS', 'AASDS', 1),
                                                    ('AASDS', 'ADSDS', 0.8),
                                                    ('AASSS', 'A-S-S', 0.6),
                                                    ('A--SS', 'A-S--', 0.2)])
-def test_alignment_sequences_identity(query, target, expected):
-    assert alignment_sequences_identity(query, target) == expected
+def test_alignment_identity(query, target, expected):
+    seq_id = alignment_identity(query, target)
+    assert round(seq_id, 6) == expected
+
+
+@pytest.mark.parametrize(
+    'query,reference,alignment,gapped_query,gapped_reference',
+    [('AACT', 'AAT', 'MMDM', 'AACT', 'AA-T'),
+     ('AAT', 'AATC', 'MMMI', 'AAT-', 'AATC'),
+     ('AAT', 'FGTC', 'XXMI', 'AAT-', 'FGTC')])
+def test_insert_gaps_full_length(query, reference, alignment, gapped_query,
+                                 gapped_reference):
+    gapped_seq, gapped_ref = insert_gaps(query, reference, alignment)
+    assert gapped_seq == gapped_query
+    assert gapped_ref == gapped_reference
+    assert len(gapped_seq) == len(gapped_ref)
diff --git a/tests/test_cli.py b/tests/test_cli.py
index 2db87e3..84d1e94 100644
--- a/tests/test_cli.py
+++ b/tests/test_cli.py
@@ -20,7 +20,6 @@ def test_cli_build_db(runner):
     result = runner.invoke(cli.main, ["build-db", "--help"])
     assert result.exit_code == 0
     assert "--help" in result.output
-    assert "-m" in result.output
     assert "-t" in result.output
 
 
diff --git a/tests/test_cpp_lib.py b/tests/test_cpp_lib.py
deleted file mode 100644
index 240f877..0000000
--- a/tests/test_cpp_lib.py
+++ /dev/null
@@ -1,66 +0,0 @@
-import filecmp
-import pickle
-import tempfile
-from typing import Callable
-
-import numpy as np
-import pytest
-
-from mDeepFRI.CPP_lib.libAtomDistanceIO import (initialize,
-                                                load_aligned_contact_map,
-                                                load_contact_map, save_atoms)
-from mDeepFRI.structure_files.parsers import parse_pdb
-
-
-@pytest.fixture()
-def reference_binary():
-    return "tests/data/1S3P-A.bin"
-
-
-@pytest.fixture()
-def expected_contact_map():
-    with open("tests/data/1S3P-A.bin.cmap.pkl", "rb") as f:
-        expected_contact_map = pickle.load(f)
-    return expected_contact_map
-
-
-def error_wrapper(func: Callable, *args, error_contains: str = None):
-    try:
-        func(*args)
-    except Exception as e:
-        if error_contains in str(e):
-            return
-        else:
-            raise e
-
-
-def test_functions():
-    initialize()
-
-    error_wrapper(save_atoms, error_contains="did not match C++ signature")
-    error_wrapper(load_contact_map,
-                  error_contains="did not match C++ signature")
-    error_wrapper(load_aligned_contact_map,
-                  error_contains="did not match C++ signature")
-
-
-def test_load_contact_map(reference_binary, expected_contact_map):
-    initialize()
-    contact_map_cpp = load_contact_map(reference_binary, 6)
-    assert np.array_equal(contact_map_cpp, expected_contact_map)
-
-
-def test_save_atoms(expected_contact_map):
-    with open("tests/data/structures/1S3P-A.pdb") as f:
-        _, positions, groups = parse_pdb(f)
-    _, group_index = np.unique(groups, return_index=True)
-    group_index = np.sort(group_index)
-    group_indexes = np.append(group_index, positions.shape[0]).astype(np.int32)
-    initialize()
-    with tempfile.TemporaryDirectory() as tmpdirname:
-        save_atoms(positions, group_indexes, tmpdirname + "/1S3P-A.bin")
-        assert filecmp.cmp(tmpdirname + "/1S3P-A.bin", "tests/data/1S3P-A.bin")
-
-
-def test_load_aligned_contact_map():
-    pass
diff --git a/tests/test_install.py b/tests/test_install.py
new file mode 100644
index 0000000..9ebf6fc
--- /dev/null
+++ b/tests/test_install.py
@@ -0,0 +1,18 @@
+import subprocess
+from pathlib import Path
+
+import mDeepFRI
+
+
+# is a tricky test, will fail if ran from a wrong directory
+def test_foldcomp_download():
+    """
+    Test if foldcomp binary was downloaded properly
+    and can be invoked from Python.
+    """
+    foldcomp_bin = Path(mDeepFRI.__path__[0]).parent / "foldcomp"
+    stdout = subprocess.run([foldcomp_bin],
+                            capture_output=True,
+                            check=True,
+                            text=True).stdout
+    assert "foldcomp compress" in stdout
diff --git a/tests/test_parsers.py b/tests/test_parsers.py
index 94d170a..a6d13f7 100644
--- a/tests/test_parsers.py
+++ b/tests/test_parsers.py
@@ -4,7 +4,7 @@
 import numpy as np
 import pytest
 
-from mDeepFRI.structure_files.parsers import parse_mmcif, parse_pdb
+from mDeepFRI.parsers import parse_mmcif, parse_pdb
 
 
 @pytest.fixture
@@ -12,56 +12,56 @@ def pdb_file():
     test_file = pathlib.Path(
         __file__
     ).parent / "data" / "structures" / "AF-A0A2Z5TJB0-F1-model_v4.pdb.gz"
-    return gzip.open(test_file, "rt")
+    return gzip.open(test_file, "rb").read().split(b"\n")
 
 
 def test_parse_pdb_sequence(pdb_file):
     sequence, _, _ = parse_pdb(pdb_file)
 
-    assert (sequence[:4] == ['MET', 'MET', 'MET', 'MET']).all()
+    assert sequence[:4] == [b'MET', b'MET', b'MET', b'MET']
 
 
 def test_parse_pdb_positions(pdb_file):
     _, positions, _ = parse_pdb(pdb_file)
-    expected = np.array([[-47.025, 49.762, -23.074]], dtype=np.float32)
-    assert np.array_equal(positions[:1], expected)
+    expected = np.array([-47.025, 49.762, -23.074], dtype=np.float32)
+    assert np.array_equal(positions[:3], expected)
 
 
 def test_parse_pdb_groups(pdb_file):
     _, _, groups = parse_pdb(pdb_file)
-    assert groups[:1] == ['A   1']
+    assert groups[:1] == [b'A   1']
 
 
 def test_parse_pdb_lengths(pdb_file):
     sequence, positions, groups = parse_pdb(pdb_file)
     assert len(sequence) > 0
-    assert len(sequence) == len(positions) == len(groups)
+    assert len(sequence) == len(positions) / 3 == len(groups)
 
 
 @pytest.fixture
 def mmcif_file():
     test_file = pathlib.Path(
         __file__).parent / "data" / "structures" / "6a0j.cif.gz"
-    return gzip.open(test_file, "rt")
+    return gzip.open(test_file, "rb").read().split(b"\n")
 
 
 def test_parse_mmcif_sequence(mmcif_file):
     sequence, _, _ = parse_mmcif(mmcif_file)
-    assert all(sequence[:4] == ['VAL', 'VAL', 'VAL', 'VAL'])
+    assert sequence[:4] == [b'VAL', b'VAL', b'VAL', b'VAL']
 
 
 def test_parse_mmcif_positions(mmcif_file):
     _, positions, _ = parse_mmcif(mmcif_file)
-    expected = np.array([[-6.123, 31.928, 92.243]], dtype=np.float32)
-    assert np.array_equal(positions[:1], expected)
+    expected = np.array([-6.123, 31.928, 92.243], dtype=np.float32)
+    assert np.array_equal(positions[:3], expected)
 
 
 def test_parse_mmcif_groups(mmcif_file):
     _, _, groups = parse_mmcif(mmcif_file)
-    assert groups[:1] == ['A22']
+    assert groups[:1] == [b'A22']
 
 
 def test_parse_mmcif_lengths(mmcif_file):
     sequence, positions, groups = parse_mmcif(mmcif_file)
     assert len(sequence) > 0
-    assert len(sequence) == len(positions) == len(groups)
+    assert len(sequence) == len(positions) / 3 == len(groups)
diff --git a/tests/test_structure_parsers.py b/tests/test_structure_parsers.py
deleted file mode 100644
index f645bb0..0000000
--- a/tests/test_structure_parsers.py
+++ /dev/null
@@ -1,29 +0,0 @@
-import pathlib
-
-import numpy as np
-import pytest
-
-from mDeepFRI.structure_files.parse_structure_file import (SeqAtoms,
-                                                           read_structure_file)
-
-
-@pytest.mark.parametrize(
-    "file_name, excepted_id, excepted_amino_group, excepted_position, excepted_group, excepted_length",
-    [("6a0j.cif.gz", "6a0j", ['VAL', 'VAL', 'VAL', 'VAL'
-                              ], [-6.123, 31.928, 92.243], ['A22'], 3443),
-     ("AF-A0A2Z5TJB0-F1-model_v4.pdb.gz", "AF-A0A2Z5TJB0-F1-model_v4", [
-         'MET', 'MET', 'MET', 'MET'
-     ], [-47.025, 49.762, -23.074], ['A   1'], 3956)])
-def test_proper_structure_files(file_name, excepted_id, excepted_amino_group,
-                                excepted_position, excepted_group,
-                                excepted_length):
-    file_path = pathlib.Path("tests") / "data" / "structures" / file_name
-    output: SeqAtoms = read_structure_file(file_path)
-
-    assert output.protein_id == excepted_id
-    assert (output.atom_amino_group[:4] == excepted_amino_group).all()
-    assert (output.positions[:1] == np.array([excepted_position],
-                                             dtype=np.float32)).all()
-    assert output.groups[:1] == excepted_group
-    assert len(output.atom_amino_group) == len(output.positions) == len(
-        output.groups) == excepted_length
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 07a2ac8..b56f6e7 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -1,6 +1,4 @@
-import pytest
-
-from mDeepFRI.utils import bio_utils, utils
+from mDeepFRI import bio_utils, utils
 
 
 def test_protein_letters():
@@ -36,25 +34,6 @@ def test_protein_letters():
     assert bio_utils.PROTEIN_LETTERS == expected
 
 
-def test_load_deepfri_config(tmp_path):
-    """Test load_deepfri_config() function."""
-    tmp_json_file = tmp_path / "model_config.json"
-    tmp_json_file.write_text(
-        '{"model1": "./trained_models/model1.h5", "model2": "./trained_models/model2.h5"}'
-    )
-
-    config = utils.load_deepfri_config(str(tmp_json_file))
-
-    assert config["model1"] == str(tmp_path / "model1.h5")
-    assert config["model2"] == str(tmp_path / "model2.h5")
-
-
-def test_load_deepfri_config_file_not_found():
-    """Test FileNotFoundError is raised if config file is not found."""
-    with pytest.raises(FileNotFoundError):
-        utils.load_deepfri_config("non_existing_file.json")
-
-
 def test_run_command():
     command = "echo 'Hello World!'"
     result = utils.run_command(command).strip()