Repository: WrightonLabCSU/DRAM (language: python)

Issue:
Circleci project setup
Add CircleCI support to merge into master


Discussion / hints:

Fix patch:
diff --git a/.circleci/config.yml b/.circleci/config.yml
new file mode 100644
index 00000000..74eeaa31
--- /dev/null
+++ b/.circleci/config.yml
@@ -0,0 +1,44 @@
+# Use the latest 2.1 version of CircleCI pipeline process engine.
+# See: https://circleci.com/docs/2.0/configuration-reference
+version: 2.1
+
+# TODO More Circleci features may be usfull
+
+# Define jobs
+jobs:
+  build-and-test:
+    docker:
+      - image: ubuntu:focal
+    steps:
+      - checkout
+      - run:
+          name: Setup Miniconda
+          command: |
+            apt update
+            apt install -y wget
+            cd $HOME
+            wget -q https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh
+            export MINICONDA_PREFIX="$HOME/miniconda"
+            bash miniconda.sh -b -p $MINICONDA_PREFIX
+            export PATH="$MINICONDA_PREFIX/bin:$PATH"
+            conda config --set always_yes yes
+            conda update -q conda
+            conda config --add channels bioconda
+            conda config --add channels conda-forge
+            conda info -a
+      - run:
+          name: Run tests in enviroment
+          # This assumes pytest is installed via the install-package step above
+          command: |
+            export PATH="$HOME/miniconda/bin:$PATH"
+            conda update -y conda
+            conda create -n myenv python=$PYTHON_VERSION
+            wget https://raw.githubusercontent.com/shafferm/DRAM/master/environment.yaml
+            conda env create -f environment.yaml -n DRAM
+            source activate DRAM
+            pytest
+# Invoke jobs via workflows
+workflows:
+  all-tests: 
+    jobs:
+      - build-and-test
diff --git a/.travis.yml b/.travis.yml
deleted file mode 100644
index 8efb0512..00000000
--- a/.travis.yml
+++ /dev/null
@@ -1,26 +0,0 @@
-language: python
-python:
-- '3.6'
-- '3.7'
-cache: pip
-before_install:
-- wget -q https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O
-  miniconda.sh
-- export MINICONDA_PREFIX="$HOME/miniconda"
-- bash miniconda.sh -b -p $MINICONDA_PREFIX
-- export PATH="$MINICONDA_PREFIX/bin:$PATH"
-- conda config --set always_yes yes
-- conda update -q conda
-- conda config --add channels bioconda
-- conda config --add channels conda-forge
-- conda info -a
-install:
-- conda create -n test-env pytest pytest-cov python-coveralls python=$TRAVIS_PYTHON_VERSION pandas scikit-bio prodigal
-  mmseqs2!=10.6d92c hmmer trnascan-se >=2 sqlalchemy barrnap altair >=4 openpyxl networkx ruby parallel
-- source activate test-env
-- conda list
-- pip install --no-cache-dir --editable .
-script:
-- py.test --cov DRAM --cov-report term-missing
-after_success:
-- coveralls
diff --git a/README.md b/README.md
index 64de80cd..f3510326 100644
--- a/README.md
+++ b/README.md
@@ -1,5 +1,5 @@
 # DRAM
-[![Build Status](https://travis-ci.com/shafferm/checkMetab.svg?branch=master)](https://travis-ci.com/shafferm/checkMetab)
+[![CircleCI](https://circleci.com/gh/WrightonLabCSU/DRAM/tree/master.svg?style=svg)](https://circleci.com/gh/WrightonLabCSU/DRAM/tree/master)
 
 DRAM (Distilled and Refined Annotation of Metabolism) is a tool for annotating metagenomic assembled genomes and [VirSorter](https://github.com/simroux/VirSorter) identified viral contigs. DRAM annotates MAGs and viral contigs using [KEGG](https://www.kegg.jp/) (if provided by the user), [UniRef90](https://www.uniprot.org/), [PFAM](https://pfam.xfam.org/), [dbCAN](http://bcb.unl.edu/dbCAN2/), [RefSeq viral](https://www.ncbi.nlm.nih.gov/genome/viruses/), [VOGDB](http://vogdb.org/) and the [MEROPS](https://www.ebi.ac.uk/merops/) peptidase database as well as custom user databases. DRAM is run in two stages. First an annotation step to assign database identifiers to gene and then a distill step to curate these annotations into useful functional categories. Additionally viral contigs are further analyzed during to identify potential AMGs. This is done via assigning an auxiliary score and flags representing the confidence that a gene is both metabolic and viral.
 
diff --git a/mag_annotator/CONFIG b/mag_annotator/CONFIG
index 02364a98..ffd7360f 100644
--- a/mag_annotator/CONFIG
+++ b/mag_annotator/CONFIG
@@ -1,1 +1,1 @@
-{"description_db": "None", "kegg": null, "kofam": null, "kofam_ko_list": null, "uniref": null, "pfam": null, "pfam_hmm_dat": null, "dbcan": null, "dbcan_fam_activities": null, "viral": null, "peptidase": null, "vogdb": null, "vog_annotations": null, "genome_summary_form": "/Users/shafferm/lab/DRAM/data/genome_summary_form.tsv", "module_step_form": "/Users/shafferm/lab/DRAM/data/module_step_form.tsv", "etc_module_database": "/Users/shafferm/lab/DRAM/data/etc_module_database.tsv", "function_heatmap_form": "/Users/shafferm/lab/DRAM/data/function_heatmap_form.tsv", "amg_database": "/Users/shafferm/lab/DRAM/data/amg_database.tsv"}
\ No newline at end of file
+{"description_db": "None", "kegg": null, "kofam": null, "kofam_ko_list": null, "uniref": null, "pfam": null, "pfam_hmm_dat": null, "dbcan": null, "dbcan_fam_activities": null, "viral": null, "peptidase": null, "vogdb": null, "vog_annotations": null, "genome_summary_form": null, "module_step_form": null, "etc_module_database": null, "function_heatmap_form": null, "amg_database": null}
\ No newline at end of file
diff --git a/mag_annotator/annotate_bins.py b/mag_annotator/annotate_bins.py
index 16966c78..c746d206 100644
--- a/mag_annotator/annotate_bins.py
+++ b/mag_annotator/annotate_bins.py
@@ -1,3 +1,11 @@
+import re
+import io
+import time
+import warnings
+from glob import glob
+from functools import partial
+from datetime import datetime
+from typing import Callable
 from skbio.io import read as read_sequence
 from skbio.io import write as write_sequence
 from skbio import Sequence
@@ -5,15 +13,9 @@
 from os import path, mkdir, stat
 from shutil import rmtree, copy2
 import pandas as pd
-from datetime import datetime
-import re
-from glob import glob
-import warnings
-from functools import partial
-import io
-import time
 
-from mag_annotator.utils import run_process, make_mmseqs_db, merge_files, get_database_locs, multigrep, remove_suffix
+
+from mag_annotator.utils import run_process, make_mmseqs_db, merge_files, multigrep, remove_suffix
 from mag_annotator.database_handler import DatabaseHandler
 
 # TODO: add ability to take into account multiple best hits as in old_code.py
@@ -23,6 +25,7 @@
 # TODO: in annotated gene faa checkout out ko id for actual kegg gene id
 # TODO: add ability to handle [] in file names
 
+MAG_DBS_TO_ANNOTATE = ('kegg', 'kofam', 'kofam_ko_list', 'uniref', 'peptidase', 'pfam', 'dbcan', 'vogdb')
 BOUTFMT6_COLUMNS = ['qId', 'tId', 'seqIdentity', 'alnLen', 'mismatchCnt', 'gapOpenCnt', 'qStart', 'qEnd', 'tStart',
                     'tEnd', 'eVal', 'bitScore']
 HMMSCAN_ALL_COLUMNS = ['query_id', 'query_ascession', 'query_length', 'target_id', 'target_ascession', 'target_length',
@@ -31,7 +34,6 @@
                        'alignment_end', 'query_start', 'query_end', 'accuracy', 'description']
 HMMSCAN_COLUMN_TYPES = [str, str, int, str, str, int, float, float, float, int, int, float, float, float, float, int,
                         int, int, int, int, int, float, str]
-MAG_DBS_TO_ANNOTATE = ('kegg', 'kofam', 'kofam_ko_list', 'uniref', 'peptidase', 'pfam', 'dbcan', 'vogdb')
 
 
 def filter_fasta(fasta_loc, min_len=5000, output_loc=None):
@@ -110,11 +112,11 @@ def process_reciprocal_best_hits(forward_output_loc, reverse_output_loc, target_
     forward_hits = forward_hits.set_index('qId')
     reverse_hits = pd.read_csv(reverse_output_loc, sep='\t', header=None, names=BOUTFMT6_COLUMNS)
     reverse_hits = reverse_hits.set_index('qId')
-    def check_hit(row:pd.Series):
+
+    def check_hit(row: pd.Series):
         rbh = False
         if row.tId in reverse_hits.index:
-            if row.name == reverse_hits.loc[row.tId].tId:
-                rbh = True
+            rbh = row.name == reverse_hits.loc[row.tId].tId
         return {'%s_hit' % target_prefix:      row.tId,
                 '%s_RBH' % target_prefix:      rbh,
                 '%s_identity' % target_prefix: row.seqIdentity,
@@ -122,7 +124,7 @@ def check_hit(row:pd.Series):
                 '%s_eVal' % target_prefix:     row.eVal,
                 'index':                       row.name
                 }
-    hits =  forward_hits.apply(check_hit, axis=1, result_type='expand')
+    hits = forward_hits.apply(check_hit, axis=1, result_type='expand')
     # NOTE these lines may not be necessary
     hits.set_index('index', drop=True, inplace=True)
     hits.index.name = None
@@ -218,8 +220,9 @@ def run_mmseqs_profile_search(query_db, pfam_profile, output_loc, output_prefix=
         return pd.Series(name='%s_hits' % output_prefix)
 
 
-def get_sig(tstart, tend, tlen, evalue):
+def get_sig_row(row):
     """Check if hmm match is significant, based on dbCAN described parameters"""
+    tstart, tend, tlen, evalue = row[['target_start', 'target_end', 'target_length', 'full_evalue']].values
     perc_cov = (tend - tstart)/tlen
     if perc_cov >= .35 and evalue <= 1e-15:
         return True
@@ -241,140 +244,128 @@ def parse_hmmsearch_domtblout(file):
     return hmmsearch_frame
 
 
-def run_hmmscan_kofam(gene_faa, kofam_hmm, output_dir, ko_list, top_hit=True, use_dbcan2_thresholds=False, threads=1,
-                      verbose=False):
-    output = path.join(output_dir, 'kofam_profile.b6')
-    run_process(['hmmsearch', '--domtblout', output, '--cpu', str(threads), kofam_hmm, gene_faa], verbose=verbose)
-    if path.isfile(output) and stat(output).st_size > 0:
-        ko_hits = parse_hmmsearch_domtblout(output)
-
-        # check what method for determining what is significant should be used
-        if use_dbcan2_thresholds:
-            ko_hits_sig = ko_hits.loc[[get_sig(row.target_start, row.target_end,
-                                               row.target_length, row.full_evalue) for _, row in ko_hits.iterrows()]]
-        else:
-            is_sig = list()
-            for ko, frame in ko_hits.groupby('target_id'):
-                ko_row = ko_list.loc[ko]
-                if ko_row['score_type'] == 'domain':
-                    score = frame.domain_score
-                elif ko_row['score_type'] == 'full':
-                    score = frame.full_score
-                elif ko_row['score_type'] == '-':
-                    continue
-                else:
-                    raise ValueError(ko_row['score_type'])
-                frame = frame.loc[score.astype(float) > float(ko_row.threshold)]
-                is_sig.append(frame)
-            if len(is_sig) > 0:
-                ko_hits_sig = pd.concat(is_sig)
-            else:
-                ko_hits_sig = []
-        # if there are any significant results then parse to dataframe
-        if len(ko_hits_sig) == 0:
-            return pd.DataFrame(columns=['ko_id', 'kegg_hit'])
+# continue exit()
+def sig_scores(hits:pd.DataFrame, score_db:pd.DataFrame) -> pd.DataFrame:
+    is_sig = list()
+    for i, frame in hits.groupby('target_id'):
+        row = score_db.loc[i]
+        if row['score_type'] == 'domain':
+            score = frame.domain_score
+        elif row['score_type'] == 'full':
+            score = frame.full_score
+        elif row['score_type'] == '-':
+            continue
         else:
-            kegg_dict = dict()
-            for gene, frame in ko_hits_sig.groupby('query_id'):
-                # TODO: take top hit for full length genes and all hits for domains?
-                # TODO: if top hit then give all e-value and bitscore info
-                if top_hit:
-                    best_hit = frame[frame.full_evalue == frame.full_evalue.min()]
-                    ko_id = best_hit['target_id'].iloc[0]
-                    kegg_dict[gene] = [ko_id, ko_list.loc[ko_id, 'definition']]
-                else:
-                    kegg_dict[gene] = [','.join([i for i in frame.target_id]),
-                                       '; '.join([ko_list.loc[i, 'definition'] for i in frame.target_id])]
-            return pd.DataFrame(kegg_dict, index=['ko_id', 'kegg_hit']).transpose()
+            raise ValueError(row['score_type'])
+        frame = frame.loc[score.astype(float) > float(row.threshold)]
+        is_sig.append(frame)
+    if len(is_sig) > 0:
+        return pd.concat(is_sig)
     else:
-        return pd.DataFrame(columns=['ko_id', 'kegg_hit'])
-
-
-# refactor hmmscan_dbcan and hmmscan_vogdb to merge
-def run_hmmscan_dbcan(genes_faa, dbcan_loc, output_loc, threads=10, db_handler=None, verbose=False):
-    """Run hmmscan of genes against dbcan, apparently I can speed it up using hmmsearch in the reverse
-    Commands this is based on:
-    hmmscan --domtblout ~/dbCAN_test_1 dbCAN-HMMdb-V7.txt ~/shale_checkMetab_test/DRAM/genes.faa
-    cat ~/dbCAN_test_1 | grep -v '^#' | awk '{print $1,$3,$4,$6,$13,$16,$17,$18,$19}' | sed 's/ /\t/g' | \
-    sort -k 3,3 -k 8n -k 9n > dbCAN_test_1.good_cols.tsv
-    """
-    # Run hmmscan
-    dbcan_output = path.join(output_loc, 'dbcan_results.unprocessed.txt')
-    run_process(['hmmsearch', '--domtblout', dbcan_output, '--cpu', str(threads), dbcan_loc, genes_faa],
-                verbose=verbose)
-
-    # Process results
-    if path.isfile(dbcan_output) and stat(dbcan_output).st_size > 0:
-        dbcan_res = parse_hmmsearch_domtblout(dbcan_output)
+        return pd.DataFrame()
 
-        significant = [row_num for row_num, row in dbcan_res.iterrows() if get_sig(row.target_start, row.target_end,
-                                                                                   row.target_length, row.full_evalue)]
-        if len(significant) == 0:  # if nothing significant then return nothing, don't get descriptions
-            return pd.Series(name='cazy_hits')
-        dbcan_res_significant = dbcan_res.loc[significant]
 
-        dbcan_dict = dict()
-        if db_handler is not None:
-            dbcan_descriptions = db_handler.get_descriptions(set([strip_endings(i, ['.hmm']).split('_')[0] for i in
-                                                                  dbcan_res_significant.target_id]),
-                                                             'dbcan_description')
-        else:
-            dbcan_descriptions = None
-        for gene, frame in dbcan_res_significant.groupby('query_id'):
-            if dbcan_descriptions is None:
-                dbcan_dict[gene] = '; '.join([i[:-4] for i in frame.target_id])
-            else:
-                dbcan_dict[gene] = '; '.join(['%s [%s]' % (dbcan_descriptions.get(accession[:-4].split('_')[0]),
-                                                           accession[:-4]) for accession in frame.target_id])
-        return pd.Series(dbcan_dict, name='cazy_hits')
+def dbcan_hmmscan_formater(hits:pd.DataFrame,  db_name:str, db_handler=None):
+    hits_sig = hits[hits.apply(get_sig_row, axis=1)]
+    if len(hits_sig) == 0:
+        # if nothing significant then return nothing, don't get descriptions
+        return pd.Series()
+    hits_df = hits_sig.groupby('query_id').apply(
+        lambda x: '; '.join(x['target_id'].apply(lambda y:y[:-4]).unique())
+    )
+    hits_df = pd.DataFrame(hits_df)
+    hits_df.columns = [f"{db_name}_id"]
+    if db_handler is not None:
+        hits_df[f"{db_name}_hits"] = hits_df[f"{db_name}_id"].apply(
+            lambda x:'; '.join(db_handler.get_descriptions(x.split('; '), 'dbcan_description').values()))
+    hits_df.rename_axis(None, inplace=True)
+    return hits_df
+
+
+#TODO decide if we need use_hmmer_thresholds:bool=False
+def generic_hmmscan_formater(hits:pd.DataFrame,  db_name:str, hmm_info_path:str=None, top_hit:bool=True):
+    if hmm_info_path is None:
+        hmm_info = None
+        hits_sig = hits[hits.apply(get_sig_row, axis=1)]
     else:
-        return pd.Series(name='cazy_hits')
+        hmm_info = pd.read_csv(hmm_info_path, sep='\t', index_col=0)
+        hits_sig = sig_scores(hits, hmm_info)
+    if len(hits_sig) == 0:
+        # if nothing significant then return nothing, don't get descriptions
+        return pd.DataFrame()
+    if top_hit:
+        # Get the best hits
+        hits_sig = hits_sig.sort_values('full_evalue').drop_duplicates(subset=["query_id"])
+    hits_df = hits_sig[['target_id', 'query_id']]
+    hits_df.set_index('query_id', inplace=True, drop=True)
+    hits_df.rename_axis(None, inplace=True)
+    hits_df.columns = [f"{db_name}_id"]
+    if hmm_info is not None:
+        hits_df = hits_df.merge(hmm_info[['definition']], how='left', left_on=f"{db_name}_id", right_index=True)
+        hits_df.rename(columns={'definition': f"{db_name}_hits"}, inplace=True)
+    return hits_df
+
+
+def kofam_hmmscan_formater(hits:pd.DataFrame, hmm_info_path:str=None, use_dbcan2_thresholds:bool=False,
+                           top_hit:bool=True):
+    hmm_info = pd.read_csv(hmm_info_path, sep='\t', index_col=0)
+    if use_dbcan2_thresholds:
+        hits_sig = hits[hits.apply(get_sig_row, axis=1)]
+    else:
+        hits_sig = sig_scores(hits, hmm_info)
+    # if there are any significant results then parse to dataframe
+    if len(hits_sig) == 0:
+        return pd.DataFrame()
+    kegg_dict = dict()
+    for gene, frame in hits_sig.groupby('query_id'):
+        # TODO: take top hit for full length genes and all hits for domains?
+        # TODO: if top hit then give all e-value and bitscore info
+        if top_hit:
+            best_hit = frame[frame.full_evalue == frame.full_evalue.min()]
+            ko_id = best_hit['target_id'].iloc[0]
+            kegg_dict[gene] = [ko_id, hmm_info.loc[ko_id, 'definition']]
+        else:
+            kegg_dict[gene] = [','.join([i for i in frame.target_id]),
+                               '; '.join([hmm_info.loc[i, 'definition'] for i in frame.target_id])]
+    return pd.DataFrame(kegg_dict, index=['ko_id', 'kegg_hit']).transpose()
 
 
-def run_hmmscan_vogdb(genes_faa, vogdb_loc, output_loc, threads=10, db_handler=None, verbose=False):
-    # run hmmscan
-    vogdb_output = path.join(output_loc, 'vogdb_results.unprocessed.txt')
-    run_process(['hmmsearch', '--domtblout', vogdb_output, '--cpu', str(threads), vogdb_loc, genes_faa],
-                verbose=verbose)
 
-    # Process Results
-    vogdb_id_dict = dict()
-    vogdb_hit_dict = dict()
-    vogdb_category_dict = dict()
-
-    if path.isfile(vogdb_output) and stat(vogdb_output).st_size > 0:
-        vogdb_res = parse_hmmsearch_domtblout(vogdb_output)
-
-        significant = [row_num for row_num, row in vogdb_res.iterrows() if get_sig(row.target_start, row.target_end,
-                                                                                   row.target_length, row.full_evalue)]
-        if len(significant) != 0:  # if nothing significant then return nothing, don't get descriptions
-            vogdb_res = vogdb_res.loc[significant].sort_values('full_evalue')
-            vogdb_res_most_sig_list = list()
-            for gene, frame in vogdb_res.groupby('query_id'):
-                vogdb_res_most_sig_list.append(frame.index[0])
-            vogdb_res_most_sig = vogdb_res.loc[vogdb_res_most_sig_list]
-
-            if db_handler is not None:
-                vogdb_descriptions = db_handler.get_descriptions(set(vogdb_res_most_sig.target_id),
-                                                                 'vogdb_description')
-            else:
-                vogdb_descriptions = None
-            for _, row in vogdb_res_most_sig.iterrows():
-                gene = row['query_id']
-                vogdb_id = row['target_id']
-                if vogdb_descriptions is None:
-                    vogdb_hit_dict[gene] = vogdb_id
-                else:
-                    description = vogdb_descriptions.get(vogdb_id)
-                    categories_str = description.split('; ')[-1]
-                    vogdb_categories = [categories_str[0 + i:2 + i] for i in range(0, len(categories_str), 2)]
-                    vogdb_hit_dict[gene] = description
-                    vogdb_id_dict[gene] = description.split(' ')[0]
-                    vogdb_category_dict[gene] = ';'.join(set(vogdb_categories))
-
-    return pd.DataFrame((pd.Series(vogdb_id_dict, name='vogdb_id'),
-                         pd.Series(vogdb_category_dict, name='vogdb_categories'),
-                         pd.Series(vogdb_hit_dict, name='vogdb_hit'))).transpose()
+def vogdb_hmmscan_formater(hits:pd.DataFrame,  db_name:str, db_handler=None):
+    hits_sig = hits[hits.apply(get_sig_row, axis=1)]
+    if len(hits_sig) == 0:
+        # if nothing significant then return nothing, don't get descriptions
+        return pd.DataFrame()
+    # Get the best hits
+    hits_best = hits_sig.sort_values('full_evalue').drop_duplicates(subset=["query_id"])
+    if db_handler is None:
+        hits_df = hits_best[['target_id', 'query_id']].copy()
+    else:
+        # get_descriptions
+        desc_col = f"{db_name}_hits"
+        descriptions = pd.DataFrame(
+            db_handler.get_descriptions(hits_best['target_id'].unique(), f"{db_name}_description"),
+            index=[desc_col]).T
+        categories = descriptions[desc_col].apply(lambda x: x.split('; ')[-1])
+        descriptions[f"{db_name}_categories"] = categories.apply(
+            lambda x: ';'.join(set([x[i:i + 2] for i in range(0, len(x), 2)])))
+        descriptions['target_id'] = descriptions.index
+        hits_df = pd.merge(hits_best[['query_id', 'target_id']], descriptions, on=f'target_id')
+    hits_df.set_index('query_id', inplace=True, drop=True)
+    hits_df.rename_axis(None, inplace=True)
+    hits_df.rename(columns={'target_id': f"{db_name}_id"}, inplace=True)
+    return hits_df
+
+
+def run_hmmscan(genes_faa:str, db_loc:str, db_name:str, output_loc:str, formater:Callable,
+                threads:int=2, db_handler=None, verbose:bool=False):
+    output = path.join(output_loc, f'{db_name}_results.unprocessed.b6')
+    run_process(['hmmsearch', '--domtblout', output, '--cpu', str(threads), db_loc, genes_faa], verbose=verbose)
+    # Parse hmmsearch output
+    if not (path.isfile(output) and stat(output).st_size > 0):
+        return pd.DataFrame()
+    hits = parse_hmmsearch_domtblout(output)
+    return formater(hits)
 
 
 def get_gene_data(fasta_loc):
@@ -502,9 +493,11 @@ def annotate_gff(input_gff, output_gff, annotations, prefix=None):
             # get annotations to add from annotations file and add to end of line
             annotations_to_add = {strip_endings(i, ['_id']): annotations.loc[old_gene_name, i]
                                   for i in annotations.columns if i.endswith('_id')}
-            database_information = ['Dbxref="%s:%s"' % (key.strip(), value.strip())
-                                    for key, value in annotations_to_add.items()
-                                    if not pd.isna(value) and value != '']
+            database_information = ['Dbxref="%s:%s"' % (key.strip(), value.strip().replace(':', '_'))
+                                    for key, values in annotations_to_add.items()
+                                    if not pd.isna(values)
+                                    for value in values.split('; ')
+                                    if value != '']
             if len(database_information) > 0:
                 line += '%s;' % ';'.join(database_information)
         o.write('%s\n' % line)
@@ -664,7 +657,7 @@ def add_intervals_to_gff(annotations_loc, gff_loc, len_dict, interval_function,
             f.write(gff_intervals.write(io.StringIO(), format='gff3', seq_id=scaffold).getvalue())
 
 
-def do_blast_style_search(query_db, target_db, working_dir, db_handler, get_description, start_time,
+def do_blast_style_search(query_db, target_db, working_dir, db_handler, formater, start_time,
                           db_name='database', bit_score_threshold=60, rbh_bit_score_threshold=350, threads=10,
                           verbose=False):
     """A convenience function to do a blast style reciprocal best hits search"""
@@ -683,7 +676,7 @@ def do_blast_style_search(query_db, target_db, working_dir, db_handler, get_desc
         header_dict = db_handler.get_descriptions(hits['%s_hit' % db_name], '%s_description' % db_name)
     else:
         header_dict = multigrep(hits['%s_hit' % db_name], '%s_h' % target_db, '\x00', working_dir)
-    hits = get_description(hits, header_dict)
+    hits = formater(hits, header_dict)
     return hits
 
 
@@ -720,6 +713,31 @@ def process_custom_dbs(custom_fasta_loc, custom_db_name, output_dir, threads=1,
     return custom_db_locs
 
 
+def process_custom_hmms(custom_hmm_loc, custom_hmm_name, verbose=False):
+    if custom_hmm_loc is None:
+        custom_hmm_loc = ()
+    if custom_hmm_name is None:
+        custom_hmm_name = ()
+    if len(custom_hmm_loc) != len(custom_hmm_name):
+        raise ValueError('Lengths of custom db hmm list and custom hmm db name list must be the same.')
+    custom_hmm_locs = dict()
+    for i in range(len(custom_hmm_name)):
+        run_process(['hmmpress', '-f', custom_hmm_loc[i]], verbose=verbose)  # all are pressed just in case
+        custom_hmm_locs[custom_hmm_name[i]] = custom_hmm_loc[i]
+    return custom_hmm_locs
+
+def process_custom_hmm_cutoffs(custom_hmm_cutoffs_loc, custom_hmm_name, verbose=False):
+    if custom_hmm_cutoffs_loc is None:
+        return {}
+    if custom_hmm_name is None:
+        raise ValueError("You can't use the custom_hmm_cutoffs_loc argument without the custom_hmm_name and"
+                         " custom_hmm_locs aguments specified.")
+    if len(custom_hmm_cutoffs_loc) != len(custom_hmm_name):
+        warnings.warn(f"Custom hmm cutoffs and descriptions were only provided to the first {len(custom_hmm_cutoffs_loc)}."
+                      " The rest of the custom hmms will use standard cutoffs and have no descriptions.")
+    return {custom_hmm_name[i]:j for i, j in enumerate(custom_hmm_cutoffs_loc)}
+
+
 class Annotation:
     def __init__(self, name, scaffolds, genes_faa, genes_fna, gff, gbk, annotations, trnas, rrnas):
         # TODO: get abspath for every input file/dir
@@ -743,9 +761,9 @@ def get_trnas(self):
     def get_rrnas(self):
         return pd.read_csv(self.rrnas_loc, sep='\t')
 
-
-def annotate_orfs(gene_faa, db_locs, tmp_dir, start_time, db_handler, custom_db_locs=(), bit_score_threshold=60,
-                  rbh_bit_score_threshold=350, kofam_use_dbcan2_thresholds=False, threads=10, verbose=False):
+def annotate_orfs(gene_faa, db_handler, tmp_dir, start_time, custom_db_locs=(), custom_hmm_locs=(),
+                  custom_hmm_cutoffs_locs=(), bit_score_threshold=60, rbh_bit_score_threshold=350,
+                  kofam_use_dbcan2_thresholds=False, threads=10, verbose=False):
     # run reciprocal best hits searches
     print('%s: Turning genes from prodigal to mmseqs2 db' % str(datetime.now() - start_time))
     query_db = path.join(tmp_dir, 'gene.mmsdb')
@@ -754,59 +772,86 @@ def annotate_orfs(gene_faa, db_locs, tmp_dir, start_time, db_handler, custom_db_
     annotation_list = list()
 
     # Get kegg hits
-    if db_locs.get('kegg') is not None:
-        annotation_list.append(do_blast_style_search(query_db, db_locs['kegg'], tmp_dir,
+    if db_handler.db_locs.get('kegg') is not None:
+        #TODO Change the get_kegg_description name in function do_blast_style_search to formater
+        #TODO think about how this can be consitent with blast and mmseqs
+        annotation_list.append(do_blast_style_search(query_db, db_handler.db_locs['kegg'], tmp_dir,
                                                      db_handler, get_kegg_description, start_time,
                                                      'kegg', bit_score_threshold, rbh_bit_score_threshold, threads,
                                                      verbose))
-    elif db_locs.get('kofam') is not None and db_locs.get('kofam_ko_list') is not None:
+    elif db_handler.db_locs.get('kofam') is not None and db_handler.db_locs.get('kofam_ko_list') is not None:
         print('%s: Getting hits from kofam' % str(datetime.now() - start_time))
-        annotation_list.append(run_hmmscan_kofam(gene_faa, db_locs['kofam'], tmp_dir,
-                                                 pd.read_csv(db_locs['kofam_ko_list'], sep='\t', index_col=0),
-                                                 use_dbcan2_thresholds=kofam_use_dbcan2_thresholds, threads=threads,
-                                                 verbose=verbose))
+        annotation_list.append(run_hmmscan(genes_faa=gene_faa,
+                                           db_loc=db_handler.db_locs['kofam'],
+                                           db_name='kofam',
+                                           output_loc=tmp_dir, #check_impliments
+                                           threads=threads, #check_impliments
+                                           verbose=verbose,
+                                           formater=partial(
+                                               kofam_hmmscan_formater,
+                                               hmm_info_path=db_handler.db_locs['kofam_ko_list'],
+                                               top_hit=True,
+                                               use_dbcan2_thresholds=kofam_use_dbcan2_thresholds
+                                           )))
     else:
         warnings.warn('No KEGG source provided so distillation will be of limited use.')
 
     # Get uniref hits
-    if db_locs.get('uniref') is not None:
-        annotation_list.append(do_blast_style_search(query_db, db_locs['uniref'], tmp_dir,
+    if db_handler.db_locs.get('uniref') is not None:
+        annotation_list.append(do_blast_style_search(query_db, db_handler.db_locs['uniref'], tmp_dir,
                                                      db_handler, get_uniref_description,
                                                      start_time, 'uniref', bit_score_threshold,
                                                      rbh_bit_score_threshold, threads, verbose))
 
     # Get viral hits
-    if db_locs.get('viral') is not None:
+    if db_handler.db_locs.get('viral') is not None:
         get_viral_description = partial(get_basic_description, db_name='viral')
-        annotation_list.append(do_blast_style_search(query_db, db_locs['viral'], tmp_dir,
+        annotation_list.append(do_blast_style_search(query_db, db_handler.db_locs['viral'], tmp_dir,
                                                      db_handler, get_viral_description,
                                                      start_time, 'viral', bit_score_threshold,
                                                      rbh_bit_score_threshold, threads, verbose))
 
     # Get peptidase hits
-    if db_locs.get('peptidase') is not None:
-        annotation_list.append(do_blast_style_search(query_db, db_locs['peptidase'], tmp_dir,
+    if db_handler.db_locs.get('peptidase') is not None:
+        annotation_list.append(do_blast_style_search(query_db, db_handler.db_locs['peptidase'], tmp_dir,
                                                      db_handler, get_peptidase_description,
                                                      start_time, 'peptidase', bit_score_threshold,
                                                      rbh_bit_score_threshold, threads, verbose))
 
     # Get pfam hits
-    if db_locs.get('pfam') is not None:
+    if db_handler.db_locs.get('pfam') is not None:
         print('%s: Getting hits from pfam' % str(datetime.now() - start_time))
-        annotation_list.append(run_mmseqs_profile_search(query_db, db_locs['pfam'], tmp_dir, output_prefix='pfam',
-                                                         db_handler=db_handler, threads=threads, verbose=verbose))
+        annotation_list.append(run_mmseqs_profile_search(query_db, db_handler.db_locs['pfam'], tmp_dir,
+                                                         output_prefix='pfam', db_handler=db_handler, threads=threads,
+                                                         verbose=verbose))
 
     # use hmmer to detect cazy ids using dbCAN
-    if db_locs.get('dbcan') is not None:
+    if db_handler.db_locs.get('dbcan') is not None:
         print('%s: Getting hits from dbCAN' % str(datetime.now() - start_time))
-        annotation_list.append(run_hmmscan_dbcan(gene_faa, db_locs['dbcan'], tmp_dir, threads, db_handler=db_handler,
-                                                 verbose=verbose))
+        annotation_list.append(run_hmmscan(genes_faa=gene_faa,
+                                           db_loc=db_handler.db_locs['dbcan'],
+                                           db_name='cazy',
+                                           output_loc=tmp_dir,
+                                           threads=threads,
+                                           formater=partial(
+                                               dbcan_hmmscan_formater,
+                                               db_name='cazy',
+                                               db_handler=db_handler
+                                           )))
 
     # use hmmer to detect vogdbs
-    if db_locs.get('vogdb') is not None:
+    if db_handler.db_locs.get('vogdb') is not None:
         print('%s: Getting hits from VOGDB' % str(datetime.now() - start_time))
-        annotation_list.append(run_hmmscan_vogdb(gene_faa, db_locs['vogdb'], tmp_dir, threads, db_handler=db_handler,
-                                                 verbose=verbose))
+        annotation_list.append(run_hmmscan(genes_faa=gene_faa,
+                                           db_loc=db_handler.db_locs['vogdb'],
+                                           db_name='vogdb',
+                                           threads=threads,
+                                           output_loc=tmp_dir,
+                                           formater=partial(
+                                               vogdb_hmmscan_formater,
+                                               db_name='vogdb',
+                                               db_handler=db_handler
+                                           )))
 
     for db_name, db_loc in custom_db_locs.items():
         print('%s: Getting hits from %s' % (str(datetime.now() - start_time), db_name))
@@ -816,6 +861,20 @@ def annotate_orfs(gene_faa, db_locs, tmp_dir, start_time, db_handler, custom_db_
                                                      bit_score_threshold, rbh_bit_score_threshold, threads,
                                                      verbose))
 
+    # get hits to hmm style custom databases
+    for hmm_name, hmm_loc in custom_hmm_locs.items():
+        annotation_list.append(run_hmmscan(genes_faa=gene_faa,
+                                           db_loc=hmm_loc,
+                                           db_name=hmm_name,
+                                           threads=threads,
+                                           output_loc=tmp_dir,
+                                           formater=partial(
+                                               generic_hmmscan_formater,
+                                               db_name=hmm_name,
+                                               hmm_info_path=custom_hmm_cutoffs_locs.get(hmm_name),
+                                               top_hit=True
+                                           )))
+
     # heme regulatory motif count
     annotation_list.append(pd.Series(count_motifs(gene_faa, '(C..CH)'), name='heme_regulatory_motif_count'))
 
@@ -829,10 +888,11 @@ def annotate_orfs(gene_faa, db_locs, tmp_dir, start_time, db_handler, custom_db_
     return annotations
 
 
-def annotate_fasta(fasta_loc, fasta_name, output_dir, db_locs, db_handler, min_contig_size=5000, prodigal_mode='meta',
-                   trans_table='11', custom_db_locs=(), bit_score_threshold=60, rbh_bit_score_threshold=350,
-                   kofam_use_dbcan2_thresholds=False, skip_trnascan=False, start_time=datetime.now(), threads=1,
-                   rename_bins=True, keep_tmp_dir=False, verbose=False):
+def annotate_fasta(fasta_loc, fasta_name, output_dir, db_handler, min_contig_size=5000, prodigal_mode='meta',
+                   trans_table='11', custom_db_locs=(), custom_hmm_locs=(), custom_hmm_cutoffs_locs=(),
+                   bit_score_threshold=60, rbh_bit_score_threshold=350, kofam_use_dbcan2_thresholds=False,
+                   skip_trnascan=False, start_time=datetime.now(), threads=1, rename_bins=True, keep_tmp_dir=False,
+                   verbose=False):
     """Annotated a single multifasta file, all file based outputs will be in output_dir"""
     # make temporary directory
     tmp_dir = path.join(output_dir, 'tmp')
@@ -852,8 +912,9 @@ def annotate_fasta(fasta_loc, fasta_name, output_dir, db_locs, db_handler, min_c
                                                 verbose=verbose)
 
     # annotate ORFs
-    annotations = annotate_orfs(gene_faa, db_locs, tmp_dir, start_time, db_handler, custom_db_locs, bit_score_threshold,
-                                rbh_bit_score_threshold, kofam_use_dbcan2_thresholds, threads, verbose)
+    annotations = annotate_orfs(gene_faa, db_handler, tmp_dir, start_time, custom_db_locs, custom_hmm_locs,
+                                custom_hmm_cutoffs_locs, bit_score_threshold, rbh_bit_score_threshold,
+                                kofam_use_dbcan2_thresholds, threads, verbose)
     annotations = pd.concat([get_gene_data(gene_faa), annotations], axis=1, sort=False)
 
     renamed_scaffolds = path.join(output_dir, 'scaffolds.annotated.fa')
@@ -914,37 +975,15 @@ def annotate_fasta(fasta_loc, fasta_name, output_dir, db_locs, db_handler, min_c
                       gff=renamed_gffs, gbk=current_gbk, annotations=annotations_loc, trnas=trna_loc, rrnas=rrna_loc)
 
 
-def filter_db_locs(db_locs, low_mem_mode=False, use_uniref=False, use_vogdb=False, master_list=MAG_DBS_TO_ANNOTATE):
-    dbs_to_use = master_list
-    # filter out dbs for low mem mode
-    if low_mem_mode:
-        if ('kofam' not in db_locs) or ('kofam_ko_list' not in db_locs):
-            raise ValueError('To run in low memory mode KOfam must be configured for use in DRAM')
-        dbs_to_use = [i for i in dbs_to_use if i not in ('uniref', 'kegg', 'vogdb')]
-    # check on uniref status
-    if use_uniref:
-        if 'uniref' not in db_locs:
-            warnings.warn('Sequences will not be annoated against uniref as it is not configured for use in DRAM')
-    else:
-        dbs_to_use = [i for i in dbs_to_use if i != 'uniref']
-    # check on vogdb status
-    if use_vogdb:
-        if 'vogdb' not in db_locs:
-            warnings.warn('Sequences will not be annoated against VOGDB as it is not configured for use in DRAM')
-    else:
-        dbs_to_use = [i for i in dbs_to_use if i != 'vogdb']
-    db_locs = {key: value for key, value in db_locs.items() if key in dbs_to_use}
-    return db_locs
-
-
 def get_fasta_name(fasta_loc):
     return path.splitext(path.basename(remove_suffix(fasta_loc, '.gz')))[0]
 
 
-def annotate_fastas(fasta_locs, output_dir, db_locs, db_handler, min_contig_size=5000, prodigal_mode='meta',
-                    trans_table='11', bit_score_threshold=60, rbh_bit_score_threshold=350, custom_db_name=(),
-                    custom_fasta_loc=(), kofam_use_dbcan2_thresholds=False, skip_trnascan=False, rename_bins=True,
-                    keep_tmp_dir=True, start_time=datetime.now(), threads=10, verbose=True):
+def annotate_fastas(fasta_locs, output_dir, db_handler, min_contig_size=5000, prodigal_mode='meta', trans_table='11',
+                    bit_score_threshold=60, rbh_bit_score_threshold=350, custom_db_name=(), custom_fasta_loc=(),
+                    custom_hmm_name=(), custom_hmm_loc=(), custom_hmm_cutoffs_loc=(), kofam_use_dbcan2_thresholds=False,
+                    skip_trnascan=False, rename_bins=True, keep_tmp_dir=True, start_time=datetime.now(), threads=10,
+                    verbose=True):
     # check for no conflicting options/configurations
     tmp_dir = path.join(output_dir, 'working_dir')
     mkdir(tmp_dir)
@@ -952,6 +991,8 @@ def annotate_fastas(fasta_locs, output_dir, db_locs, db_handler, min_contig_size
     # setup custom databases to be searched
     custom_db_locs = process_custom_dbs(custom_fasta_loc, custom_db_name, path.join(tmp_dir, 'custom_dbs'), threads,
                                         verbose)
+    custom_hmm_locs = process_custom_hmms(custom_hmm_loc, custom_hmm_name)
+    custom_hmm_cutoffs_locs= process_custom_hmm_cutoffs(custom_hmm_cutoffs_loc, custom_hmm_name)
     print('%s: Retrieved database locations and descriptions' % (str(datetime.now() - start_time)))
 
     # iterate over list of fastas and annotate each individually
@@ -962,10 +1003,11 @@ def annotate_fastas(fasta_locs, output_dir, db_locs, db_handler, min_contig_size
         print('%s: Annotating %s' % (str(datetime.now() - start_time), fasta_name))
         fasta_dir = path.join(tmp_dir, fasta_name)
         mkdir(fasta_dir)
-        annotations_list.append(annotate_fasta(fasta_loc, fasta_name, fasta_dir, db_locs, db_handler, min_contig_size,
-                                               prodigal_mode, trans_table, custom_db_locs, bit_score_threshold,
-                                               rbh_bit_score_threshold, kofam_use_dbcan2_thresholds, skip_trnascan,
-                                               start_time, threads, rename_bins, keep_tmp_dir, verbose))
+        annotations_list.append(
+            annotate_fasta(fasta_loc, fasta_name, fasta_dir, db_handler, min_contig_size, prodigal_mode, trans_table,
+                           custom_db_locs, custom_hmm_locs, custom_hmm_cutoffs_locs, bit_score_threshold,
+                           rbh_bit_score_threshold, kofam_use_dbcan2_thresholds, skip_trnascan, start_time, threads,
+                           rename_bins, keep_tmp_dir, verbose))
     print('%s: Annotations complete, processing annotations' % str(datetime.now() - start_time))
 
     all_annotations = merge_annotations(annotations_list, output_dir)
@@ -978,9 +1020,9 @@ def annotate_fastas(fasta_locs, output_dir, db_locs, db_handler, min_contig_size
 
 def annotate_bins_cmd(input_fasta, output_dir='.', min_contig_size=5000, prodigal_mode='meta', trans_table='11',
                       bit_score_threshold=60, rbh_bit_score_threshold=350, custom_db_name=(), custom_fasta_loc=(),
-                      use_uniref=False, use_vogdb=False, kofam_use_dbcan2_thresholds=False, skip_trnascan=False,
-                      gtdb_taxonomy=(), checkm_quality=(), keep_tmp_dir=True, low_mem_mode=False, threads=10,
-                      verbose=True):
+                      custom_hmm_name=(), custom_hmm_loc=(), custom_hmm_cutoffs_loc=(), use_uniref=False, use_vogdb=False,
+                      kofam_use_dbcan2_thresholds=False, skip_trnascan=False, gtdb_taxonomy=(), checkm_quality=(),
+                      keep_tmp_dir=True, low_mem_mode=False, threads=10, verbose=True):
     fasta_locs = [j for i in input_fasta for j in glob(i)]
     if len(fasta_locs) == 0:
         raise ValueError('Given fasta locations return no paths: %s' % input_fasta)
@@ -990,9 +1032,9 @@ def annotate_bins_cmd(input_fasta, output_dir='.', min_contig_size=5000, prodiga
     print('%s fastas found' % len(fasta_locs))
     rename_bins = True
     annotate_bins(list(set(fasta_locs)), output_dir, min_contig_size, prodigal_mode, trans_table, bit_score_threshold,
-                  rbh_bit_score_threshold, custom_db_name, custom_fasta_loc, use_uniref, use_vogdb,
-                  kofam_use_dbcan2_thresholds, skip_trnascan, gtdb_taxonomy, checkm_quality, rename_bins, keep_tmp_dir,
-                  low_mem_mode, threads, verbose)
+                  rbh_bit_score_threshold, custom_db_name, custom_fasta_loc, custom_hmm_name, custom_hmm_loc,
+                  custom_hmm_cutoffs_loc, use_uniref, use_vogdb, kofam_use_dbcan2_thresholds, skip_trnascan,
+                  gtdb_taxonomy, checkm_quality, rename_bins, keep_tmp_dir, low_mem_mode, threads, verbose)
 
 
 # TODO: Add force flag to remove output dir if it already exists
@@ -1000,9 +1042,9 @@ def annotate_bins_cmd(input_fasta, output_dir='.', min_contig_size=5000, prodiga
 # TODO: make fasta loc either a string or list to remove annotate_bins_cmd and annotate_called_genes_cmd?
 def annotate_bins(fasta_locs, output_dir='.', min_contig_size=2500, prodigal_mode='meta', trans_table='11',
                   bit_score_threshold=60, rbh_bit_score_threshold=350, custom_db_name=(), custom_fasta_loc=(),
-                  use_uniref=False, use_vogdb=False, kofam_use_dbcan2_thresholds=False, skip_trnascan=False,
-                  gtdb_taxonomy=(), checkm_quality=(), rename_bins=True, keep_tmp_dir=True, low_mem_mode=False,
-                  threads=10, verbose=True):
+                  custom_hmm_name=(), custom_hmm_loc=(), custom_hmm_cutoffs_loc=(), use_uniref=False, use_vogdb=False,
+                  kofam_use_dbcan2_thresholds=False, skip_trnascan=False, gtdb_taxonomy=(), checkm_quality=(),
+                  rename_bins=True, keep_tmp_dir=True, low_mem_mode=False, threads=10, verbose=True):
     # set up
     start_time = datetime.now()
     print('%s: Annotation started' % str(datetime.now()))
@@ -1023,16 +1065,16 @@ def annotate_bins(fasta_locs, output_dir='.', min_contig_size=2500, prodigal_mod
         raise ValueError('Prodigal translation table must be 1-25')
 
     # get database locations
-    db_locs = get_database_locs()
-    db_handler = DatabaseHandler(db_locs['description_db'])
-    db_locs = filter_db_locs(db_locs, low_mem_mode, use_uniref, use_vogdb)
+    db_handler = DatabaseHandler()
+    db_handler.filter_db_locs(low_mem_mode, use_uniref, use_vogdb, master_list=MAG_DBS_TO_ANNOTATE)
 
     mkdir(output_dir)
 
-    all_annotations = annotate_fastas(fasta_locs, output_dir, db_locs, db_handler, min_contig_size, prodigal_mode,
-                                      trans_table, bit_score_threshold, rbh_bit_score_threshold, custom_db_name,
-                                      custom_fasta_loc, kofam_use_dbcan2_thresholds, skip_trnascan, rename_bins,
-                                      keep_tmp_dir, start_time, threads, verbose)
+    all_annotations = annotate_fastas(fasta_locs, output_dir, db_handler, min_contig_size, prodigal_mode, trans_table,
+                                      bit_score_threshold, rbh_bit_score_threshold, custom_db_name, custom_fasta_loc,
+                                      custom_hmm_name, custom_hmm_loc, custom_hmm_cutoffs_loc,
+                                      kofam_use_dbcan2_thresholds, skip_trnascan, rename_bins, keep_tmp_dir, start_time,
+                                      threads, verbose)
     # if given add taxonomy information
     if len(gtdb_taxonomy) > 0:
         gtdb_taxonomy = pd.concat([pd.read_csv(i, sep='\t', index_col=0) for i in gtdb_taxonomy])
@@ -1076,7 +1118,8 @@ def annotate_bins(fasta_locs, output_dir='.', min_contig_size=2500, prodigal_mod
 
 
 def annotate_called_genes_cmd(input_faa, output_dir='.', bit_score_threshold=60, rbh_bit_score_threshold=350,
-                              custom_db_name=(), custom_fasta_loc=(), use_uniref=False, use_vogdb=False,
+                              custom_db_name=(), custom_fasta_loc=(), custom_hmm_loc=(), custom_hmm_name=(),
+                              custom_hmm_cutoffs_loc=(), use_uniref=False, use_vogdb=False,
                               kofam_use_dbcan2_thresholds=False, rename_genes=True, keep_tmp_dir=True,
                               low_mem_mode=False, threads=10, verbose=True):
     fasta_locs = glob(input_faa)
@@ -1084,12 +1127,13 @@ def annotate_called_genes_cmd(input_faa, output_dir='.', bit_score_threshold=60,
         raise ValueError('Given fasta locations returns no paths: %s' % input_faa)
     print('%s fastas found' % len(fasta_locs))
     annotate_called_genes(fasta_locs, output_dir, bit_score_threshold, rbh_bit_score_threshold, custom_db_name,
-                          custom_fasta_loc, use_uniref, use_vogdb, kofam_use_dbcan2_thresholds, rename_genes,
-                          keep_tmp_dir, low_mem_mode, threads, verbose)
+                          custom_fasta_loc, custom_hmm_loc, custom_hmm_name, custom_hmm_cutoffs_loc, use_uniref, use_vogdb,
+                          kofam_use_dbcan2_thresholds, rename_genes, keep_tmp_dir, low_mem_mode, threads, verbose)
 
 
 def annotate_called_genes(fasta_locs, output_dir='.', bit_score_threshold=60, rbh_bit_score_threshold=350,
-                          custom_db_name=(), custom_fasta_loc=(), use_uniref=False, use_vogdb=False,
+                          custom_db_name=(), custom_fasta_loc=(), custom_hmm_loc=(), custom_hmm_name=(),
+                          custom_hmm_cutoffs_loc=(), use_uniref=False, use_vogdb=False,
                           kofam_use_dbcan2_thresholds=False, rename_genes=True, keep_tmp_dir=True, low_mem_mode=False,
                           threads=10, verbose=True):
     # set up
@@ -1097,9 +1141,8 @@ def annotate_called_genes(fasta_locs, output_dir='.', bit_score_threshold=60, rb
     print('%s: Annotation started' % str(datetime.now()))
 
     # get database locations
-    db_locs = get_database_locs()
-    db_handler = DatabaseHandler(db_locs['description_db'])
-    db_locs = filter_db_locs(db_locs, low_mem_mode, use_uniref, use_vogdb)
+    db_handler = DatabaseHandler()
+    db_handler.filter_db_locs(low_mem_mode, use_uniref, use_vogdb, master_list=MAG_DBS_TO_ANNOTATE)
 
     mkdir(output_dir)
     tmp_dir = path.join(output_dir, 'working_dir')
@@ -1108,6 +1151,8 @@ def annotate_called_genes(fasta_locs, output_dir='.', bit_score_threshold=60, rb
     # setup custom databases to be searched
     custom_db_locs = process_custom_dbs(custom_fasta_loc, custom_db_name, path.join(tmp_dir, 'custom_dbs'), threads,
                                         verbose)
+    custom_hmm_locs = process_custom_hmms(custom_hmm_loc, custom_hmm_name)
+    custom_hmm_cutoffs_locs= process_custom_hmm_cutoffs(custom_hmm_cutoffs_loc, custom_hmm_name)
     print('%s: Retrieved database locations and descriptions' % (str(datetime.now() - start_time)))
 
     # annotate
@@ -1120,9 +1165,9 @@ def annotate_called_genes(fasta_locs, output_dir='.', bit_score_threshold=60, rb
         mkdir(fasta_dir)
 
         # annotate
-        annotations = annotate_orfs(fasta_loc, db_locs, fasta_dir, start_time, db_handler, custom_db_locs,
-                                    bit_score_threshold, rbh_bit_score_threshold, kofam_use_dbcan2_thresholds, threads,
-                                    verbose)
+        annotations = annotate_orfs(fasta_loc, db_handler, fasta_dir, start_time, custom_db_locs, custom_hmm_locs,
+                                    custom_hmm_cutoffs_locs, bit_score_threshold, rbh_bit_score_threshold,
+                                    kofam_use_dbcan2_thresholds, threads, verbose)
 
         annotated_faa = path.join(fasta_dir, 'genes.faa')
         create_annotated_fasta(fasta_loc, annotations, annotated_faa, name=fasta_name)
diff --git a/mag_annotator/annotate_vgfs.py b/mag_annotator/annotate_vgfs.py
index 4b293ede..b1cfd272 100644
--- a/mag_annotator/annotate_vgfs.py
+++ b/mag_annotator/annotate_vgfs.py
@@ -9,8 +9,8 @@
 from skbio.io import write as write_sequence
 
 from mag_annotator.database_handler import DatabaseHandler
-from mag_annotator.annotate_bins import filter_db_locs, annotate_fastas
-from mag_annotator.utils import get_database_locs, get_ids_from_annotation
+from mag_annotator.annotate_bins import annotate_fastas
+from mag_annotator.utils import get_ids_from_annotation
 
 VMAG_DBS_TO_ANNOTATE = ('kegg', 'kofam', 'kofam_ko_list', 'uniref', 'peptidase', 'pfam', 'dbcan', 'viral', 'vogdb')
 VIRSORTER_COLUMN_NAMES = ['gene_name', 'start_position', 'end_position', 'length', 'strandedness',
@@ -43,7 +43,7 @@
                            'T01B', 'T01A', 'T03', 'U32', 'U40'}
 
 
-def is_affi_tab_not_fasta(input_file:str) -> bool:
+def is_affi_tab_not_fasta(input_file: str) -> bool:
     """
     Checks that a file is a VIRSorter affi tab file.
 
@@ -54,7 +54,7 @@ def is_affi_tab_not_fasta(input_file:str) -> bool:
     """
     with open(input_file) as input_io:
         num_delimiter = input_io.readlines()[1].count('|')
-    if num_delimiter == 11 :
+    if num_delimiter == 11:
         return True
     return False
 
@@ -77,7 +77,7 @@ def remove_bad_chars_fasta(fasta):
     return new_seqs
 
 
-def remove_bad_chars_virsorter_affi_contigs(virsorter_in:str) -> str:
+def remove_bad_chars_virsorter_affi_contigs(virsorter_in: str) -> str:
     """
     Determine and call the appropriate function to cleans incompatible characters from VIRSorter
     or VIRSorter2 affi headers.
@@ -108,7 +108,7 @@ def remove_bad_chars_virsorter_affi_contigs(virsorter_in:str) -> str:
     return '%s\n' % '\n'.join(new_lines)
 
 
-def remove_bad_chars(input_fasta:str=None, input_virsorter_affi_contigs:str=None, output:str=None):
+def remove_bad_chars(input_fasta: str = None, input_virsorter_affi_contigs: str = None, output: str = None):
     """
     Determine and call the function to cleans incompatible characters from headers of virsorter
     affi or fasta file.
@@ -376,10 +376,10 @@ def get_virsorter_affi_contigs_name(scaffold):
     return get_virsorter_original_affi_contigs_name(scaffold)
 
 
-def add_dramv_scores_and_flags(annotations, db_locs=None, virsorter_hits=None, input_fasta=None):
+def add_dramv_scores_and_flags(annotations, db_handler, virsorter_hits=None, input_fasta=None):
     # setting up scoring viral genes
-    amg_database_frame = pd.read_csv(db_locs['amg_database'], sep='\t')
-    genome_summary_form = pd.read_csv(db_locs['genome_summary_form'], sep='\t', index_col=0)
+    amg_database_frame = pd.read_csv(db_handler.dram_sheet_locs['amg_database'], sep='\t')
+    genome_summary_form = pd.read_csv(db_handler.dram_sheet_locs['genome_summary_form'], sep='\t', index_col=0)
     genome_summary_form = genome_summary_form.loc[genome_summary_form.potential_amg]
 
     # add auxiliary score
@@ -420,8 +420,9 @@ def add_dramv_scores_and_flags(annotations, db_locs=None, virsorter_hits=None, i
 
 def annotate_vgfs(input_fasta, virsorter_affi_contigs=None, output_dir='.', min_contig_size=2500, split_contigs=False,
                   prodigal_mode='meta', trans_table='11', bit_score_threshold=60, rbh_bit_score_threshold=350,
-                  custom_db_name=(), custom_fasta_loc=(), use_uniref=False, kofam_use_dbcan2_thresholds=False,
-                  skip_trnascan=False, keep_tmp_dir=True, low_mem_mode=False, threads=10, verbose=True):
+                  custom_db_name=(), custom_fasta_loc=(), custom_hmm_loc=(), custom_hmm_cutoffs_loc=(),
+                  custom_hmm_name=(), use_uniref=False, kofam_use_dbcan2_thresholds=False, skip_trnascan=False,
+                  keep_tmp_dir=True, low_mem_mode=False, threads=10, verbose=True):
     # set up
     start_time = datetime.now()
     print('%s: Viral annotation started' % str(datetime.now()))
@@ -436,18 +437,17 @@ def annotate_vgfs(input_fasta, virsorter_affi_contigs=None, output_dir='.', min_
                       'training to work well.')
 
     # get database locations
-    db_locs = get_database_locs()
-    db_handler = DatabaseHandler(db_locs['description_db'])
-    db_locs_anno = filter_db_locs(db_locs, low_mem_mode, use_uniref, True, VMAG_DBS_TO_ANNOTATE)
+    db_handler = DatabaseHandler()
+    db_handler.filter_db_locs(low_mem_mode, use_uniref, True, VMAG_DBS_TO_ANNOTATE)
 
     if virsorter_affi_contigs is not None:
         virsorter_hits = get_virsorter_hits(virsorter_affi_contigs)
     else:
         virsorter_hits = None
 
+    mkdir(output_dir)
     if split_contigs:
         # split sequences into separate fastas
-        mkdir(output_dir)
         contig_dir = path.join(output_dir, 'vMAGs')
         mkdir(contig_dir)
 
@@ -455,8 +455,8 @@ def annotate_vgfs(input_fasta, virsorter_affi_contigs=None, output_dir='.', min_
         for seq in read_sequence(input_fasta, format='fasta'):
             if len(seq) >= min_contig_size:
                 if '=' in seq.metadata['id'] or ';' in seq.metadata['id']:
-                    raise ValueError('FASTA headers must not have = or ; before the first space (%s). To run DRAM-v you '
-                                     'must rerun VIRSorter with = and ; removed from the headers or run DRAM-v.py '
+                    raise ValueError('FASTA headers must not have = or ; before the first space (%s). To run DRAM-v '
+                                     'you must rerun VIRSorter with = and ; removed from the headers or run DRAM-v.py '
                                      'remove_bad_characters and then rerun DRAM-v' % seq.metadata['id'])
                 if virsorter_hits is not None:
                     if get_virsorter_affi_contigs_name(seq.metadata['id']) not in virsorter_hits['name'].values:
@@ -470,15 +470,16 @@ def annotate_vgfs(input_fasta, virsorter_affi_contigs=None, output_dir='.', min_
 
     # annotate vMAGs
     rename_bins = False
-    annotations = annotate_fastas(contig_locs, output_dir, db_locs_anno, db_handler, min_contig_size, prodigal_mode,
-                                  trans_table, bit_score_threshold, rbh_bit_score_threshold, custom_db_name,
-                                  custom_fasta_loc, kofam_use_dbcan2_thresholds, skip_trnascan, rename_bins,
-                                  keep_tmp_dir, start_time, threads, verbose)
+    annotations = annotate_fastas(contig_locs, output_dir, db_handler, min_contig_size, prodigal_mode, trans_table,
+                                  bit_score_threshold, rbh_bit_score_threshold, custom_db_name, custom_fasta_loc,
+                                  custom_hmm_name, custom_hmm_loc, custom_hmm_cutoffs_loc, kofam_use_dbcan2_thresholds,
+                                  skip_trnascan, rename_bins, keep_tmp_dir, start_time, threads, verbose)
     print('%s: Annotations complete, assigning auxiliary scores and flags' % str(datetime.now() - start_time))
 
-    annotations = add_dramv_scores_and_flags(annotations, db_locs, virsorter_hits, input_fasta)
+    annotations = add_dramv_scores_and_flags(annotations, db_handler, virsorter_hits, input_fasta)
 
     # write annotations
     annotations.to_csv(path.join(output_dir, 'annotations.tsv'), sep='\t')
 
     print("%s: Completed annotations" % str(datetime.now() - start_time))
+
diff --git a/mag_annotator/database_handler.py b/mag_annotator/database_handler.py
index 05acd6af..f0063672 100644
--- a/mag_annotator/database_handler.py
+++ b/mag_annotator/database_handler.py
@@ -1,21 +1,60 @@
-from os import path
+from os import path, remove
 from warnings import warn
+from pkg_resources import resource_filename
+import json
+import gzip
+from shutil import copy2
+import warnings
 
 from sqlalchemy import create_engine
 from sqlalchemy.orm import sessionmaker
 
-from mag_annotator.database_setup import TABLE_NAME_TO_CLASS_DICT
+import pandas as pd
+
+from mag_annotator.database_setup import TABLE_NAME_TO_CLASS_DICT, create_description_db
 from mag_annotator.utils import divide_chunks
 
+SEARCH_DATABASES = ('kegg', 'kofam', 'kofam_ko_list', 'uniref', 'pfam', 'dbcan', 'viral', 'peptidase', 'vogdb')
+DRAM_SHEETS = ('genome_summary_form', 'module_step_form', 'etc_module_database', 'function_heatmap_form',
+               'amg_database')
+DATABASE_DESCRIPTIONS = ('pfam_hmm_dat', 'dbcan_fam_activities', 'vog_annotations')
 
 # TODO: store all sequence db locs within database handler class
+# TODO: store scoring information here e.g. bitscore_threshold, hmm cutoffs
+# TODO: set up custom databases here
+# TODO: in advanced config separate search databases, search database description files, description db, DRAM sheets
+# TODO: ko_list should be parsed into the DB and stored as a database description file and not a search database
+
+
+def get_config_loc():
+    return path.abspath(resource_filename('mag_annotator', 'CONFIG'))
 
 
 class DatabaseHandler:
-    def __init__(self, database_loc):
-        if not path.exists(database_loc):
-            raise ValueError('Database does not exist at path %s' % database_loc)
-        engine = create_engine('sqlite:///%s' % database_loc)
+    def __init__(self, config_loc=None):
+        if config_loc is None:
+            config_loc = get_config_loc()
+
+        # read in configuration # TODO: validate config file after reading it in
+        self.config_loc = config_loc
+        config = json.loads(open(config_loc).read())
+        self.db_locs = {key: value for key, value in config.items() if key in SEARCH_DATABASES}
+        self.db_description_locs = {key: value for key, value in config.items() if key in DATABASE_DESCRIPTIONS}
+        self.dram_sheet_locs = {key: value for key, value in config.items() if key in DRAM_SHEETS}
+
+        # set up description database connection
+        self.description_loc = config.get('description_db')
+        if self.description_loc is None:
+            self.session = None
+            warnings.warn('Database does not exist at path %s' % self.description_loc)
+        elif not path.exists(self.description_loc):
+            self.session = None
+            warnings.warn('Database does not exist at path %s' % self.description_loc)
+        else:
+            self.start_db_session()
+
+    def start_db_session(self):
+        engine = create_engine('sqlite:///%s' % self.description_loc)
         db_session = sessionmaker(bind=engine)
         self.session = db_session()
 
@@ -24,7 +63,7 @@ def add_descriptions_to_database(self, description_list, db_name, clear_table=Tr
         description_class = TABLE_NAME_TO_CLASS_DICT[db_name]
         if clear_table:
             self.session.query(description_class).delete()
-        self.session.bulk_save_objects([description_class(**i) for i in description_list])
+        self.session.bulk_save_objects([description_class(**i) for i in description_list])  # TODO: try batching
         self.session.commit()
         self.session.expunge_all()
 
@@ -46,3 +85,260 @@ def get_descriptions(self, ids, db_name):
     @staticmethod
     def get_database_names():
         return TABLE_NAME_TO_CLASS_DICT.keys()
+
+    def set_database_paths(self, kegg_db_loc=None, kofam_hmm_loc=None, kofam_ko_list_loc=None, uniref_db_loc=None,
+                           pfam_db_loc=None, pfam_hmm_dat=None, dbcan_db_loc=None, dbcan_fam_activities=None,
+                           viral_db_loc=None, peptidase_db_loc=None, vogdb_db_loc=None, vog_annotations=None,
+                           description_db_loc=None, genome_summary_form_loc=None, module_step_form_loc=None,
+                           etc_module_database_loc=None, function_heatmap_form_loc=None, amg_database_loc=None,
+                           write_config=True):
+        def check_exists_and_add_to_location_dict(loc, old_value):
+            if loc is None:  # if location is none then return the old value
+                return old_value
+            elif path.isfile(loc):  # if location exists return full path
+                return path.realpath(loc)
+            else:  # if the location doesn't exist then raise error
+                raise ValueError("Database location does not exist: %s" % loc)
+        self.db_locs['kegg'] = check_exists_and_add_to_location_dict(kegg_db_loc, self.db_locs.get('kegg'))
+        self.db_locs['kofam'] = check_exists_and_add_to_location_dict(kofam_hmm_loc, self.db_locs.get('kofam'))
+        self.db_locs['kofam_ko_list'] = check_exists_and_add_to_location_dict(kofam_ko_list_loc,
+                                                                              self.db_locs.get('kofam_ko_list'))
+        self.db_locs['uniref'] = check_exists_and_add_to_location_dict(uniref_db_loc, self.db_locs.get('uniref'))
+        self.db_locs['pfam'] = check_exists_and_add_to_location_dict(pfam_db_loc, self.db_locs.get('pfam'))
+
+        self.db_locs['dbcan'] = check_exists_and_add_to_location_dict(dbcan_db_loc, self.db_locs.get('dbcan'))
+
+        self.db_locs['viral'] = check_exists_and_add_to_location_dict(viral_db_loc, self.db_locs.get('viral'))
+        self.db_locs['peptidase'] = check_exists_and_add_to_location_dict(peptidase_db_loc,
+                                                                          self.db_locs.get('peptidase'))
+        self.db_locs['vogdb'] = check_exists_and_add_to_location_dict(vogdb_db_loc, self.db_locs.get('vogdb'))
+
+        self.db_description_locs['pfam_hmm_dat'] = \
+            check_exists_and_add_to_location_dict(pfam_hmm_dat, self.db_description_locs.get('pfam_hmm_dat'))
+        self.db_description_locs['dbcan_fam_activities'] = \
+            check_exists_and_add_to_location_dict(dbcan_fam_activities,
+                                                  self.db_description_locs.get('dbcan_fam_activities'))
+        self.db_description_locs['vog_annotations'] = \
+            check_exists_and_add_to_location_dict(vog_annotations, self.db_description_locs.get('vog_annotations'))
+
+        self.dram_sheet_locs['genome_summary_form'] = \
+            check_exists_and_add_to_location_dict(genome_summary_form_loc,
+                                                  self.dram_sheet_locs.get('genome_summary_form'))
+        self.dram_sheet_locs['module_step_form'] = \
+            check_exists_and_add_to_location_dict(module_step_form_loc, self.dram_sheet_locs.get('module_step_form'))
+        self.dram_sheet_locs['etc_module_database'] = \
+            check_exists_and_add_to_location_dict(etc_module_database_loc,
+                                                  self.dram_sheet_locs.get('etc_module_database'))
+        self.dram_sheet_locs['function_heatmap_form'] = \
+            check_exists_and_add_to_location_dict(function_heatmap_form_loc,
+                                                  self.dram_sheet_locs.get('function_heatmap_form'))
+        self.dram_sheet_locs['amg_database'] = \
+            check_exists_and_add_to_location_dict(amg_database_loc, self.dram_sheet_locs.get('amg_database'))
+
+        self.description_loc = check_exists_and_add_to_location_dict(description_db_loc, self.description_loc)
+        self.start_db_session()
+
+        if write_config:
+            self.write_config()
+
+    def write_config(self, config_loc=None):
+        if config_loc is None:
+            config_loc = self.config_loc
+        with open(config_loc, 'w') as f:
+            config = {key: value for dict_ in (self.db_locs, self.db_description_locs, self.dram_sheet_locs)
+                      for key, value in dict_.items()}
+            config['description_db'] = self.description_loc
+            f.write(json.dumps(config))
+
+    @staticmethod
+    def make_header_dict_from_mmseqs_db(mmseqs_db):
+        mmseqs_headers_handle = open('%s_h' % mmseqs_db, 'rb')
+        mmseqs_headers = mmseqs_headers_handle.read().decode(errors='ignore')
+        mmseqs_headers = [i.strip() for i in mmseqs_headers.strip().split('\n\x00') if len(i) > 0]
+        return [{'id': i.split(' ')[0], 'description': i} for i in mmseqs_headers]
+
+    @staticmethod
+    def process_pfam_descriptions(pfam_hmm_dat):
+        if pfam_hmm_dat.endswith('.gz'):
+            f = gzip.open(pfam_hmm_dat, 'r').read().decode('utf-8')
+        else:
+            f = open(pfam_hmm_dat).read()
+        entries = f.strip().split('//')
+        description_list = list()
+        for i, entry in enumerate(entries):
+            if len(entry) > 0:
+                entry = entry.split('\n')
+                ascession = None
+                description = None
+                for line in entry:
+                    line = line.strip()
+                    if line.startswith('#=GF AC'):
+                        ascession = line.split('   ')[-1]
+                    if line.startswith('#=GF DE'):
+                        description = line.split('   ')[-1]
+                description_list.append({'id': ascession, 'description': description})
+        return description_list
+
+    @staticmethod
+    def process_dbcan_descriptions(dbcan_fam_activities):
+        f = open(dbcan_fam_activities)
+        description_list = list()
+        for line in f.readlines():
+            if not line.startswith('#') and len(line.strip()) != 0:
+                line = line.strip().split()
+                if len(line) == 1:
+                    description = line[0]
+                elif line[0] == line[1]:
+                    description = ' '.join(line[1:])
+                else:
+                    description = ' '.join(line)
+                description_list.append({'id': line[0], 'description': description.replace('\n', ' ')})
+        return description_list
+
+    @staticmethod
+    def process_vogdb_descriptions(vog_annotations):
+        annotations_table = pd.read_csv(vog_annotations, sep='\t', index_col=0)
+        annotations_list = [{'id': vog, 'description': '%s; %s' % (row['ConsensusFunctionalDescription'],
+                                                                   row['FunctionalCategory'])}
+                            for vog, row in annotations_table.iterrows()]
+        return annotations_list
+
+    # TODO: Make option to build on description database that already exists?
+    def populate_description_db(self, output_loc=None, update_config=True):
+        if self.description_loc is None and output_loc is None:  # description db location must be set somewhere
+            raise ValueError('Must provide output location if description db location is not set in configuration')
+        if output_loc is not None:  # if new description db location is set then save it there
+            self.description_loc = output_loc
+            self.start_db_session()
+        if path.exists(self.description_loc):
+            remove(self.description_loc)
+        create_description_db(self.description_loc)
+
+        # fill database
+        if self.db_locs.get('kegg') is not None:
+            self.add_descriptions_to_database(self.make_header_dict_from_mmseqs_db(self.db_locs['kegg']), 'kegg',
+                                              clear_table=True)
+        if self.db_locs.get('uniref') is not None:
+            self.add_descriptions_to_database(self.make_header_dict_from_mmseqs_db(self.db_locs['uniref']), 'uniref',
+                                              clear_table=True)
+        if self.db_description_locs.get('pfam_hmm_dat') is not None:
+            self.add_descriptions_to_database(self.process_pfam_descriptions(self.db_description_locs['pfam_hmm_dat']),
+                                              'pfam_description', clear_table=True)
+        if self.db_description_locs.get('dbcan_fam_activities') is not None:
+            self.add_descriptions_to_database(self.process_dbcan_descriptions(
+                self.db_description_locs['dbcan_fam_activities']), 'dbcan_description', clear_table=True)
+        if self.db_locs.get('viral') is not None:
+            self.add_descriptions_to_database(self.make_header_dict_from_mmseqs_db(self.db_locs['viral']),
+                                              'viral_description', clear_table=True)
+        if self.db_locs.get('peptidase') is not None:
+            self.add_descriptions_to_database(self.make_header_dict_from_mmseqs_db(self.db_locs['peptidase']),
+                                              'peptidase_description', clear_table=True)
+        if self.db_description_locs.get('vog_annotations') is not None:
+            self.add_descriptions_to_database(
+                self.process_vogdb_descriptions(self.db_description_locs['vog_annotations']), 'vogdb_description',
+                clear_table=True)
+
+        if update_config:  # if new description db is set then save it
+            self.write_config()
+
+    def print_database_locations(self):
+        # search databases
+        print('Processed search databases')
+        print('KEGG db: %s' % self.db_locs.get('kegg'))
+        print('KOfam db: %s' % self.db_locs.get('kofam'))
+        print('KOfam KO list: %s' % self.db_locs.get('kofam_ko_list'))
+        print('UniRef db: %s' % self.db_locs.get('uniref'))
+        print('Pfam db: %s' % self.db_locs.get('pfam'))
+        print('dbCAN db: %s' % self.db_locs.get('dbcan'))
+        print('RefSeq Viral db: %s' % self.db_locs.get('viral'))
+        print('MEROPS peptidase db: %s' % self.db_locs.get('peptidase'))
+        print('VOGDB db: %s' % self.db_locs.get('vogdb'))
+        print()
+        # database descriptions used during description db population
+        print('Descriptions of search database entries')
+        print('Pfam hmm dat: %s' % self.db_description_locs.get('pfam_hmm_dat'))
+        print('dbCAN family activities: %s' % self.db_description_locs.get('dbcan_fam_activities'))
+        print('VOG annotations: %s' % self.db_description_locs.get('vog_annotations'))
+        print()
+        # description database
+        print('Description db: %s' % self.description_loc)
+        print()
+        # DRAM sheets
+        print('DRAM distillation sheets')
+        print('Genome summary form: %s' % self.dram_sheet_locs.get('genome_summary_form'))
+        print('Module step form: %s' % self.dram_sheet_locs.get('module_step_form'))
+        print('ETC module database: %s' % self.dram_sheet_locs.get('etc_module_database'))
+        print('Function heatmap form: %s' % self.dram_sheet_locs.get('function_heatmap_form'))
+        print('AMG database: %s' % self.dram_sheet_locs.get('amg_database'))
+
+    def filter_db_locs(self, low_mem_mode=False, use_uniref=True, use_vogdb=True, master_list=None):
+        if master_list is None:
+            dbs_to_use = self.db_locs.keys()
+        else:
+            dbs_to_use = master_list
+        # filter out dbs for low mem mode
+        if low_mem_mode:
+            if ('kofam' not in self.db_locs) or ('kofam_ko_list' not in self.db_locs):
+                raise ValueError('To run in low memory mode KOfam must be configured for use in DRAM')
+            dbs_to_use = [i for i in dbs_to_use if i not in ('uniref', 'kegg', 'vogdb')]
+        # check on uniref status
+        if use_uniref:
+            if 'uniref' not in self.db_locs:
+                warnings.warn('Sequences will not be annoated against uniref as it is not configured for use in DRAM')
+        else:
+            dbs_to_use = [i for i in dbs_to_use if i != 'uniref']
+        # check on vogdb status
+        if use_vogdb:
+            if 'vogdb' not in self.db_locs:
+                warnings.warn('Sequences will not be annoated against VOGDB as it is not configured for use in DRAM')
+        else:
+            dbs_to_use = [i for i in dbs_to_use if i != 'vogdb']
+        self.db_locs = {key: value for key, value in self.db_locs.items() if key in dbs_to_use}
+
+    def clear_config(self):
+        self.db_locs = {}
+        self.db_description_locs = {}
+        self.dram_sheet_locs = {}
+        self.description_loc = None
+
+
+def set_database_paths(kegg_db_loc=None, kofam_hmm_loc=None, kofam_ko_list_loc=None, uniref_db_loc=None,
+                       pfam_db_loc=None, pfam_hmm_dat=None, dbcan_db_loc=None, dbcan_fam_activities=None,
+                       viral_db_loc=None, peptidase_db_loc=None, vogdb_db_loc=None, vog_annotations=None,
+                       description_db_loc=None, genome_summary_form_loc=None, module_step_form_loc=None,
+                       etc_module_database_loc=None, function_heatmap_form_loc=None, amg_database_loc=None,
+                       clear_config=False, update_description_db=False):
+    db_handler = DatabaseHandler()
+    if clear_config:
+        db_handler.clear_config()
+    db_handler.set_database_paths(kegg_db_loc, kofam_hmm_loc, kofam_ko_list_loc, uniref_db_loc,
+                                  pfam_db_loc, pfam_hmm_dat, dbcan_db_loc, dbcan_fam_activities,
+                                  viral_db_loc, peptidase_db_loc, vogdb_db_loc, vog_annotations,
+                                  description_db_loc, genome_summary_form_loc, module_step_form_loc,
+                                  etc_module_database_loc, function_heatmap_form_loc, amg_database_loc,
+                                  write_config=True)
+    if update_description_db:
+        db_handler.populate_description_db()
+
+
+def print_database_locations(config_loc=None):
+    db_handler = DatabaseHandler(config_loc)
+    db_handler.print_database_locations()
+
+
+def populate_description_db(output_loc=None, config_loc=None):
+    db_handler = DatabaseHandler(config_loc)
+    db_handler.populate_description_db(output_loc)
+
+
+def export_config(output_file=None):
+    config_loc = get_config_loc()
+    if output_file is None:
+        print(open(config_loc).read())
+    else:
+        copy2(config_loc, output_file)
+
+
+def import_config(config_loc):
+    system_config = get_config_loc()
+    copy2(config_loc, system_config)
diff --git a/mag_annotator/database_processing.py b/mag_annotator/database_processing.py
index 4fd787ab..fa58347d 100644
--- a/mag_annotator/database_processing.py
+++ b/mag_annotator/database_processing.py
@@ -1,20 +1,15 @@
-from os import path, mkdir, remove
-from shutil import copy2
+from os import path, mkdir
 from datetime import datetime
 from shutil import move, rmtree
 from glob import glob
-import json
 import gzip
 import tarfile
-import pandas as pd
 from collections import defaultdict
 from skbio import read as read_sequence
 from skbio import write as write_sequence
 
-from mag_annotator.utils import run_process, make_mmseqs_db, get_database_locs, download_file, merge_files, \
-    get_config_loc, remove_prefix
 from mag_annotator.database_handler import DatabaseHandler
-from mag_annotator.database_setup import create_description_db
+from mag_annotator.utils import run_process, make_mmseqs_db, download_file, merge_files, remove_prefix
 
 
 # TODO: check if dbcan or pfam is down, raise appropriate error
@@ -25,22 +20,6 @@ def get_iso_date():
     return datetime.today().strftime('%Y%m%d')
 
 
-def check_file_exists(db_loc):
-    if db_loc is None:
-        return True
-    elif path.isfile(db_loc):
-        return True
-    else:
-        raise ValueError("Database location does not exist: %s" % db_loc)
-
-
-def make_header_dict_from_mmseqs_db(mmseqs_db):
-    mmseqs_headers_handle = open('%s_h' % mmseqs_db, 'rb')
-    mmseqs_headers = mmseqs_headers_handle.read().decode(errors='ignore')
-    mmseqs_headers = [i.strip() for i in mmseqs_headers.strip().split('\n\x00') if len(i) > 0]
-    return [{'id': i.split(' ')[0], 'description': i} for i in mmseqs_headers]
-
-
 def generate_modified_kegg_fasta(kegg_fasta, gene_ko_link_loc=None):
     """Takes kegg fasta file and gene ko link file, adds kos not already in headers to headers"""
     genes_ko_dict = defaultdict(list)
@@ -141,28 +120,6 @@ def download_pfam_descriptions(output_dir='.', verbose=True):
     return pfam_hmm_dat
 
 
-def process_pfam_descriptions(pfam_hmm_dat):
-    if pfam_hmm_dat.endswith('.gz'):
-        f = gzip.open(pfam_hmm_dat, 'r').read().decode('utf-8')
-    else:
-        f = open(pfam_hmm_dat).read()
-    entries = f.strip().split('//')
-    description_list = list()
-    for i, entry in enumerate(entries):
-        if len(entry) > 0:
-            entry = entry.split('\n')
-            ascession = None
-            description = None
-            for line in entry:
-                line = line.strip()
-                if line.startswith('#=GF AC'):
-                    ascession = line.split('   ')[-1]
-                if line.startswith('#=GF DE'):
-                    description = line.split('   ')[-1]
-            description_list.append({'id': ascession, 'description': description})
-    return description_list
-
-
 def download_and_process_dbcan(dbcan_hmm=None, output_dir='.', dbcan_release='9', verbose=True):
     if dbcan_hmm is None:  # download database if not provided
         dbcan_hmm = path.join(output_dir, 'dbCAN-HMMdb-V%s.txt' % dbcan_release)
@@ -183,22 +140,6 @@ def download_dbcan_descriptions(output_dir='.', upload_date='07302020', verbose=
     return dbcan_fam_activities
 
 
-def process_dbcan_descriptions(dbcan_fam_activities):
-    f = open(dbcan_fam_activities)
-    description_list = list()
-    for line in f.readlines():
-        if not line.startswith('#') and len(line.strip()) != 0:
-            line = line.strip().split()
-            if len(line) == 1:
-                description = line[0]
-            elif line[0] == line[1]:
-                description = ' '.join(line[1:])
-            else:
-                description = ' '.join(line)
-            description_list.append({'id': line[0], 'description': description.replace('\n', ' ')})
-    return description_list
-
-
 def download_and_process_viral_refseq(merged_viral_faas=None, output_dir='.', viral_files=2, threads=10, verbose=True):
     """Can only download newest version"""
     # download all of the viral protein files, need to know the number of files
@@ -255,14 +196,6 @@ def download_vog_annotations(output_dir, vogdb_version='latest', verbose=True):
     return vog_annotations
 
 
-def process_vogdb_descriptions(vog_annotations):
-    annotations_table = pd.read_csv(vog_annotations, sep='\t', index_col=0)
-    annotations_list = [{'id': vog, 'description': '%s; %s' % (row['ConsensusFunctionalDescription'],
-                                                               row['FunctionalCategory'])}
-                        for vog, row in annotations_table.iterrows()]
-    return annotations_list
-
-
 def download_and_process_genome_summary_form(output_dir, branch='master', verbose=True):
     genome_summary_form = path.join(output_dir, 'genome_summary_form.%s.tsv' % get_iso_date())
     download_file('https://raw.githubusercontent.com/shafferm/DRAM/%s/data/genome_summary_form.tsv' % branch,
@@ -298,114 +231,13 @@ def download_and_process_amg_database(output_dir, branch='master', verbose=True)
     return amg_database
 
 
-def check_exists_and_add_to_location_dict(loc, name, dict_to_update):
-    if loc is not None:  # if location give and exists then add to dict, else raise ValueError
-        if check_file_exists(loc):
-            dict_to_update[name] = path.realpath(loc)
-    else:  # if location not given and is not in dict then set to none, else leave previous value
-        if name not in dict_to_update:
-            dict_to_update[name] = None
-    return dict_to_update
-
-
-def add_to_description_db(loc, name, get_description_list, db_handler):
-    if loc is not None:  # if location give and exists then add to dict, else raise ValueError
-        description_list = get_description_list(loc)
-        db_handler.add_descriptions_to_database(description_list, name, clear_table=True)
-
-
-def set_database_paths(kegg_db_loc=None, kofam_hmm_loc=None, kofam_ko_list_loc=None, uniref_db_loc=None,
-                       pfam_db_loc=None, pfam_hmm_dat=None, dbcan_db_loc=None, dbcan_fam_activities=None,
-                       viral_db_loc=None, peptidase_db_loc=None, vogdb_db_loc=None, vog_annotations=None,
-                       description_db_loc=None, genome_summary_form_loc=None, module_step_form_loc=None,
-                       etc_module_database_loc=None, function_heatmap_form_loc=None, amg_database_loc=None,
-                       start_time=None, config_loc=None, use_current_locs=True, update_description_db=False):
-    if start_time is None:
-        start_time = datetime.now()
-    print('%s: Setting database paths' % str(datetime.now() - start_time))
-    if use_current_locs:
-        db_dict = get_database_locs()
+def check_file_exists(db_loc):
+    if db_loc is None:
+        return True
+    elif path.isfile(db_loc):
+        return True
     else:
-        db_dict = {}
-
-    db_dict = check_exists_and_add_to_location_dict(kegg_db_loc, 'kegg', db_dict)
-    db_dict = check_exists_and_add_to_location_dict(kofam_hmm_loc, 'kofam', db_dict)
-    db_dict = check_exists_and_add_to_location_dict(kofam_ko_list_loc, 'kofam_ko_list', db_dict)
-    db_dict = check_exists_and_add_to_location_dict(uniref_db_loc, 'uniref', db_dict)
-    db_dict = check_exists_and_add_to_location_dict(pfam_db_loc, 'pfam', db_dict)
-    db_dict = check_exists_and_add_to_location_dict(pfam_hmm_dat, 'pfam_hmm_dat', db_dict)
-    db_dict = check_exists_and_add_to_location_dict(dbcan_db_loc, 'dbcan', db_dict)
-    db_dict = check_exists_and_add_to_location_dict(dbcan_fam_activities, 'dbcan_fam_activities', db_dict)
-    db_dict = check_exists_and_add_to_location_dict(viral_db_loc, 'viral', db_dict)
-    db_dict = check_exists_and_add_to_location_dict(peptidase_db_loc, 'peptidase', db_dict)
-    db_dict = check_exists_and_add_to_location_dict(vogdb_db_loc, 'vogdb', db_dict)
-    db_dict = check_exists_and_add_to_location_dict(vog_annotations, 'vog_annotations', db_dict)
-
-    db_dict = check_exists_and_add_to_location_dict(genome_summary_form_loc, 'genome_summary_form', db_dict)
-    db_dict = check_exists_and_add_to_location_dict(module_step_form_loc, 'module_step_form', db_dict)
-    db_dict = check_exists_and_add_to_location_dict(etc_module_database_loc, 'etc_module_database', db_dict)
-    db_dict = check_exists_and_add_to_location_dict(function_heatmap_form_loc, 'function_heatmap_form', db_dict)
-    db_dict = check_exists_and_add_to_location_dict(amg_database_loc, 'amg_database', db_dict)
-    print('%s: Database locations added to CONFIG' % str(datetime.now() - start_time))
-
-    if update_description_db:
-        if description_db_loc is None:
-            description_db_loc = db_dict['description_db']
-        populate_description_db(description_db_loc, db_dict, start_time)
-        print('%s: Database descriptions updated' % str(datetime.now() - start_time))
-    db_dict = check_exists_and_add_to_location_dict(description_db_loc, 'description_db', db_dict)
-
-    # change data paths
-    if config_loc is None:
-        config_loc = get_config_loc()
-    with open(config_loc, 'w') as f:
-        f.write(json.dumps(db_dict))
-    print('%s: Database locations set' % str(datetime.now() - start_time))
-
-
-def populate_description_db(output_loc=None, db_dict=None, start_time=None):
-    if start_time is None:
-        start_time = datetime.now()
-        print('%s: Populating description database' % str(datetime.now() - start_time))
-    # setup
-    if db_dict is None:
-        db_dict = get_database_locs()
-
-    if db_dict.get('description_db') is None and output_loc is not None:
-        db_dict['description_db'] = output_loc
-    elif db_dict.get('description_db') is None and output_loc is None:
-        raise ValueError('Must provide output location if description db location is not set in configuration')
-    elif path.exists(db_dict['description_db']):
-        remove(db_dict['description_db'])
-
-    create_description_db(db_dict['description_db'])
-    db_handler = DatabaseHandler(db_dict['description_db'])
-    print('%s: Database connection established' % str(datetime.now() - start_time))
-
-    # fill database
-    add_to_description_db(db_dict['kegg'], 'kegg_description', make_header_dict_from_mmseqs_db,
-                          db_handler)
-    print('%s: KEGG descriptions added to description database' % str(datetime.now() - start_time))
-    add_to_description_db(db_dict['uniref'], 'uniref_description', make_header_dict_from_mmseqs_db,
-                          db_handler)
-    print('%s: UniRef descriptions added to description database' % str(datetime.now() - start_time))
-    add_to_description_db(db_dict['pfam_hmm_dat'], 'pfam_description', process_pfam_descriptions,
-                          db_handler)
-    print('%s: PFAM descriptions added to description database' % str(datetime.now() - start_time))
-    add_to_description_db(db_dict['dbcan_fam_activities'], 'dbcan_description',
-                          process_dbcan_descriptions, db_handler)
-    print('%s: dbCAN descriptions added to description database' % str(datetime.now() - start_time))
-    add_to_description_db(db_dict['viral'], 'viral_description', make_header_dict_from_mmseqs_db,
-                          db_handler)
-    print('%s: RefSeq viral descriptions added to description database' % str(datetime.now() - start_time))
-    add_to_description_db(db_dict['peptidase'], 'peptidase_description',
-                          make_header_dict_from_mmseqs_db, db_handler)
-    print('%s: MEROPS descriptions added to description database' % str(datetime.now() - start_time))
-    add_to_description_db(db_dict['vog_annotations'], 'vogdb_description', process_vogdb_descriptions,
-                          db_handler)
-    print('%s: VOGdb descriptions added to description database' % str(datetime.now() - start_time))
-    print('%s: Description database populated' % str(datetime.now() - start_time))
-    # TODO: if output loc is given then set in config file where the description database is
+        raise ValueError("Database location does not exist: %s" % db_loc)
 
 
 def prepare_databases(output_dir, kegg_loc=None, gene_ko_link_loc=None, kofam_hmm_loc=None, kofam_ko_list_loc=None,
@@ -523,37 +355,16 @@ def prepare_databases(output_dir, kegg_loc=None, gene_ko_link_loc=None, kofam_hm
 
     output_dbs['description_db_loc'] = path.realpath(path.join(output_dir, 'description_db.sqlite'))
 
-    set_database_paths(**output_dbs, use_current_locs=False, update_description_db=True, start_time=start_time)
+    db_handler = DatabaseHandler()
+    db_handler.populate_description_db(output_dbs['description_db_loc'], update_config=False)
+    db_handler.set_database_paths(**output_dbs)
+    print('%s: DRAM description database populated' % str(datetime.now() - start_time))
 
     if not keep_database_files:
         rmtree(temporary)
     print('%s: Database preparation completed' % str(datetime.now() - start_time))
 
 
-def print_database_locations(db_locs=None):
-    if db_locs is None:
-        db_locs = get_database_locs()
-
-    print('KEGG db: %s' % db_locs.get('kegg'))
-    print('KOfam db: %s' % db_locs.get('kofam'))
-    print('KOfam KO list: %s' % db_locs.get('kofam_ko_list'))
-    print('UniRef db: %s' % db_locs.get('uniref'))
-    print('Pfam db: %s' % db_locs.get('pfam'))
-    print('Pfam hmm dat: %s' % db_locs.get('pfam_hmm_dat'))
-    print('dbCAN db: %s' % db_locs.get('dbcan'))
-    print('dbCAN family activities: %s' % db_locs.get('dbcan_fam_activities'))
-    print('RefSeq Viral db: %s' % db_locs.get('viral'))
-    print('MEROPS peptidase db: %s' % db_locs.get('peptidase'))
-    print('VOGDB db: %s' % db_locs.get('vogdb'))
-    print('VOG annotations: %s' % db_locs.get('vog_annotations'))
-    print('Description db: %s' % db_locs.get('description_db'))
-    print('Genome summary form: %s' % db_locs.get('genome_summary_form'))
-    print('Module step form: %s' % db_locs.get('module_step_form'))
-    print('ETC module database: %s' % db_locs.get('etc_module_database'))
-    print('Function heatmap form: %s' % db_locs.get('function_heatmap_form'))
-    print('AMG database: %s' % db_locs.get('amg_database'))
-
-
 def update_dram_forms(output_dir, branch='master'):
     if not path.isdir(output_dir):
         mkdir(output_dir)
@@ -564,15 +375,5 @@ def update_dram_forms(output_dir, branch='master'):
     form_locs['etc_module_database_loc'] = download_and_process_etc_module_database(output_dir, branch)
     form_locs['function_heatmap_form_loc'] = download_and_process_function_heatmap_form(output_dir, branch)
     form_locs['amg_database_loc'] = download_and_process_amg_database(output_dir, branch)
-    set_database_paths(**form_locs, update_description_db=False)
-
-
-def export_config(output_file=None):
-    if output_file is None:
-        print(open(get_config_loc()).read())
-    else:
-        copy2(get_config_loc(), output_file)
-
-
-def import_config(config_loc):
-    copy2(config_loc, get_config_loc())
+    db_handler = DatabaseHandler()
+    db_handler.set_database_paths(**form_locs)
diff --git a/mag_annotator/pull_sequences.py b/mag_annotator/pull_sequences.py
index cbdccd16..80166502 100644
--- a/mag_annotator/pull_sequences.py
+++ b/mag_annotator/pull_sequences.py
@@ -4,27 +4,80 @@
 from os import mkdir, path
 import warnings
 
-from mag_annotator.utils import get_genes_from_identifiers
 from mag_annotator.summarize_vgfs import filter_to_amgs
+from mag_annotator.utils import get_ids_from_row
+from mag_annotator.database_handler import DatabaseHandler
 
 # TODO: filter by taxonomic level, completeness, contamination
 # TODO: filter scaffolds file, gff file
 # TODO: add negate, aka pull not from given list
 
 
-# def pull_sequences(input_annotations, input_fasta, output_fasta, fastas=None, scaffolds=None, genes=None,
-#                    identifiers=None, categories=None, taxonomy=None, completeness=None, contamination=None,
-#                    amg_flags=None, aux_scores=None, virsorter_category=None, putative_amgs=False, max_auxiliary_score=3,
-#                    remove_transposons=False, remove_fs=False, remove_js=False):
+def get_genes_from_identifiers(annotations, genes=None, fastas=None, scaffolds=None, identifiers=None, categories=None,
+                               custom_distillate=None):
+    specific_genes_to_keep = list()
+    # filter fastas
+    if fastas is not None:
+        for fasta in fastas:
+            specific_genes_to_keep += list(annotations.loc[annotations['fasta'] == fasta].index)
+    # filter scaffolds
+    # TODO: remove this functionality or modify since scaffolds are guaranteed unique
+    if scaffolds is not None:
+        for scaffold in scaffolds:
+            specific_genes_to_keep += list(annotations.loc[annotations['scaffold'] == scaffold].index)
+    # filter genes
+    if genes is not None:
+        specific_genes_to_keep += genes
+    # filter down annotations based on specific genes
+    if len(specific_genes_to_keep) > 0:
+        annotations_to_keep = annotations.loc[specific_genes_to_keep]
+    else:
+        annotations_to_keep = annotations
+
+    # filter based on annotations
+    if (identifiers is not None) or (categories is not None):
+        annotation_genes_to_keep = list()
+
+        # make a dictionary of genes to
+        gene_to_ids = dict()
+        for i, row in annotations_to_keep.iterrows():
+            row_ids = get_ids_from_row(row)
+            if len(row_ids) > 0:
+                gene_to_ids[i] = set(row_ids)
+
+        # get genes with ids
+        if identifiers is not None:
+            identifiers = set(identifiers)
+            for gene, ids in gene_to_ids.items():
+                if len(set(ids) & set(identifiers)) > 0:
+                    annotation_genes_to_keep.append(gene)
+
+        # get genes from distillate categories
+        if categories is not None:
+            database_handler = DatabaseHandler()
+            genome_summary_form = pd.read_csv(database_handler.db_locs['genome_summary_form'], sep='\t')
+            if custom_distillate is not None:
+                genome_summary_form = pd.concat([genome_summary_form, pd.read_csv(custom_distillate, sep='\t')])
+            for level in ['module', 'sheet', 'header', 'subheader']:
+                for category, frame in genome_summary_form.loc[~pd.isna(genome_summary_form[level])].groupby(level):
+                    if category in categories:
+                        for gene, ids in gene_to_ids.items():
+                            if len(ids & set(frame['gene_id'])) > 0:
+                                annotation_genes_to_keep.append(gene)
+    else:
+        annotation_genes_to_keep = list(annotations_to_keep.index)
+    return annotation_genes_to_keep
+
+
 def pull_sequences(input_annotations, input_fasta, output_fasta, fastas=None, scaffolds=None, genes=None,
                    identifiers=None, categories=None, taxonomy=None, completeness=None, contamination=None,
                    amg_flags=None, aux_scores=None, virsorter_category=None, putative_amgs=False,
-                   max_auxiliary_score=3, remove_transposons=False, remove_fs=False):
+                   max_auxiliary_score=3, remove_transposons=False, remove_fs=False, custom_distillate=None):
     annotations = pd.read_csv(input_annotations, sep='\t', index_col=0)
 
     # first filter based on specific names
     annotation_genes_to_keep = get_genes_from_identifiers(annotations, genes, fastas, scaffolds, identifiers,
-                                                          categories)
+                                                          categories, custom_distillate)
     annotations = annotations.loc[annotation_genes_to_keep]
     if len(annotations) == 0:
         raise ValueError("Categories or identifiers provided yielded no annotations")
@@ -99,7 +152,7 @@ def find_neighborhoods(annotations, genes_from_ids, distance_bp=None, distance_g
 
 
 def get_gene_neighborhoods(input_file, output_dir, genes=None, identifiers=None, categories=None, genes_loc=None,
-                           scaffolds_loc=None, distance_genes=None, distance_bp=None):
+                           scaffolds_loc=None, distance_genes=None, distance_bp=None, custom_distillate=None):
     # check inputs, make output
     if distance_genes is None and distance_bp is None:
         raise ValueError("Must provide distance away in bp, genes or both.")
@@ -107,7 +160,7 @@ def get_gene_neighborhoods(input_file, output_dir, genes=None, identifiers=None,
     # get data
     annotations = pd.read_csv(input_file, sep='\t', index_col=0)
     genes_from_ids = get_genes_from_identifiers(annotations, genes=genes, identifiers=identifiers,
-                                                categories=categories)
+                                                categories=categories, custom_distillate=custom_distillate)
     if len(genes_from_ids) == 0:
         raise ValueError("No genes were found based on your filtering parameters. No neighborhoods will be generated.")
 
diff --git a/mag_annotator/summarize_genomes.py b/mag_annotator/summarize_genomes.py
index 75bdac8e..3de203bb 100644
--- a/mag_annotator/summarize_genomes.py
+++ b/mag_annotator/summarize_genomes.py
@@ -8,7 +8,8 @@
 import numpy as np
 from datetime import datetime
 
-from mag_annotator.utils import get_database_locs, get_ids_from_annotation, get_ids_from_row
+from mag_annotator.database_handler import DatabaseHandler
+from mag_annotator.utils import get_ids_from_annotation, get_ids_from_row, get_ordered_uniques
 
 # TODO: add RBH information to output
 # TODO: add flag to output table and not xlsx
@@ -26,12 +27,6 @@
 TAXONOMY_LEVELS = ['d', 'p', 'c', 'o', 'f', 'g', 's']
 
 
-def get_ordered_uniques(seq):
-    seen = set()
-    seen_add = seen.add
-    return [x for x in seq if not (x in seen or seen_add(x))]
-
-
 def fill_genome_summary_frame(annotations, genome_summary_frame, groupby_column):
     genome_summary_id_sets = [set([k.strip() for k in j.split(',')]) for j in genome_summary_frame['gene_id']]
     for genome, frame in annotations.groupby(groupby_column, sort=False):
@@ -574,22 +569,22 @@ def summarize_genomes(input_file, trna_path=None, rrna_path=None, output_dir='.'
         rrna_frame = pd.read_csv(rrna_path, sep='\t')
 
     # get db_locs and read in dbs
-    db_locs = get_database_locs()
-    if 'genome_summary_form' not in db_locs:
+    database_handler = DatabaseHandler()
+    if 'genome_summary_form' not in database_handler.dram_sheet_locs:
         raise ValueError('Genome summary form location must be set in order to summarize genomes')
-    if 'module_step_form' not in db_locs:
+    if 'module_step_form' not in database_handler.dram_sheet_locs:
         raise ValueError('Module step form location must be set in order to summarize genomes')
-    if 'function_heatmap_form' not in db_locs:
+    if 'function_heatmap_form' not in database_handler.dram_sheet_locs:
         raise ValueError('Functional heat map form location must be set in order to summarize genomes')
 
     # read in dbs
-    genome_summary_form = pd.read_csv(db_locs['genome_summary_form'], sep='\t')
+    genome_summary_form = pd.read_csv(database_handler.dram_sheet_locs['genome_summary_form'], sep='\t')
     if custom_distillate is not None:
         genome_summary_form = pd.concat([genome_summary_form, pd.read_csv(custom_distillate, sep='\t')])
     genome_summary_form = genome_summary_form.drop('potential_amg', axis=1)
-    module_steps_form = pd.read_csv(db_locs['module_step_form'], sep='\t')
-    function_heatmap_form = pd.read_csv(db_locs['function_heatmap_form'], sep='\t')
-    etc_module_df = pd.read_csv(db_locs['etc_module_database'], sep='\t')
+    module_steps_form = pd.read_csv(database_handler.dram_sheet_locs['module_step_form'], sep='\t')
+    function_heatmap_form = pd.read_csv(database_handler.dram_sheet_locs['function_heatmap_form'], sep='\t')
+    etc_module_df = pd.read_csv(database_handler.dram_sheet_locs['etc_module_database'], sep='\t')
     print('%s: Retrieved database locations and descriptions' % (str(datetime.now() - start_time)))
 
     # make output folder
diff --git a/mag_annotator/summarize_vgfs.py b/mag_annotator/summarize_vgfs.py
index 8b393267..07403ddb 100644
--- a/mag_annotator/summarize_vgfs.py
+++ b/mag_annotator/summarize_vgfs.py
@@ -7,8 +7,8 @@
 from datetime import datetime
 import warnings
 
-from mag_annotator.utils import get_database_locs, get_ids_from_row, get_ids_from_annotation
-from mag_annotator.summarize_genomes import get_ordered_uniques
+from mag_annotator.database_handler import DatabaseHandler
+from mag_annotator.utils import get_ids_from_row, get_ids_from_annotation, get_ordered_uniques
 
 VOGDB_TYPE_NAMES = {'Xr': 'Viral replication genes', 'Xs': 'Viral structure genes',
                     'Xh': 'Viral genes with host benefits', 'Xp': 'Viral genes with viral benefits',
@@ -217,11 +217,11 @@ def summarize_vgfs(input_file, output_dir, groupby_column='scaffold', max_auxili
 
     # set up
     annotations = pd.read_csv(input_file, sep='\t', index_col=0).fillna('')
-    db_locs = get_database_locs()
-    if 'genome_summary_form' not in db_locs:
+    database_handler = DatabaseHandler()
+    if database_handler.dram_sheet_locs.get('genome_summary_form') is None:
         raise ValueError('Genome summary form location must be set in order to summarize genomes')
     mkdir(output_dir)
-    genome_summary_form = pd.read_csv(db_locs['genome_summary_form'], sep='\t', index_col=0)
+    genome_summary_form = pd.read_csv(database_handler.dram_sheet_locs['genome_summary_form'], sep='\t', index_col=0)
     if custom_distillate is not None:
         custom_distillate_form = pd.read_csv(custom_distillate, sep='\t', index_col=0)
         genome_summary_form = pd.concat([genome_summary_form, custom_distillate_form])
diff --git a/mag_annotator/utils.py b/mag_annotator/utils.py
index 20f614bc..80633343 100644
--- a/mag_annotator/utils.py
+++ b/mag_annotator/utils.py
@@ -2,22 +2,10 @@
 import subprocess
 from collections import Counter
 from os import path
-from pkg_resources import resource_filename
-import json
 from urllib.request import urlopen
 import pandas as pd
 
 
-def get_config_loc():
-    return path.abspath(resource_filename('mag_annotator', 'CONFIG'))
-
-
-def get_database_locs(config_loc=None):
-    if config_loc is None:
-        config_loc = get_config_loc()
-    return json.loads(open(config_loc).read())
-
-
 def download_file(url, output_file=None, verbose=True):
     # TODO: catching error 4 and give error message to retry or retry automatically
     if verbose:
@@ -91,11 +79,12 @@ def get_ids_from_annotation(frame):
     if 'peptidase_family' in frame:
         id_list += [j.strip() for i in frame.peptidase_family.dropna() for j in i.split(';')]
     # get cazy ids
+    if 'cazy_id' in frame:
+        id_list += [j for i in frame.cazy_id.dropna() for j in i.split('; ')]
+    # get cazy ec numbers
     if 'cazy_hits' in frame:
-        id_list += [j.split(' ')[0] for i in frame.cazy_hits.dropna() for j in i.split(';')]
-        # get cazy ec numbers
-        for cazy_hit in frame.cazy_hits.dropna():
-            id_list += [i[1:-1].split('_')[0] for i in re.findall(r'\[[A-Z]*\d*?\]', cazy_hit)]
+        id_list += [f"{j[1:3]}:{j[4:-1]}" for i in frame.cazy_hits.dropna()
+                    for j in re.findall(r'\(EC [\d+\.]+[\d-]\)', i)]
     # get pfam ids
     if 'pfam_hits' in frame:
         id_list += [j[1:-1].split('.')[0] for i in frame.pfam_hits.dropna()
@@ -103,7 +92,7 @@ def get_ids_from_annotation(frame):
     return Counter(id_list)
 
 
-# unify this with get_ids_from_annotation
+#TODO unify this with get_ids_from_annotation
 def get_ids_from_row(row):
     id_list = list()
     # get kegg gene ids
@@ -146,54 +135,7 @@ def remove_suffix(text, suffix):
     return text  # or whatever
 
 
-def get_genes_from_identifiers(annotations, genes=None, fastas=None, scaffolds=None, identifiers=None, categories=None):
-    specific_genes_to_keep = list()
-    # filter fastas
-    if fastas is not None:
-        for fasta in fastas:
-            specific_genes_to_keep += list(annotations.loc[annotations['fasta'] == fasta].index)
-    # filter scaffolds
-    # TODO: remove this functionality or modify since scaffolds are guaranteed unique
-    if scaffolds is not None:
-        for scaffold in scaffolds:
-            specific_genes_to_keep += list(annotations.loc[annotations['scaffold'] == scaffold].index)
-    # filter genes
-    if genes is not None:
-        specific_genes_to_keep += genes
-    # filter down annotations based on specific genes
-    if len(specific_genes_to_keep) > 0:
-        annotations_to_keep = annotations.loc[specific_genes_to_keep]
-    else:
-        annotations_to_keep = annotations
-
-    # filter based on annotations
-    if (identifiers is not None) or (categories is not None):
-        annotation_genes_to_keep = list()
-
-        # make a dictionary of genes to
-        gene_to_ids = dict()
-        for i, row in annotations_to_keep.iterrows():
-            row_ids = get_ids_from_row(row)
-            if len(row_ids) > 0:
-                gene_to_ids[i] = set(row_ids)
-
-        # get genes with ids
-        if identifiers is not None:
-            identifiers = set(identifiers)
-            for gene, ids in gene_to_ids.items():
-                if len(set(ids) & set(identifiers)) > 0:
-                    annotation_genes_to_keep.append(gene)
-
-        # get genes from distillate categories
-        if categories is not None:
-            db_locs = get_database_locs()
-            genome_summary_form = pd.read_csv(db_locs['genome_summary_form'], sep='\t')
-            for level in ['module', 'sheet', 'header', 'subheader']:
-                for category, frame in genome_summary_form.loc[~pd.isna(genome_summary_form[level])].groupby(level):
-                    if category in categories:
-                        for gene, ids in gene_to_ids.items():
-                            if len(ids & set(frame['gene_id'])) > 0:
-                                annotation_genes_to_keep.append(gene)
-    else:
-        annotation_genes_to_keep = list(annotations_to_keep.index)
-    return annotation_genes_to_keep
+def get_ordered_uniques(seq):
+    seen = set()
+    seen_add = seen.add
+    return [x for x in seq if not (x in seen or seen_add(x))]
diff --git a/scripts/DRAM-setup.py b/scripts/DRAM-setup.py
index 625c47dc..ca7ee567 100644
--- a/scripts/DRAM-setup.py
+++ b/scripts/DRAM-setup.py
@@ -2,8 +2,9 @@
 
 import argparse
 
-from mag_annotator.database_processing import prepare_databases, set_database_paths, print_database_locations,\
-                                              populate_description_db, update_dram_forms, export_config, import_config
+from mag_annotator.database_processing import prepare_databases, update_dram_forms
+from mag_annotator.database_handler import set_database_paths, print_database_locations, populate_description_db, \
+    export_config, import_config
 from mag_annotator import __version__ as version
 
 
@@ -123,6 +124,9 @@ def print_version():
     # parser for updating database descriptions
     update_description_db_parser.add_argument('--output_loc', help="Location to store desciption database, will be "
                                                                    "stored in location set in CONFIG if not given")
+    update_description_db_parser.add_argument('--config_loc', help="Location of CONFIG file to use for finding "
+                                                                   "databases, by default the location in the built in "
+                                                                   "CONFIG will be used")
     update_description_db_parser.set_defaults(func=populate_description_db)
 
     # parser for updating DRAM databases only
@@ -132,6 +136,9 @@ def print_version():
     update_dram_forms_parser.set_defaults(func=update_dram_forms)
 
     # parser for printing out database configuration information
+    print_db_locs_parser.add_argument('--config_loc', help="Location of CONFIG to print locations from, by "
+                                                                   "default the locations from the built in CONFIG "
+                                                                   "will be used")
     print_db_locs_parser.set_defaults(func=print_database_locations)
 
     # parser for printing out or saving CONFIG to file
@@ -140,7 +147,9 @@ def print_version():
     export_config_parser.set_defaults(func=export_config)
 
     # parser for importing CONFIG file
-    import_config_parser.add_argument('--config_loc', help="CONFIG file to replace current CONFIG with")
+    import_config_parser.add_argument('--config_loc', help="CONFIG file to replace current CONFIG with. This will "
+                                                           "overwrite your previous configuration, export it if you "
+                                                           "would like to save it for future use.")
     import_config_parser.set_defaults(func=import_config)
 
     args = parser.parse_args()
diff --git a/scripts/DRAM-v.py b/scripts/DRAM-v.py
index 02575348..cba67f98 100644
--- a/scripts/DRAM-v.py
+++ b/scripts/DRAM-v.py
@@ -48,11 +48,19 @@
     annotate_parser.add_argument('--kofam_use_dbcan2_thresholds', action='store_true', default=False,
                                  help='Use dbcan2 suggested HMM cutoffs for KOfam annotation instead of KOfam '
                                       'recommended cutoffs. This will be ignored if annotating with KEGG Genes.')
-    annotate_parser.add_argument('--custom_db_name', action='append', help="Names of custom databases, can be used"
-                                                                           "multiple times.")
-    annotate_parser.add_argument('--custom_fasta_loc', action='append',
-                                 help="Location of fastas to annotated against, can be used multiple times but"
+    annotate_parser.add_argument('--custom_db_name', action='append', default=(),
+                                 help="Names of custom databases, can be used multiple times.")
+    annotate_parser.add_argument('--custom_fasta_loc', action='append', default=(),
+                                 help="Location of fastas to annotate against, can be used multiple times but"
                                       "must match nubmer of custom_db_name's")
+    annotate_parser.add_argument('--custom_hmm_name', action='append',  default=(),
+                                 help="Names of custom hmm databases, can be used multiple times.")
+    annotate_parser.add_argument('--custom_hmm_loc', action='append',  default=(),
+                                 help="Location of hmms to annotate against, can be used multiple times but"
+                                      "must match nubmer of custom_hmm_name's")
+    annotate_parser.add_argument('--custom_hmm_cutoffs_loc', action='append',
+                                       help="Location of file with custom HMM cutoffs and descriptions, can be used "
+                                            "multiple times.")
     annotate_parser.add_argument('--use_uniref', action='store_true', default=False,
                                  help='Annotate these fastas against UniRef, drastically increases run time and memory '
                                       'requirements')
diff --git a/scripts/DRAM.py b/scripts/DRAM.py
index 6bccba27..82d5353c 100644
--- a/scripts/DRAM.py
+++ b/scripts/DRAM.py
@@ -53,13 +53,21 @@
                                       'recommended cutoffs. This will be ignored if annotating with KEGG Genes.')
     annotate_parser.add_argument('--custom_db_name', action='append', help="Names of custom databases, can be used"
                                                                            "multiple times.")
-    annotate_parser.add_argument('--custom_fasta_loc', action='append',
-                                 help="Location of fastas to annotated against, can be used multiple times but"
+    annotate_parser.add_argument('--custom_fasta_loc', action='append', default=(),
+                                 help="Location of fastas to annotate against, can be used multiple times but"
                                       "must match nubmer of custom_db_name's")
-    annotate_parser.add_argument('--gtdb_taxonomy', action='append', default=[],
+    annotate_parser.add_argument('--custom_hmm_name', action='append',  default=(),
+                                 help="Names of custom hmm databases, can be used multiple times.")
+    annotate_parser.add_argument('--custom_hmm_loc', action='append', default=(),
+                                 help="Location of hmms to annotate against, can be used multiple times but"
+                                      "must match nubmer of custom_hmm_name's")
+    annotate_parser.add_argument('--custom_hmm_cutoffs_loc', action='append', default=(),
+                                 help="Location of file with custom HMM cutoffs and descriptions, can be used "
+                                      "multiple times.")
+    annotate_parser.add_argument('--gtdb_taxonomy', action='append', default=(),
                                  help='Summary file from gtdbtk taxonomy assignment from bins, can be used multiple'
                                       'times')
-    annotate_parser.add_argument('--checkm_quality', action='append', default=[],
+    annotate_parser.add_argument('--checkm_quality', action='append', default=(),
                                  help='Summary of of checkM quality assessment from bins, can be used multiple times')
     annotate_parser.add_argument('--use_uniref', action='store_true', default=False,
                                  help='Annotate these fastas against UniRef, drastically increases run time and memory '
@@ -86,11 +94,19 @@
     annotate_genes_parser.add_argument('--kofam_use_dbcan2_thresholds', action='store_true', default=False,
                                        help='Use dbcan2 suggested HMM cutoffs for KOfam annotation instead of KOfam '
                                             'recommended cutoffs. This will be ignored if annotating with KEGG Genes.')
-    annotate_genes_parser.add_argument('--custom_db_name', action='append', help="Names of custom databases, can be "
-                                                                                 "used multiple times.")
-    annotate_genes_parser.add_argument('--custom_fasta_loc', action='append',
-                                       help="Location of fastas to annotated against, can be used multiple times but"
+    annotate_genes_parser.add_argument('--custom_db_name', action='append', default=(),
+                                       help="Names of custom databases, can be used multiple times.")
+    annotate_genes_parser.add_argument('--custom_fasta_loc', action='append', default=(),
+                                       help="Location of fastas to annotate against, can be used multiple times but"
                                             "must match nubmer of custom_db_name's")
+    annotate_genes_parser.add_argument('--custom_hmm_name', action='append', default=(),
+                                       help="Names of custom hmm databases, can be used multiple times.")
+    annotate_genes_parser.add_argument('--custom_hmm_loc', action='append', default=(),
+                                       help="Location of hmms to annotate against, can be used multiple times but"
+                                            "must match nubmer of custom_hmm_name's")
+    annotate_genes_parser.add_argument('--custom_hmm_cutoffs_loc', action='append', default=(),
+                                       help="Location of file with custom HMM cutoffs and descriptions, can be used "
+                                            "multiple times.")
     annotate_genes_parser.add_argument('--use_uniref', action='store_true', default=False,
                                        help='Annotate these fastas against UniRef, drastically increases run time and '
                                             'memory requirements')
@@ -133,6 +149,7 @@
     annotation_group.add_argument('--identifiers', nargs='*', default=None, help='database identifiers to keep')
     annotation_group.add_argument('--categories', nargs='*', default=None,
                                   help='distillate categories to keep genes from')
+    annotation_group.add_argument('--custom_distillate', help="Custom distillate form to add your own modules")
     dram_group = strainer_parser.add_argument_group('DRAM based filters')
     dram_group.add_argument('--taxonomy', nargs='*', default=None,
                             help='Level of GTDBTk taxonomy to keep (e.g. c__Clostridia), space separated list')
@@ -148,6 +165,7 @@
     neighborhood_parser.add_argument("--genes", nargs='*', help="Gene names from DRAM to find neighborhoods around")
     neighborhood_parser.add_argument("--identifiers", nargs='*',
                                      help="Database identifiers assigned by DRAM to find neighborhoods around")
+    neighborhood_parser.add_argument('--custom_distillate', help="Custom distillate form to add your own modules")
     neighborhood_parser.add_argument("--categories", help="Distillate categories to build gene neighborhoods around.")
     neighborhood_parser.add_argument("--genes_loc", help="Location of genes.fna/genes.faa file to filter to "
                                                          "neighborhoods")


Test patch:
diff --git a/tests/data/hmm_thresholds.txt b/tests/data/hmm_thresholds.txt
new file mode 100644
index 00000000..4083eaf6
--- /dev/null
+++ b/tests/data/hmm_thresholds.txt
@@ -0,0 +1,10 @@
+hmm_name	threshold	score_type	definition
+K00001	550	full	KO1; description 1
+K00002	450	full	KO2; description 2
+K00003	450	full	DZR; description 3
+K00004	450	full	DZR; description 4
+K00005	450	full	DZR; description 5
+K00006	450	full	DZR; description 6
+K00007	450	full	DZR; description 7
+K00008	450	full	DZR; description 8
+K00009	450	full	DZR; description 9
diff --git a/tests/data/hmmsearch_hit.txt b/tests/data/hmmsearch_hit.txt
index 3c13de02..82874540 100644
--- a/tests/data/hmmsearch_hit.txt
+++ b/tests/data/hmmsearch_hit.txt
@@ -8,7 +8,7 @@ NP_040710.1          -             38 Microvir_J           PF04726.13    37   6.
 # Pipeline mode:   SEARCH
 # Query file:      Microvir_J.hmm
 # Target file:     NC_001422.faa
-# Option settings: hmmsearch --domtblout hits.txt Microvir_J.hmm NC_001422.faa 
+# Option settings: hmmsearch --domtblout hits.txt Microvir_J.hmm NC_001422.faa
 # Current dir:     /Users/shafferm/lab/DRAM/tests/data
 # Date:            Tue Mar  3 11:53:34 2020
 # [ok]
diff --git a/tests/data/test_CONFIG b/tests/data/test_CONFIG
new file mode 100644
index 00000000..3ab8cac2
--- /dev/null
+++ b/tests/data/test_CONFIG
@@ -0,0 +1,1 @@
+{"kegg": null, "kofam": null, "kofam_ko_list": null, "uniref": null, "pfam": "tests/data/Pfam-A_subset.full.gz", "pfam_hmm_dat": "tests/data/Pfam-A_subset.hmm.dat.gz", "dbcan": null, "dbcan_fam_activities": "tests/data/CAZyDB.07312019.fam-activities.subset.txt", "viral": null, "peptidase": null, "vogdb": null, "vog_annotations": "tests/data/vog_annotations_latest.subset.tsv.gz", "genome_summary_form": "data/genome_summary_form.tsv", "module_step_form": "data/module_step_form.tsv", "etc_module_database": "data/etc_module_database.tsv", "function_heatmap_form": "data/function_heatmap_form.tsv", "amg_database": "data/amg_database.tsv"}
\ No newline at end of file
diff --git a/tests/data/unformatted_cazy.b6 b/tests/data/unformatted_cazy.b6
new file mode 100644
index 00000000..d4ec5773
--- /dev/null
+++ b/tests/data/unformatted_cazy.b6
@@ -0,0 +1,26 @@
+#                                                                            --- full sequence --- -------------- this domain -------------   hmm coord   ali coord   env coord
+# target name        accession   tlen query name           accession   qlen   E-value  score  bias   #  of  c-Evalue  i-Evalue  score  bias  from    to  from    to  from    to  acc description of target
+#------------------- ---------- ----- -------------------- ---------- ----- --------- ------ ----- --- --- --------- --------- ------ ----- ----- ----- ----- ----- ----- ----- ---- ---------------------
+bin_1.scaffold_2 -            463 GH1.hmm              -            429    6e-115  377.0  10.2   1   1  7.3e-117  7.3e-115  376.7  10.2     5   427     6   459     3   461 0.92 rank: A; erb:A4V01_06220 glycosyl hydrolase family 1 (db=kegg)
+bin_1.scaffold_2 -            463 GH5.hmm              -            275     6e-08   25.2   7.3   1   1   4.1e-09   4.1e-07   22.5   7.3    34   261    45   443    14   456 0.57 rank: A; erb:A4V01_06220 glycosyl hydrolase family 1 (db=kegg)
+bin_1.scaffold_1 -            478 GT35.hmm             -            674   1.1e-07   23.4   0.0   1   1   2.5e-09   2.5e-07   22.1   0.0   377   527   251   385   245   476 0.76 rank: A; erb:A4V01_04670 starch synthase; K00703 (db=kegg)
+bin_1.scaffold_1 -            478 GT4.hmm              -            160   1.6e-16   53.1   0.0   1   1   3.1e-18   3.1e-16   52.2   0.0    11   142   296   431   288   440 0.87 rank: A; erb:A4V01_04670 starch synthase; K00703 (db=kegg)
+bin_1.scaffold_1 -            478 GT5.hmm              -            472  1.2e-182  600.3   0.6   1   1  1.4e-184  1.4e-182  600.2   0.6     1   471     2   474     2   475 0.98 rank: A; erb:A4V01_04670 starch synthase; K00703 (db=kegg)
+bin_1.other -            155 CBM76.hmm            -            161   0.00079   11.9   0.0   1   1   1.1e-05    0.0011   11.5   0.0    81   153    21    91    11    97 0.82 rank: C; erb:A4V01_09140 hypothetical protein (db=kegg)
+bin_1.scaffold_2 -            463 GH5_19.hmm           -            273   2.6e-05   16.4   0.1   1   1   4.3e-07   4.3e-05   15.6   0.1    24   121    61   155    42   160 0.79 rank: A; erb:A4V01_06220 glycosyl hydrolase family 1 (db=kegg)
+bin_1.scaffold_2 -            463 GH5_36.hmm           -            287   0.00029   12.9   0.6   1   1   4.6e-06   0.00046   12.2   0.6    18   107    61   147    54   208 0.79 rank: A; erb:A4V01_06220 glycosyl hydrolase family 1 (db=kegg)
+bin_1.scaffold_2 -            463 GH5_39.hmm           -            293   0.00054   12.1   0.0   1   1   8.7e-06   0.00087   11.4   0.0     9   105    61   158    57   181 0.77 rank: A; erb:A4V01_06220 glycosyl hydrolase family 1 (db=kegg)
+bin_1.scaffold_2 -            463 GH5_46.hmm           -            257   7.1e-06   18.2   0.0   1   1   1.6e-07   1.6e-05   17.1   0.0    17    66    60   111    53   120 0.80 rank: A; erb:A4V01_06220 glycosyl hydrolase family 1 (db=kegg)
+bin_1.scaffold_2 -            463 GH5_4.hmm            -            276   0.00024   13.5   0.7   1   2   7.9e-06   0.00079   11.8   0.0    13    74    62   124    58   153 0.85 rank: A; erb:A4V01_06220 glycosyl hydrolase family 1 (db=kegg)
+bin_1.scaffold_2 -            463 GH5_4.hmm            -            276   0.00024   13.5   0.7   2   2     0.059       5.9   -0.9   0.1   190   269   321   409   305   416 0.67 rank: A; erb:A4V01_06220 glycosyl hydrolase family 1 (db=kegg)
+bin_1.scaffold_2 -            463 GH5_5.hmm            -            282   0.00019   13.3   0.0   1   1   3.4e-06   0.00034   12.6   0.0    24    80    56   113    45   121 0.83 rank: A; erb:A4V01_06220 glycosyl hydrolase family 1 (db=kegg)
+#
+# Program:         hmmsearch
+# Version:         3.3.2 (Nov 2020)
+# Pipeline mode:   SEARCH
+# Query file:      /home/projects/DRAM/dbs_20210218/dbCAN-HMMdb-V9.txt
+# Target file:     ./genes_small
+# Option settings: hmmsearch --domtblout orfs_test/cazy_results.unprocessed.b6 --cpu 10 /home/projects/DRAM/dbs_20210218/dbCAN-HMMdb-V9.txt ./genes_small
+# Current dir:     /home/projects/DRAM/scratch_space_rory/nov19_custom_hmm/function_tests
+# Date:            Tue Nov 30 11:28:59 2021
+# [ok]
diff --git a/tests/data/unformatted_kofam.b6 b/tests/data/unformatted_kofam.b6
new file mode 100644
index 00000000..246559a3
--- /dev/null
+++ b/tests/data/unformatted_kofam.b6
@@ -0,0 +1,36 @@
+#                                                                                             --- full sequence --- -------------- this domain -------------   hmm coord   ali coord   env coord
+# target name                         accession   tlen query name           accession   qlen   E-value  score  bias   #  of  c-Evalue  i-Evalue  score  bias  from    to  from    to  from    to  acc description of target
+#                 ------------------- ---------- ----- -------------------- ---------- ----- --------- ------ ----- --- --- --------- --------- ------ ----- ----- ----- ----- ----- ----- ----- ---- ---------------------
+bin_1.scaffold_1 -            387 K00003               -            742   0.00058   11.0   0.1   1   2    0.0018      0.18    2.8   0.0   111   149   119   157    27   214 0.62 rank: A; rix:RO1_05310 Phosphoglycerate dehydrogenase and related dehydrogenases; K00001 (db=kegg)
+bin_1.scaffold_1 -            387 K00003               -            742   0.00058   11.0   0.1   2   2   0.00013     0.013    6.5   0.0   462   517   325   382   290   386 0.83 rank: A; rix:RO1_05310 Phosphoglycerate dehydrogenase and related dehydrogenases; K00001 (db=kegg)
+bin_1.scaffold_1 -            387 K00008               -            344   3.9e-10   31.9   0.1   1   1   5.6e-12   5.6e-10   31.4   0.1    22   133   110   222    48   247 0.81 rank: A; rix:RO1_05310 Phosphoglycerate dehydrogenase and related dehydrogenases; K00001 (db=kegg)
+bin_1.scaffold_1 -            387 K00008               -            315   8.3e-38  122.7   0.1   1   1   1.2e-39   1.2e-37  122.1   0.1    68   314    50   305    29   306 0.87 rank: A; rix:RO1_05310 Phosphoglycerate dehydrogenase and related dehydrogenases; K00001 (db=kegg)
+bin_1.scaffold_1 -            387 K00008               -            385   6.2e-44  143.0   0.0   1   1   8.6e-46   8.6e-44  142.6   0.0    80   373    26   317    13   328 0.83 rank: A; rix:RO1_05310 Phosphoglycerate dehydrogenase and related dehydrogenases; K00001 (db=kegg)
+bin_1.scaffold_1 -            387 K00001               -            509  1.2e-106  349.9   0.8   1   1  1.4e-108  1.4e-106  349.7   0.8    61   440     9   381     1   385 0.86 rank: A; rix:RO1_05310 Phosphoglycerate dehydrogenase and related dehydrogenases; K00001 (db=kegg)
+bin_1.scaffold_1 -            387 K00008               -            332   3.4e-57  186.1   0.1   1   1   4.8e-59   4.8e-57  185.6   0.1    45   330    22   305     1   306 0.83 rank: A; rix:RO1_05310 Phosphoglycerate dehydrogenase and related dehydrogenases; K00001 (db=kegg)
+bin_1.scaffold_1 -            387 K00008               -            462   6.8e-64  209.1   0.2   1   1   8.9e-66   8.9e-64  208.8   0.2    93   381    24   313     6   380 0.80 rank: A; rix:RO1_05310 Phosphoglycerate dehydrogenase and related dehydrogenases; K00001 (db=kegg)
+bin_1.scaffold_1 -            387 K00008               -            476   4.7e-09   28.2   0.0   1   1   6.8e-11   6.8e-09   27.7   0.0    50   153   120   228     9   251 0.71 rank: A; rix:RO1_05310 Phosphoglycerate dehydrogenase and related dehydrogenases; K00001 (db=kegg)
+bin_1.scaffold_1 -            387 K00008               -            324   4.1e-47  153.2   0.0   1   1   5.9e-49   5.9e-47  152.7   0.0    40   319    26   305    11   310 0.85 rank: A; rix:RO1_05310 Phosphoglycerate dehydrogenase and related dehydrogenases; K00001 (db=kegg)
+bin_1.scaffold_1 -            387 K00008               -            326   1.9e-51  167.7   0.0   1   1   2.9e-53   2.9e-51  167.1   0.0    36   320    22   305    11   313 0.86 rank: A; rix:RO1_05310 Phosphoglycerate dehydrogenase and related dehydrogenases; K00001 (db=kegg)
+bin_1.scaffold_2 -            338 K00008               -            342   7.1e-05   14.3   0.2   2   2     0.052       2.6   -0.7   0.0   112   138   229   255   194   271 0.69 rank: A; erb:A4V01_20675 type I glyceraldehyde-3-phosphate dehydrogenase; K00002 (db=kegg)
+bin_1.scaffold_2 -            338 K00002               -            617   1.3e-09   29.9   0.2   1   2   3.1e-10   3.1e-08   25.2   0.1   140   270    52   184    38   210 0.81 rank: A; erb:A4V01_20675 type I glyceraldehyde-3-phosphate dehydrogenase; K00002 (db=kegg)
+bin_1.scaffold_2 -            338 K00008               -            346   0.00016   13.2   0.2   1   1   9.8e-06   0.00049   11.5   0.2   169   271     4   127     1   142 0.65 rank: A; erb:A4V01_20675 type I glyceraldehyde-3-phosphate dehydrogenase; K00002 (db=kegg)
+bin_1.scaffold_2 -            338 K00008               -            478   2.1e-07   22.5   0.2   1   1     3e-09     3e-07   22.0   0.2   222   331     4   140     1   187 0.68 rank: A; erb:A4V01_20675 type I glyceraldehyde-3-phosphate dehydrogenase; K00002 (db=kegg)
+bin_1.scaffold_2 -            338 K00002               -            976   9.9e-21   66.6   0.2   1   1   9.9e-23   9.9e-21   66.6   0.2   714   877   139   318    22   337 0.73 rank: A; erb:A4V01_20675 type I glyceraldehyde-3-phosphate dehydrogenase; K00002 (db=kegg)
+bin_1.scaffold_2 -            338 K00002               -            343   1.7e-48  158.1   2.3   1   1   2.1e-50   2.1e-48  157.8   2.3     3   303     2   328     1   332 0.88 rank: A; erb:A4V01_20675 type I glyceraldehyde-3-phosphate dehydrogenase; K00002 (db=kegg)
+bin_1.scaffold_2 -            338 K00007               -            704   0.00018   12.8   0.1   1   1   2.3e-06   0.00023   12.4   0.1   200   301     1   129     1   181 0.70 rank: A; erb:A4V01_20675 type I glyceraldehyde-3-phosphate dehydrogenase; K00002 (db=kegg)
+bin_1.scaffold_2 -            338 K00008               -            342   7.1e-05   14.3   0.2   1   2     1e-05   0.00051   11.5   0.1   173   275     4   134     2   166 0.71 rank: A; erb:A4V01_20675 type I glyceraldehyde-3-phosphate dehydrogenase; K00002 (db=kegg)
+bin_1.scaffold_2 -            338 K00002               -            617   1.3e-09   29.9   0.2   2   2    0.0021      0.21    2.7   0.0   342   378   233   269   221   289 0.85 rank: A; erb:A4V01_20675 type I glyceraldehyde-3-phosphate dehydrogenase; K00002 (db=kegg)
+bin_1.scaffold_2 -            338 K00002               -            488  1.8e-138  455.0   2.6   1   1    2e-140    2e-138  454.9   2.6    96   424     1   334     1   337 0.95 rank: A; erb:A4V01_20675 type I glyceraldehyde-3-phosphate dehydrogenase; K00002 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_49 -            148 K00008               -            431     2e-05   16.6   0.0   1   1   2.4e-07   2.4e-05   16.3   0.0   160   211    50   101    31   109 0.94 rank: C; cthd:CDO33_15135 pyridoxamine 5'-phosphate oxidase (db=kegg)
+bin_1.scaffold_1 -            387 K00002               -            494   4.7e-32  103.8   0.1   1   1   6.7e-34   6.7e-32  103.3   0.1   160   354    49   252    44   258 0.94 rank: A; rix:RO1_05310 Phosphoglycerate dehydrogenase and related dehydrogenases; K00001 (db=kegg)
+#
+# Program:         hmmsearch
+# Version:         3.3.2 (Nov 2020)
+# Pipeline mode:   SEARCH
+# Query file:      /home/projects/DRAM/dbs_20210218/kofam_profiles.hmm
+# Target file:     ./genes_small
+# Option settings: hmmsearch --domtblout orfs_test/kofam_results.unprocessed.b6 --cpu 10 /home/projects/DRAM/dbs_20210218/kofam_profiles.hmm ./genes_small
+# Current dir:     /home/projects/DRAM/scratch_space_rory/nov19_custom_hmm/function_tests
+# Date:            Tue Nov 30 11:28:28 2021
+# [ok]
diff --git a/tests/data/unformatted_vogdb.b6 b/tests/data/unformatted_vogdb.b6
new file mode 100644
index 00000000..5af50a86
--- /dev/null
+++ b/tests/data/unformatted_vogdb.b6
@@ -0,0 +1,39 @@
+#                                                                                             --- full sequence --- -------------- this domain -------------   hmm coord   ali coord   env coord
+# target name                         accession   tlen query name           accession   qlen   E-value  score  bias   #  of  c-Evalue  i-Evalue  score  bias  from    to  from    to  from    to  acc description of target
+#                 ------------------- ---------- ----- -------------------- ---------- ----- --------- ------ ----- --- --- --------- --------- ------ ----- ----- ----- ----- ----- ----- ----- ---- ---------------------
+bin_1.scaffold_1 -            445 VOG00001             -            509   1.1e-56  185.1   2.1   1   1     5e-58     5e-56  182.9   2.1    11   335    19   316    10   344 0.87 rank: A; erb:A4V01_20660 cysteine--tRNA ligase; K01883 (db=kegg)
+bin_1.scaffold_2 -           1283 VOG00002             -           1308     2e-79  260.0   1.7   1   3     8e-64     8e-62  201.7   0.1   182   783   213   849   203   854 0.83 rank: A; erb:A4V01_05270 DNA-directed RNA polymerase subunit beta'; K03046 (db=kegg)
+bin_1.scaffold_2 -           1283 VOG00002             -           1308     2e-79  260.0   1.7   2   3   2.1e-12   2.1e-10   31.3   0.0   937   994   947  1006   923  1042 0.78 rank: A; erb:A4V01_05270 DNA-directed RNA polymerase subunit beta'; K03046 (db=kegg)
+bin_1.scaffold_2 -           1283 VOG00002             -           1308     2e-79  260.0   1.7   3   3   4.2e-08   4.2e-06   17.1   0.0  1218  1273  1164  1220  1157  1248 0.90 rank: A; erb:A4V01_05270 DNA-directed RNA polymerase subunit beta'; K03046 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_60 -            382 VOG00001             -            354   5.4e-19   61.4   0.0   1   1   5.4e-21   5.4e-19   61.4   0.0    99   176     1   132     1   369 0.67 rank: C; cls:CXIVA_18970 HipB; predicted transcriptional regulator (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_74 -            113 VOG00003             -            394   4.5e-05   14.9   0.2   1   1     1e-06   5.1e-05   14.8   0.2   144   196    28    91    15   108 0.72 rank: C; erb:A4V01_17195 PadR family transcriptional regulator; K10947 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_60 -            382 VOG00003             -            394   0.00022   12.6   0.0   1   2   1.7e-05   0.00084   10.8   0.0   141   170    16    45     3    63 0.88 rank: C; cls:CXIVA_18970 HipB; predicted transcriptional regulator (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_60 -            382 VOG00003             -            394   0.00022   12.6   0.0   2   2      0.07       3.5   -1.2   0.0   174   218   326   367   313   379 0.69 rank: C; cls:CXIVA_18970 HipB; predicted transcriptional regulator (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_11 -           1335 VOG00015             -            182   1.4e-06   21.4   0.4   1   3     0.046       4.6    0.1   0.1    97   133   189   231    99   249 0.73 rank: A; ebm:SG0102_29580 cas9; CRISPR-associated endonuclease Cas9; K09952 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_11 -           1335 VOG00015             -            182   1.4e-06   21.4   0.4   2   3      0.15        15   -1.6   3.9    64   138   364   442   329   455 0.56 rank: A; ebm:SG0102_29580 cas9; CRISPR-associated endonuclease Cas9; K09952 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_11 -           1335 VOG00015             -            182   1.4e-06   21.4   0.4   3   3   1.4e-08   1.4e-06   21.4   0.4    62   157   786   888   765   914 0.65 rank: A; ebm:SG0102_29580 cas9; CRISPR-associated endonuclease Cas9; K09952 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_81 -            196 VOG00044             -            291   0.00092   11.6   0.0   1   1   1.2e-05    0.0012   11.3   0.0   100   156   140   191    84   194 0.86 rank: C; erb:A4V01_20645 RNA polymerase subunit sigma; K03091 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_81 -            196 VOG00048             -            907   1.1e-06   20.1   9.0   1   1   6.1e-08   6.1e-06   17.7   9.0   314   528    18   191     3   193 0.76 rank: C; erb:A4V01_20645 RNA polymerase subunit sigma; K03091 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_25 -            235 VOG00062             -            920   1.4e-11   36.0   0.1   1   2   6.7e-08   1.7e-06   19.2   0.3   279   307    19    48     1    50 0.75 rank: A; erb:A4V01_05015 ABC transporter ATP-binding protein; K01996 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_25 -            235 VOG00062             -            920   1.4e-11   36.0   0.1   2   2   1.5e-06   3.8e-05   14.7   0.0   762   851   119   194    91   217 0.83 rank: A; erb:A4V01_05015 ABC transporter ATP-binding protein; K01996 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_26 -            259 VOG00062             -            920   2.6e-10   31.8   0.3   1   2   2.1e-06   5.3e-05   14.2   0.3   269   307    11    50     6    70 0.85 rank: A; erb:A4V01_05010 ABC transporter ATP-binding protein; K01995 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_26 -            259 VOG00062             -            920   2.6e-10   31.8   0.3   2   2   7.9e-07     2e-05   15.7   0.0   807   855   173   220   158   233 0.80 rank: A; erb:A4V01_05010 ABC transporter ATP-binding protein; K01995 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_20 -            337 VOG00062             -            920   3.7e-08   24.7   0.3   1   2     0.017      0.43    1.3   0.0   276   300    25    47    18    54 0.77 rank: A; csep:CP523_14490 ABC transporter ATP-binding protein; K15583 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_20 -            337 VOG00062             -            920   3.7e-08   24.7   0.3   2   2   1.7e-08   4.4e-07   21.2   0.0   779   849   161   220   158   225 0.94 rank: A; csep:CP523_14490 ABC transporter ATP-binding protein; K15583 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_19 -            324 VOG00062             -            920   9.5e-07   20.0   0.4   1   2     0.059       1.5   -0.5   0.1   291   306    46    61    29    68 0.84 rank: A; erb:A4V01_05040 ABC transporter ATP-binding protein; K10823 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_19 -            324 VOG00062             -            920   9.5e-07   20.0   0.4   2   2   1.5e-07   3.6e-06   18.1   0.0   779   847   155   212   154   222 0.89 rank: A; erb:A4V01_05040 ABC transporter ATP-binding protein; K10823 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_84 -            445 VOG00068             -            155    0.0017   11.7   1.6   1   2   7.6e-05    0.0076    9.6   0.0    42   126    74   151    29   158 0.69 rank: A; erb:A4V01_20660 cysteine--tRNA ligase; K01883 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_84 -            445 VOG00068             -            155    0.0017   11.7   1.6   2   2     0.048       4.8    0.5   0.4    64   103   377   425   296   437 0.64 rank: A; erb:A4V01_20660 cysteine--tRNA ligase; K01883 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_47 -            451 VOG00069             -            717     3e-07   21.9   0.1   1   2   5.9e-07   5.9e-05   14.3   0.0   340   405   290   351   289   360 0.82 rank: A; erb:A4V01_16625 23S rRNA (uracil-5-)-methyltransferase RumA; K03215 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_47 -            451 VOG00069             -            717     3e-07   21.9   0.1   2   2   6.9e-05    0.0069    7.4   0.0   168   252   345   411   341   425 0.68 rank: A; erb:A4V01_16625 23S rRNA (uracil-5-)-methyltransferase RumA; K03215 (db=kegg)
+LO_bin.102_LO_bin.102_scaffold_191_26 -            259 VOG00074             -            365   0.00054   11.9   0.0   1   1   9.7e-06   0.00097   11.0   0.0    98   121    30    53     6    66 0.82 rank: A; erb:A4V01_05010 ABC transporter ATP-binding protein; K01995 (db=kegg)
+#
+# Program:         hmmsearch
+# Version:         3.3.2 (Nov 2020)
+# Pipeline mode:   SEARCH
+# Query file:      /home/projects/DRAM/dbs_20210218/vog_latest_hmms.txt
+# Target file:     ./genes_small
+# Option settings: hmmsearch --domtblout orfs_test/vogdb_results.unprocessed.b6 --cpu 10 /home/projects/DRAM/dbs_20210218/vog_latest_hmms.txt ./genes_small
+# Current dir:     /home/projects/DRAM/scratch_space_rory/nov19_custom_hmm/function_tests
+# Date:            Tue Nov 30 11:29:40 2021
+# [ok]
diff --git a/tests/test_annotate_bins.py b/tests/test_annotate_bins.py
index a084f273..514a4377 100644
--- a/tests/test_annotate_bins.py
+++ b/tests/test_annotate_bins.py
@@ -14,11 +14,12 @@
 from mag_annotator.utils import make_mmseqs_db
 from mag_annotator.annotate_bins import filter_fasta, run_prodigal, get_best_hits, \
     get_reciprocal_best_hits, process_reciprocal_best_hits, get_kegg_description, get_uniref_description, \
-    get_basic_description, get_peptidase_description, get_sig, get_gene_data, get_unannotated, assign_grades, \
+    get_basic_description, get_peptidase_description, get_sig_row, get_gene_data, get_unannotated, assign_grades, \
     generate_annotated_fasta, create_annotated_fasta, generate_renamed_fasta, rename_fasta, run_trna_scan, \
     run_barrnap, do_blast_style_search, count_motifs, strip_endings, process_custom_dbs, get_dups, \
     parse_hmmsearch_domtblout, annotate_gff, make_gbk_from_gff_and_fasta, make_trnas_interval, make_rrnas_interval,\
-    add_intervals_to_gff, filter_db_locs
+    add_intervals_to_gff, vogdb_hmmscan_formater, generic_hmmscan_formater, dbcan_hmmscan_formater, \
+    kofam_hmmscan_formater
 
 
 @pytest.fixture()
@@ -185,11 +186,11 @@ def test_get_peptidase_description():
                                                                            ' YP_611907~'
 
 
-def test_get_sig():
-    assert not get_sig(1, 85, 100, 1)
-    assert get_sig(1, 86, 100, 1e-16)
-    assert not get_sig(1, 29, 100, 1e-20)
-
+def test_get_sig_row():
+    names = ['target_start', 'target_end', 'target_length', 'full_evalue']
+    assert not get_sig_row(pd.Series([1, 85, 100, 1], index=names))
+    assert get_sig_row(pd.Series([1, 86, 100, 1e-16], index=names))
+    assert not get_sig_row(pd.Series([1, 29, 100, 1e-20], index=names))
 
 @pytest.fixture()
 def phix_prodigal_genes():
@@ -463,7 +464,7 @@ def fake_phix_annotations():
 def fake_gff_loc():
     return os.path.join('tests', 'data', 'fake_gff.gff')
 
-
+#TODO
 def test_annotate_gff(annotated_fake_gff_loc, fake_phix_annotations, fake_gff_loc, tmpdir):
     gff_output = tmpdir.mkdir('gff_annotate')
     test_annotated_gff_loc = os.path.join(gff_output, 'annotated.gff')
@@ -471,6 +472,99 @@ def test_annotate_gff(annotated_fake_gff_loc, fake_phix_annotations, fake_gff_lo
     assert cmp(annotated_fake_gff_loc, test_annotated_gff_loc)
 
 
+def test_kofam_hmmscan_formater_dbcan():
+    output_expt = pd.DataFrame({
+        "bin_1.scaffold_1": ["K00001", "KO1; description 1"],
+        "bin_1.scaffold_2": ["K00002", "KO2; description 2"]
+    }, index=["ko_id", "kegg_hit"]).T
+    input_b6 = os.path.join('tests', 'data', 'unformatted_kofam.b6')
+    hits = parse_hmmsearch_domtblout(input_b6)
+    output_rcvd = kofam_hmmscan_formater(
+        hits=hits,
+        hmm_info_path=os.path.join('tests', 'data', 'hmm_thresholds.txt'),
+        top_hit=True,
+        use_dbcan2_thresholds=True)
+    output_rcvd.sort_index(inplace=True)
+    output_expt.sort_index(inplace=True)
+    assert output_rcvd.equals(output_expt), "Error in kofam_hmmscan_formater with dbcam cuts"
+
+
+def test_kofam_hmmscan_formater():
+    output_expt = pd.DataFrame({
+        "bin_1.scaffold_2": ["K00002", "KO2; description 2"]
+    }, index=["ko_id", "kegg_hit"]).T
+    input_b6 = os.path.join('tests', 'data', 'unformatted_kofam.b6')
+    hits = parse_hmmsearch_domtblout(input_b6)
+    output_rcvd = kofam_hmmscan_formater(
+        hits=hits,
+        hmm_info_path=os.path.join('tests', 'data', 'hmm_thresholds.txt'),
+        top_hit=True,
+        use_dbcan2_thresholds=False)
+    output_rcvd.sort_index(inplace=True)
+    output_expt.sort_index(inplace=True)
+    assert output_rcvd.equals(output_expt), "Error in kofam_hmmscan_formater"
+
+
+def test_generic_hmmscan_formater_custom_cuts():
+    output_expt = pd.DataFrame({
+        "bin_1.scaffold_2": ["K00002", "KO2; description 2"]
+    }, index=["test_id", "test_hits"]).T
+    input_b6 = os.path.join('tests', 'data', 'unformatted_kofam.b6')
+    hits = parse_hmmsearch_domtblout(input_b6)
+    output_rcvd = generic_hmmscan_formater(
+        hits,
+        hmm_info_path=os.path.join('tests', 'data', 'hmm_thresholds.txt'),
+        db_name='test',
+        top_hit=True)
+    output_rcvd.sort_index(inplace=True)
+    output_expt.sort_index(inplace=True)
+    assert output_rcvd.equals(output_expt), "Error in generic_hmmscan_formater, with custom cuts"
+
+
+def test_generic_hmmscan_formater():
+    output_expt = pd.DataFrame({
+        "bin_1.scaffold_1": ["K00001"],
+        "bin_1.scaffold_2": ["K00002"]
+    }, index=["test_id"]).T
+    input_b6 = os.path.join('tests', 'data', 'unformatted_kofam.b6')
+    hits = parse_hmmsearch_domtblout(input_b6)
+    output_rcvd = generic_hmmscan_formater(
+        hits,
+        db_name='test',
+        top_hit=True)
+    output_rcvd.sort_index(inplace=True)
+    output_expt.sort_index(inplace=True)
+    assert output_rcvd.equals(output_expt), "Error in generic_hmmscan_formater"
+
+
+def test_vogdb_hmmscan_formater():
+    #TODO Not Done
+    output_expt = pd.DataFrame({
+        "bin_1.scaffold_1": "VOG00001",
+        "bin_1.scaffold_2": "VOG00002"
+    }, index=["vogdb_id"]).T
+    input_b6 = os.path.join('tests', 'data', 'unformatted_vogdb.b6')
+    hits = parse_hmmsearch_domtblout(input_b6)
+    output_rcvd = vogdb_hmmscan_formater(hits=hits, db_name='vogdb')
+    output_rcvd.sort_index(inplace=True)
+    output_expt.sort_index(inplace=True)
+    assert output_rcvd.equals(output_expt), "Error in vogdb_hmmscan_formater"
+
+
+def test_dbcan_hmmscan_formater():
+    # TODO can we test with db-handler?
+    output_expt = pd.DataFrame({"bin_1.scaffold_1": "GT4; GT5",
+                           "bin_1.scaffold_2": "GH1"}, index=["cazy_id"]).T
+    input_b6 = os.path.join('tests', 'data', 'unformatted_cazy.b6')
+    hits = parse_hmmsearch_domtblout(input_b6)
+    output_rcvd = dbcan_hmmscan_formater(hits=hits, db_name='cazy')
+    #hits.sort_values('full_evalue').drop_duplicates(subset=["query_id"])
+    # output_rcvd
+    output_rcvd.sort_index(inplace=True)
+    output_expt.sort_index(inplace=True)
+    assert output_rcvd.equals(output_expt), "Error in dbcan_hmmscan_formater"
+
+
 test_gbk = """LOCUS       NC_001422.1   5386 bp   DNA   linear   ENV   %s
 FEATURES             Location/Qualifiers
      CDS             51..221
@@ -674,16 +768,3 @@ def test_add_intervals_to_gff(annotated_fake_gff_loc, tmpdir):
     gff = list(read_sequence(annotate_fake_gff_loc_w_rna, format='gff3'))
     assert type(gff) is list
     assert open(annotate_fake_gff_loc_w_rna).read() == annotated_fake_gff_w_rna
-
-
-def test_filter_db_locs():
-    test1 = filter_db_locs({'uniref': 'a fake locations', 'kegg': '/a/fake/loc', 'viral': ''})
-    assert test1 == {'kegg': '/a/fake/loc'}
-    with pytest.raises(ValueError):
-        test2 = filter_db_locs({'uniref': 'a fake locations', 'kegg': '/a/fake/loc', 'viral': ''}, low_mem_mode=True)
-    test3 = filter_db_locs({'kofam': '1', 'kofam_ko_list': '2', 'kegg': '/a/fake/loc', 'viral': ''}, low_mem_mode=True)
-    assert test3 == {'kofam': '1', 'kofam_ko_list': '2'}
-    test4 = filter_db_locs({'uniref': 'a fake locations', 'kegg': '/a/fake/loc', 'viral': ''}, use_uniref=True)
-    assert test4 == {'uniref': 'a fake locations', 'kegg': '/a/fake/loc'}
-    test5 = filter_db_locs({'kegg': '/a/fake/loc', 'viral': ''}, use_uniref=True)
-    assert test5 == {'kegg': '/a/fake/loc'}
diff --git a/tests/test_database_handler.py b/tests/test_database_handler.py
index 423e438e..b4d75e5e 100644
--- a/tests/test_database_handler.py
+++ b/tests/test_database_handler.py
@@ -1,32 +1,28 @@
 import pytest
+
 from os import path
-from mag_annotator.database_setup import create_description_db, KeggDescription
-from mag_annotator.database_handler import DatabaseHandler
 from sqlalchemy.orm.exc import NoResultFound
+import json
 
-
-@pytest.fixture()
-def db_loc(tmpdir):
-    db_loc = str(tmpdir.mkdir('test_db').join('test_db.sqlite'))
-    create_description_db(db_loc)
-    return db_loc
-
-
-def test_create_description_db(db_loc):
-    assert path.isfile(db_loc)
+from mag_annotator.database_setup import KeggDescription, create_description_db
+from mag_annotator.database_handler import DatabaseHandler
 
 
 @pytest.fixture()
-def db_handler(db_loc):
-    return DatabaseHandler(db_loc)
+def db_handler():
+    return DatabaseHandler(path.join('tests', 'data', 'test_CONFIG'))
 
 
-def test_init_db_handler():
-    with pytest.raises(ValueError):
-        DatabaseHandler('my.fake.path.123.db.sqlite')
+def test_db_hander_init():  # TODO: this
+    pass
 
+def test_add_descriptions_to_database(tmpdir):
+    db_handler = DatabaseHandler(path.join('tests', 'data', 'test_CONFIG'))
+    test_db = path.join(tmpdir.mkdir('test_db'), 'test_db.sqlite')
+    create_description_db(test_db)
+    db_handler.description_loc = test_db
+    db_handler.start_db_session()
 
-def test_add_descriptions_to_database(db_handler):
     kegg_entries = [{'id': 'K00001', 'description': 'The first KO'},
                     {'id': 'K00002', 'description': 'The second KO'}]
     db_handler.add_descriptions_to_database(kegg_entries, 'kegg_description')
@@ -38,7 +34,12 @@ def test_add_descriptions_to_database(db_handler):
 
 
 @pytest.fixture()
-def db_w_entries(db_handler):
+def db_w_entries(tmpdir):
+    db_handler = DatabaseHandler(path.join('tests', 'data', 'test_CONFIG'))
+    test_db = path.join(tmpdir.mkdir('test_db_w_entries'), 'test_db.sqlite')
+    create_description_db(test_db)
+    db_handler.description_loc = test_db
+    db_handler.start_db_session()
     db_handler.add_descriptions_to_database([{'id': 'K00001', 'description': 'The first KO'},
                                              {'id': 'K00002', 'description': 'The second KO'}], 'kegg_description',
                                             clear_table=True)
@@ -67,3 +68,147 @@ def test_get_descriptions(db_w_entries):
 def test_get_database_names(db_w_entries):
     names = db_w_entries.get_database_names()
     assert len(names) == 7
+
+
+def test_filter_db_locs():
+    db_handler = DatabaseHandler(path.join('tests', 'data', 'test_CONFIG'))
+    db_handler.db_locs = {'kegg': 'a/fake/location/kegg.fasta', 'kofam': 'a/fake/location/kofam.hmm',
+                          'kofam_ko_list': 'a/fake/location/kofam_ko_list.tsv.gz',
+                          'vogdb': 'a/fake/location/vogdb.hmm', 'pfam': 'tests/data/Pfam-A_subset.full.gz',
+                          'dbcan': None, 'peptidase': None, 'uniref': 'a/fake/location/uniref.fasta'}
+    db_handler.filter_db_locs(low_mem_mode=True)
+    assert db_handler.db_locs == {'kofam': 'a/fake/location/kofam.hmm',
+                                  'kofam_ko_list': 'a/fake/location/kofam_ko_list.tsv.gz',
+                                  'pfam': 'tests/data/Pfam-A_subset.full.gz',
+                                  'dbcan': None, 'peptidase': None}
+    db_handler.db_locs = {'kegg': 'a/fake/location/kegg.fasta', 'kofam': 'a/fake/location/kofam.hmm',
+                          'kofam_ko_list': 'a/fake/location/kofam_ko_list.tsv.gz',
+                          'vogdb': 'a/fake/location/vogdb.hmm', 'pfam': 'tests/data/Pfam-A_subset.full.gz',
+                          'dbcan': None, 'peptidase': None}
+    db_handler.filter_db_locs(use_vogdb=False)
+    assert db_handler.db_locs == {'kegg': 'a/fake/location/kegg.fasta', 'kofam': 'a/fake/location/kofam.hmm',
+                                  'kofam_ko_list': 'a/fake/location/kofam_ko_list.tsv.gz',
+                                  'pfam': 'tests/data/Pfam-A_subset.full.gz', 'dbcan': None, 'peptidase': None}
+    db_handler.db_locs = {'kegg': 'a/fake/location/kegg.fasta', 'kofam': 'a/fake/location/kofam.hmm',
+                          'kofam_ko_list': 'a/fake/location/kofam_ko_list.tsv.gz',
+                          'vogdb': 'a/fake/location/vogdb.hmm', 'pfam': 'tests/data/Pfam-A_subset.full.gz',
+                          'dbcan': None, 'peptidase': None, 'uniref': 'a/fake/location/uniref.fasta'}
+    db_handler.filter_db_locs(use_uniref=False)
+    assert db_handler.db_locs == {'kegg': 'a/fake/location/kegg.fasta', 'kofam': 'a/fake/location/kofam.hmm',
+                                  'kofam_ko_list': 'a/fake/location/kofam_ko_list.tsv.gz',
+                                  'vogdb': 'a/fake/location/vogdb.hmm', 'pfam': 'tests/data/Pfam-A_subset.full.gz',
+                                  'dbcan': None, 'peptidase': None}
+    # old version
+    db_handler.db_locs = {'uniref': 'a fake locations', 'kegg': '/a/fake/loc', 'viral': ''}
+    db_handler.filter_db_locs(master_list=('uniref', 'kegg'))
+    assert db_handler.db_locs == {'uniref': 'a fake locations', 'kegg': '/a/fake/loc'}
+    # test
+    db_handler.db_locs = {'uniref': 'a fake locations', 'kegg': '/a/fake/loc', 'viral': ''}
+    with pytest.raises(ValueError):
+        db_handler.filter_db_locs(low_mem_mode=True)
+    # test
+    db_handler.db_locs = {'kofam': '1', 'kofam_ko_list': '2', 'kegg': '/a/fake/loc', 'viral': ''}
+    db_handler.filter_db_locs(low_mem_mode=True)
+    assert db_handler.db_locs == {'kofam': '1', 'kofam_ko_list': '2', 'viral': ''}
+    # test
+    db_handler.db_locs = {'uniref': 'a fake locations', 'kegg': '/a/fake/loc', 'viral': ''}
+    db_handler.filter_db_locs(use_uniref=True)
+    assert db_handler.db_locs == {'uniref': 'a fake locations', 'kegg': '/a/fake/loc', 'viral': ''}
+    # test
+    db_handler.db_locs = {'uniref': 'a fake locations', 'kegg': '/a/fake/loc', 'viral': ''}
+    db_handler.filter_db_locs(use_uniref=False)
+    assert db_handler.db_locs == {'kegg': '/a/fake/loc', 'viral': ''}
+
+
+def test_make_header_dict_from_mmseqs_db():
+    test_headers = DatabaseHandler.make_header_dict_from_mmseqs_db(path.join('tests', 'data', 'NC_001422.mmsdb'))
+    assert test_headers == [{'id': 'NP_040705.1', 'description': 'NP_040705.1 B [Escherichia virus phiX174]'},
+                            {'id': 'NP_040704.1', 'description': 'NP_040704.1 A* [Escherichia virus phiX174]'},
+                            {'id': 'NP_040703.1', 'description': 'NP_040703.1 A [Escherichia virus phiX174]'},
+                            {'id': 'NP_040713.1', 'description': 'NP_040713.1 H [Escherichia virus phiX174]'},
+                            {'id': 'NP_040712.1', 'description': 'NP_040712.1 G [Escherichia virus phiX174]'},
+                            {'id': 'NP_040711.1', 'description': 'NP_040711.1 F [Escherichia virus phiX174]'},
+                            {'id': 'NP_040710.1', 'description': 'NP_040710.1 J [Escherichia virus phiX174]'},
+                            {'id': 'NP_040709.1', 'description': 'NP_040709.1 E [Escherichia virus phiX174]'},
+                            {'id': 'NP_040708.1', 'description': 'NP_040708.1 D [Escherichia virus phiX174]'},
+                            {'id': 'NP_040707.1', 'description': 'NP_040707.1 C [Escherichia virus phiX174]'},
+                            {'id': 'NP_040706.1', 'description': 'NP_040706.1 K [Escherichia virus phiX174]'}]
+
+
+def test_process_pfam_descriptions():
+    description_list = DatabaseHandler.process_pfam_descriptions(path.join('tests', 'data', 'Pfam-A_subset.hmm.dat.gz'))
+    assert description_list == [{'id': 'PF10417.9', 'description': 'C-terminal domain of 1-Cys peroxiredoxin'},
+                                {'id': 'PF12574.8', 'description': '120 KDa Rickettsia surface antigen'},
+                                {'id': 'PF09847.9', 'description': 'Membrane protein of 12 TMs'},
+                                {'id': 'PF00244.20', 'description': '14-3-3 protein'},
+                                {'id': 'PF16998.5', 'description': '17 kDa outer membrane surface antigen'}]
+
+
+def test_process_dbcan_descriptions():
+    description_list = DatabaseHandler.process_dbcan_descriptions(
+        path.join('tests', 'data', 'CAZyDB.07312019.fam-activities.subset.txt'))
+    assert description_list == [{'id': 'AA0', 'description': 'AA0'},
+                                {'id': 'AA10', 'description': 'AA10 (formerly CBM33) proteins are copper-dependent '
+                                                              'lytic polysaccharide monooxygenases (LPMOs); some '
+                                                              'proteins have been shown to act on chitin, others on '
+                                                              'cellulose; lytic cellulose monooxygenase '
+                                                              '(C1-hydroxylating) (EC 1.14.99.54); lytic cellulose '
+                                                              'monooxygenase (C4-dehydrogenating)(EC 1.14.99.56); '
+                                                              'lytic chitin monooxygenase (EC 1.14.99.53)'},
+                                {'id': 'AA11', 'description': 'AA11 proteins are copper-dependent lytic polysaccharide '
+                                                              'monooxygenases (LPMOs); cleavage of chitin chains with '
+                                                              'oxidation of C-1 has been demonstrated for a AA11 LPMO '
+                                                              'from Aspergillus oryzae;'},
+                                {'id': 'AA12', 'description': 'AA12 The pyrroloquinoline quinone-dependent oxidoreductase '
+                                                              'activity was demonstrated for the CC1G_09525 protein of '
+                                                              'Coprinopsis cinerea.'}]
+
+
+def test_process_vogdb_descriptions():
+    description_list = DatabaseHandler.process_vogdb_descriptions(
+        path.join('tests', 'data', 'vog_annotations_latest.subset.tsv.gz'))
+    assert description_list == [{'id': 'VOG00001', 'description': 'sp|Q5UQJ2|YR863_MIMIV Putative ankyrin repeat '
+                                                                  'protein R863; Xh'},
+                                {'id': 'VOG00002', 'description': 'sp|Q9J4Z6|V244_FOWPN Putative ankyrin repeat '
+                                                                  'protein FPV244; Xh'},
+                                {'id': 'VOG00003', 'description': 'sp|Q91FD6|388R_IIV6 Putative MSV199 '
+                                                                  'domain-containing protein 388R; Xu'},
+                                {'id': 'VOG00004', 'description': 'sp|P51704|RPC1_BPHC1 Repressor protein CI; Xu'},
+                                {'id': 'VOG00005', 'description': 'sp|P17766|POLG_PPVNA Genome polyprotein; XhXrXs'}]
+
+
+def test_export_config(tmpdir, db_handler):
+    test_exported_config = path.join(tmpdir.mkdir('test_export'), 'test_CONFIG')
+    db_handler.write_config(test_exported_config)
+    assert type(json.loads(open(test_exported_config).read())) is dict
+
+
+def test_print_database_locations(capsys, db_handler):
+    db_handler.print_database_locations()
+    out, err = capsys.readouterr()
+    assert 'Description db: ' in out
+    assert len(err) == 0
+
+
+def test_set_database_paths():
+    db_handler = DatabaseHandler()
+    # first test that adding nothing doesn't change CONFIG
+    pretest_db_dict = db_handler.db_locs
+    db_handler.set_database_paths(write_config=False)
+    assert pretest_db_dict == db_handler.db_locs
+    # test that adding something that doesn't exist throws error
+    test_fake_database = path.join('tests', 'data', 'fake_database.mmsdb')
+    with pytest.raises(ValueError):
+        db_handler.set_database_paths(kegg_db_loc=test_fake_database, write_config=False)
+    # test that adding something real is really added
+    kegg_loc = path.join('tests', 'data', 'fake_gff.gff')
+    db_handler.set_database_paths(kegg_db_loc=kegg_loc, write_config=False)
+    assert db_handler.db_locs['kegg'] == path.realpath(kegg_loc)
+
+
+def test_clear_config():
+    db_handler = DatabaseHandler()
+    db_handler.clear_config()
+    assert db_handler.db_locs == {}
+    assert db_handler.db_description_locs == {}
+    assert db_handler.description_loc is None
diff --git a/tests/test_database_processing.py b/tests/test_database_processing.py
index ffd79c9a..2a5e4772 100644
--- a/tests/test_database_processing.py
+++ b/tests/test_database_processing.py
@@ -1,36 +1,7 @@
 import pytest
-
 import os
-import json
-
-from mag_annotator.utils import get_database_locs
-from mag_annotator.database_processing import check_file_exists, make_header_dict_from_mmseqs_db, \
-    generate_modified_kegg_fasta, process_kegg_db, process_mmspro, process_pfam_descriptions, \
-    process_dbcan_descriptions, process_vogdb_descriptions, check_exists_and_add_to_location_dict, export_config, \
-    print_database_locations, set_database_paths
-
-
-def test_check_file_exists():
-    assert check_file_exists(None)
-    assert check_file_exists(os.path.join('tests', 'data', 'fake_gff.gff'))
-    with pytest.raises(ValueError):
-        check_file_exists(os.path.join('tests', 'data', 'a_fake_nonexistent_file.txt'))
-
-
-def test_make_header_dict_from_mmseqs_db():
-    test_headers = make_header_dict_from_mmseqs_db(os.path.join('tests', 'data', 'NC_001422.mmsdb'))
-    assert test_headers == [{'id': 'NP_040705.1', 'description': 'NP_040705.1 B [Escherichia virus phiX174]'},
-                            {'id': 'NP_040704.1', 'description': 'NP_040704.1 A* [Escherichia virus phiX174]'},
-                            {'id': 'NP_040703.1', 'description': 'NP_040703.1 A [Escherichia virus phiX174]'},
-                            {'id': 'NP_040713.1', 'description': 'NP_040713.1 H [Escherichia virus phiX174]'},
-                            {'id': 'NP_040712.1', 'description': 'NP_040712.1 G [Escherichia virus phiX174]'},
-                            {'id': 'NP_040711.1', 'description': 'NP_040711.1 F [Escherichia virus phiX174]'},
-                            {'id': 'NP_040710.1', 'description': 'NP_040710.1 J [Escherichia virus phiX174]'},
-                            {'id': 'NP_040709.1', 'description': 'NP_040709.1 E [Escherichia virus phiX174]'},
-                            {'id': 'NP_040708.1', 'description': 'NP_040708.1 D [Escherichia virus phiX174]'},
-                            {'id': 'NP_040707.1', 'description': 'NP_040707.1 C [Escherichia virus phiX174]'},
-                            {'id': 'NP_040706.1', 'description': 'NP_040706.1 K [Escherichia virus phiX174]'}]
 
+from mag_annotator.database_processing import generate_modified_kegg_fasta, process_kegg_db, process_mmspro
 
 @pytest.fixture()
 def phix_proteins():
@@ -64,97 +35,3 @@ def test_process_mmspro(tmpdir):
     assert os.path.isfile(os.path.join(processed_mmspro, 'fake.mmspro_h'))
     assert os.path.isfile(os.path.join(processed_mmspro, 'fake.mmspro.idx'))
     assert os.path.isfile(os.path.join(processed_mmspro, 'fake.mmspro_h.index'))
-
-
-def test_process_pfam_descriptions():
-    description_list = process_pfam_descriptions(os.path.join('tests', 'data', 'Pfam-A_subset.hmm.dat.gz'))
-    assert description_list == [{'id': 'PF10417.9', 'description': 'C-terminal domain of 1-Cys peroxiredoxin'},
-                                {'id': 'PF12574.8', 'description': '120 KDa Rickettsia surface antigen'},
-                                {'id': 'PF09847.9', 'description': 'Membrane protein of 12 TMs'},
-                                {'id': 'PF00244.20', 'description': '14-3-3 protein'},
-                                {'id': 'PF16998.5', 'description': '17 kDa outer membrane surface antigen'}]
-
-
-def test_process_dbcan_descriptions():
-    description_list = process_dbcan_descriptions(os.path.join('tests', 'data', 'CAZyDB.07312019.fam-activities.subset.txt'))
-    assert description_list == [{'id': 'AA0', 'description': 'AA0'},
-                                {'id': 'AA10', 'description': 'AA10 (formerly CBM33) proteins are copper-dependent '
-                                                              'lytic polysaccharide monooxygenases (LPMOs); some '
-                                                              'proteins have been shown to act on chitin, others on '
-                                                              'cellulose; lytic cellulose monooxygenase '
-                                                              '(C1-hydroxylating) (EC 1.14.99.54); lytic cellulose '
-                                                              'monooxygenase (C4-dehydrogenating)(EC 1.14.99.56); '
-                                                              'lytic chitin monooxygenase (EC 1.14.99.53)'},
-                                {'id': 'AA11', 'description': 'AA11 proteins are copper-dependent lytic polysaccharide '
-                                                              'monooxygenases (LPMOs); cleavage of chitin chains with '
-                                                              'oxidation of C-1 has been demonstrated for a AA11 LPMO '
-                                                              'from Aspergillus oryzae;'},
-                                {'id': 'AA12', 'description': 'AA12 The pyrroloquinoline quinone-dependent oxidoreductase '
-                                                              'activity was demonstrated for the CC1G_09525 protein of '
-                                                              'Coprinopsis cinerea.'}]
-
-
-def test_process_vogdb_descriptions():
-    description_list = process_vogdb_descriptions(os.path.join('tests', 'data', 'vog_annotations_latest.subset.tsv.gz'))
-    assert description_list == [{'id': 'VOG00001', 'description': 'sp|Q5UQJ2|YR863_MIMIV Putative ankyrin repeat '
-                                                                  'protein R863; Xh'},
-                                {'id': 'VOG00002', 'description': 'sp|Q9J4Z6|V244_FOWPN Putative ankyrin repeat '
-                                                                  'protein FPV244; Xh'},
-                                {'id': 'VOG00003', 'description': 'sp|Q91FD6|388R_IIV6 Putative MSV199 '
-                                                                  'domain-containing protein 388R; Xu'},
-                                {'id': 'VOG00004', 'description': 'sp|P51704|RPC1_BPHC1 Repressor protein CI; Xu'},
-                                {'id': 'VOG00005', 'description': 'sp|P17766|POLG_PPVNA Genome polyprotein; XhXrXs'}]
-
-
-def test_check_exists_and_add_to_location_dict(phix_proteins):
-    test_dict1 = check_exists_and_add_to_location_dict(phix_proteins, 'fake', {})
-    assert test_dict1 == {'fake': os.path.realpath(phix_proteins)}
-    test_dict2 = check_exists_and_add_to_location_dict(None, 'fake', {})
-    assert test_dict2 == {'fake': None}
-    test_dict3 = check_exists_and_add_to_location_dict(None, 'fake', {'fake': 'abc123'})
-    assert test_dict3 == {'fake': 'abc123'}
-
-
-def test_export_config(capsys, tmpdir):
-    export_config()
-    out, err = capsys.readouterr()
-    assert len(err) == 0
-    assert 'description_db' in out
-    assert out.startswith('{')
-    assert out.strip().endswith('}')
-
-    test_exported_config = os.path.join(tmpdir.mkdir('test_export'), 'test_CONFIG')
-    export_config(test_exported_config)
-    assert type(json.loads(open(test_exported_config).read())) is dict
-
-
-def test_print_database_locations(capsys):
-    print_database_locations()
-    out, err = capsys.readouterr()
-    assert 'Description db: ' in out
-    assert len(err) == 0
-
-
-def test_set_database_paths(tmpdir):
-    test_config_dir = tmpdir.mkdir('test_config')
-    # first test that adding nothing doesn't change CONFIG
-    test_config = os.path.join(test_config_dir, 'CONFIG')
-    pretest_db_dict = get_database_locs()
-    set_database_paths(config_loc=test_config)
-    test_db_dict = get_database_locs(test_config)
-    assert type(test_db_dict) is dict
-    assert pretest_db_dict == test_db_dict
-    # test that adding something that doesn't exist throws error
-    test_fake_database = os.path.join(test_config_dir, 'fake_database.mmsdb')
-    with pytest.raises(ValueError):
-        set_database_paths(kegg_db_loc=test_fake_database)
-    # test that adding something real is really added
-    kegg_loc = os.path.join('tests', 'data', 'fake_gff.gff')
-    set_database_paths(kegg_db_loc=kegg_loc, config_loc=test_config)
-    test_db_dict = get_database_locs(test_config)
-    assert test_db_dict['kegg'] == os.path.realpath(kegg_loc)
-    # test that adding something with use_current_locs False works
-    set_database_paths(kegg_db_loc=kegg_loc, config_loc=test_config, use_current_locs=False)
-    test_db_dict = get_database_locs(test_config)
-    assert test_db_dict['kegg'] == os.path.realpath(kegg_loc)
-    assert test_db_dict['description_db'] is None
diff --git a/tests/test_database_setup.py b/tests/test_database_setup.py
new file mode 100644
index 00000000..d8bcf212
--- /dev/null
+++ b/tests/test_database_setup.py
@@ -0,0 +1,16 @@
+import pytest
+
+from os import path
+
+from mag_annotator.database_setup import create_description_db
+
+
+@pytest.fixture()
+def db_loc(tmpdir):
+    db_loc = str(tmpdir.mkdir('test_db').join('test_db.sqlite'))
+    create_description_db(db_loc)
+    return db_loc
+
+
+def test_create_description_db(db_loc):
+    assert path.isfile(db_loc)
\ No newline at end of file
diff --git a/tests/test_summarize_genomes.py b/tests/test_summarize_genomes.py
index 91029160..41bfb9f7 100644
--- a/tests/test_summarize_genomes.py
+++ b/tests/test_summarize_genomes.py
@@ -5,12 +5,13 @@
 import networkx as nx
 import altair as alt
 
-from mag_annotator.summarize_genomes import get_ordered_uniques, fill_genome_summary_frame, summarize_rrnas, \
+from mag_annotator.summarize_genomes import fill_genome_summary_frame, summarize_rrnas, \
     summarize_trnas, build_module_net, get_module_step_coverage, make_module_coverage_df, make_module_coverage_frame, \
     make_module_coverage_heatmap, pairwise, first_open_paren_is_all, split_into_steps, is_ko, make_module_network, \
     get_module_coverage, make_etc_coverage_df, make_etc_coverage_heatmap, make_functional_df, make_functional_heatmap, \
     fill_liquor_dfs, make_liquor_heatmap, make_liquor_df, make_genome_summary, write_summarized_genomes_to_xlsx,\
     get_phylum_and_most_specific
+from mag_annotator.utils import get_ordered_uniques
 
 
 def test_get_ordered_uniques():
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 652465c4..94da8c74 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -4,8 +4,9 @@
 import json
 import pandas as pd
 
-from mag_annotator.utils import run_process, make_mmseqs_db, merge_files, multigrep, get_database_locs, \
-    get_config_loc, remove_prefix, remove_suffix, get_ids_from_row, get_genes_from_identifiers
+from mag_annotator.utils import run_process, make_mmseqs_db, merge_files, multigrep, remove_prefix, remove_suffix, \
+    get_ids_from_row
+from mag_annotator.pull_sequences import get_genes_from_identifiers
 
 
 def test_run_process():
@@ -85,18 +86,6 @@ def test_multigrep(multigrep_inputs):
     assert dict_['gene5'] == 'gene5 data including gene5'
 
 
-def test_get_config_loc():
-    test_config_loc = get_config_loc()
-    assert os.path.isfile(test_config_loc)
-    assert type(json.loads(open(test_config_loc).read())) is dict
-
-
-def test_get_database_locs():
-    test_database_locs = get_database_locs()
-    assert type(test_database_locs) is dict
-    assert 'description_db' in test_database_locs
-
-
 def test_remove_prefix():
     assert remove_prefix('prefix', 'pre') == 'fix'
     assert remove_prefix('postfix', 'pre') == 'postfix'