Repository: WrightonLabCSU/DRAM (language: python)

Issue:
DRAM-setup.py set_database_locations issue
Hello

I installed DRAM (v 1.4.1) using miniconda and imported existing databases using `DRAM-setup.py import_config --config_loc  my_old_config.txt`. 

Two of the databases (kofam_profiles.hmm and Pfam-A.hmm.dat.gz) were not found, so I tried `DRAM-setup.py set_database_locations --kofam_hmm_loc /work/ebg_lab/referenceDatabases/DRAM2/kofam_profiles.hmm --pfam_hmm_loc /work/ebg_lab/referenceDatabases/DRAM2/Pfam-A.hmm.dat.gz --update_description_db` to add them. However, I got the following error:

```
2022-12-08 16:42:12,932 - Logging to console
Traceback (most recent call last):
  File "/home/daniel.gittins/miniconda3/envs/dram/bin/DRAM-setup.py", line 184, in <module>
    args.func(**args_dict)
  File "/home/daniel.gittins/miniconda3/envs/dram/lib/python3.10/site-packages/mag_annotator/database_handler.py", line 434, in set_database_paths
    db_handler.set_database_paths(**kargs, write_config=True)
TypeError: DatabaseHandler.set_database_paths() got an unexpected keyword argument 'dbcan_fam_activities'
```

Do you know what is causing this issue?

For reference, the output of `DRAM-setup.py print_config`:

```
2022-12-09 08:26:36,326 - Logging to console
Processed search databases
KEGG db: None
KOfam db: None
KOfam KO list: /work/ebg_lab/referenceDatabases/DRAM2/kofam_ko_list.tsv
UniRef db: /work/ebg_lab/referenceDatabases/DRAM2/uniref90.20220122.mmsdb
Pfam db: /work/ebg_lab/referenceDatabases/DRAM2/pfam.mmspro
dbCAN db: /work/ebg_lab/referenceDatabases/DRAM2/dbCAN-HMMdb-V9.txt
RefSeq Viral db: /work/ebg_lab/referenceDatabases/DRAM2/refseq_viral.20220122.mmsdb
MEROPS peptidase db: /work/ebg_lab/referenceDatabases/DRAM2/peptidases.20220122.mmsdb
VOGDB db: /work/ebg_lab/referenceDatabases/DRAM2/vog_latest_hmms.txt
Descriptions of search database entries
Pfam hmm dat: None
dbCAN family activities: /work/ebg_lab/referenceDatabases/DRAM2/CAZyDB.07302020.fam-activities.txt
VOG annotations: /work/ebg_lab/referenceDatabases/DRAM2/vog_annotations_latest.tsv.gz

Description db: /work/ebg_lab/referenceDatabases/DRAM2/description_db.sqlite

DRAM distillation sheets
Genome summary form: /work/ebg_lab/referenceDatabases/DRAM2/genome_summary_form.20220122.tsv
Module step form: /work/ebg_lab/referenceDatabases/DRAM2/module_step_form.20220122.tsv
ETC module database: /work/ebg_lab/referenceDatabases/DRAM2/etc_mdoule_database.20220122.tsv
Function heatmap form: /work/ebg_lab/referenceDatabases/DRAM2/function_heatmap_form.20220122.tsv
AMG database: /work/ebg_lab/referenceDatabases/DRAM2/amg_database.20220122.tsv
```

Thank you





Discussion / hints:
Yes I do, this is a true bug and I do apologize. I will have it fixed in a moment. 

Fix patch:
diff --git a/mag_annotator/__init__.py b/mag_annotator/__init__.py
index 8e3c933c..98d186be 100644
--- a/mag_annotator/__init__.py
+++ b/mag_annotator/__init__.py
@@ -1,1 +1,1 @@
-__version__ = '1.4.1'
+__version__ = '1.4.2'
diff --git a/mag_annotator/annotate_bins.py b/mag_annotator/annotate_bins.py
index adb61e08..33e49ec6 100644
--- a/mag_annotator/annotate_bins.py
+++ b/mag_annotator/annotate_bins.py
@@ -16,13 +16,34 @@
 from shutil import rmtree, copy2
 import pandas as pd
 
-from mag_annotator.utils import run_process, make_mmseqs_db, merge_files, \
-    multigrep, remove_suffix, setup_logger, run_hmmscan, sig_scores, get_sig_row, \
-    generic_hmmscan_formater, get_reciprocal_best_hits, get_best_hits, BOUTFMT6_COLUMNS
+from mag_annotator.utils import (
+    run_process,
+    make_mmseqs_db,
+    merge_files,
+    multigrep,
+    remove_suffix,
+    setup_logger,
+    run_hmmscan,
+    sig_scores,
+    get_sig_row,
+    generic_hmmscan_formater,
+    get_reciprocal_best_hits,
+    get_best_hits,
+    BOUTFMT6_COLUMNS,
+)
 from mag_annotator.database_handler import DatabaseHandler
 from mag_annotator.fasta_dup_name_test import fastas_dup_check
 
-MAG_DBS_TO_ANNOTATE = ('kegg', 'kofam_hmm', 'kofam_ko_list', 'uniref', 'peptidase', 'pfam', 'dbcan', 'vogdb')
+MAG_DBS_TO_ANNOTATE = (
+    "kegg",
+    "kofam_hmm",
+    "kofam_ko_list",
+    "uniref",
+    "peptidase",
+    "pfam",
+    "dbcan",
+    "vogdb",
+)
 
 # TODO: Bind verbose to logging
 # TODO Exceptions are not fully handled by logging
@@ -30,47 +51,78 @@
 
 def filter_fasta(fasta_loc, min_len=5000, output_loc=None):
     """Removes sequences shorter than a set minimum from fasta files, outputs an object or to a file"""
-    kept_seqs = (seq for seq in read_sequence(fasta_loc, format='fasta') if len(seq) >= min_len)
+    kept_seqs = (
+        seq for seq in read_sequence(fasta_loc, format="fasta") if len(seq) >= min_len
+    )
     if output_loc is None:
         return list(kept_seqs)
     else:
-        write_sequence(kept_seqs, format='fasta', into=output_loc)
+        write_sequence(kept_seqs, format="fasta", into=output_loc)
 
 
-def run_prodigal(fasta_loc, output_dir, logger, mode='meta', trans_table='11', verbose=False):
+def run_prodigal(
+    fasta_loc, output_dir, logger, mode="meta", trans_table="11", verbose=False
+):
     """Runs the prodigal gene caller on a given fasta file, outputs resulting files to given directory"""
-    output_gff = path.join(output_dir, 'genes.gff')
-    output_fna = path.join(output_dir, 'genes.fna')
-    output_faa = path.join(output_dir, 'genes.faa')
-
-    run_process(['prodigal', '-i', fasta_loc, '-p', mode, '-g', trans_table, '-f', 'gff', '-o', output_gff, '-a',
-                 output_faa, '-d', output_fna], logger, verbose=verbose)
+    output_gff = path.join(output_dir, "genes.gff")
+    output_fna = path.join(output_dir, "genes.fna")
+    output_faa = path.join(output_dir, "genes.faa")
+
+    run_process(
+        [
+            "prodigal",
+            "-i",
+            fasta_loc,
+            "-p",
+            mode,
+            "-g",
+            trans_table,
+            "-f",
+            "gff",
+            "-o",
+            output_gff,
+            "-a",
+            output_faa,
+            "-d",
+            output_fna,
+        ],
+        logger,
+        verbose=verbose,
+    )
     return output_gff, output_fna, output_faa
 
 
-def process_reciprocal_best_hits(forward_output_loc, reverse_output_loc, target_prefix='target'):
+def process_reciprocal_best_hits(
+    forward_output_loc, reverse_output_loc, target_prefix="target"
+):
     """Process the forward and reverse best hits results to find reverse best hits
     Returns the query gene, target gene, if it was a reverse best hit, % identity, bit score and e-value
     """
-    forward_hits = pd.read_csv(forward_output_loc, sep='\t', header=None, names=BOUTFMT6_COLUMNS)
-    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')
+    forward_hits = pd.read_csv(
+        forward_output_loc, sep="\t", header=None, names=BOUTFMT6_COLUMNS
+    )
+    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:
             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,
-                '%s_bitScore' % target_prefix: row.bitScore,
-                '%s_eVal' % target_prefix:     row.eVal,
-                'index':                       row.name
-                }
-    hits = forward_hits.apply(check_hit, axis=1, result_type='expand')
+        return {
+            "%s_hit" % target_prefix: row.tId,
+            "%s_RBH" % target_prefix: rbh,
+            "%s_identity" % target_prefix: row.seqIdentity,
+            "%s_bitScore" % target_prefix: row.bitScore,
+            "%s_eVal" % target_prefix: row.eVal,
+            "index": row.name,
+        }
+
+    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.set_index("index", drop=True, inplace=True)
     hits.index.name = None
     return hits
 
@@ -82,15 +134,20 @@ def get_kegg_description(kegg_hits, header_dict):
     for kegg_hit in kegg_hits.kegg_hit:
         header = header_dict[kegg_hit]
         gene_description.append(header)
-        kos = re.findall(r'(K\d\d\d\d\d)', header)
+        kos = re.findall(r"(K\d\d\d\d\d)", header)
         if len(kos) == 0:
-            ko_list.append('')
+            ko_list.append("")
         else:
-            ko_list.append(','.join(kos))
+            ko_list.append(",".join(kos))
     # TODO: change kegg_id to kegg_genes_id so that people get an error and not the wrong identifier
-    new_df = pd.DataFrame([kegg_hits['kegg_hit'].values, ko_list, gene_description],
-                          index=['kegg_genes_id', 'ko_id', 'kegg_hit'], columns=kegg_hits.index)
-    return pd.concat([new_df.transpose(), kegg_hits.drop('kegg_hit', axis=1)], axis=1, sort=False)
+    new_df = pd.DataFrame(
+        [kegg_hits["kegg_hit"].values, ko_list, gene_description],
+        index=["kegg_genes_id", "ko_id", "kegg_hit"],
+        columns=kegg_hits.index,
+    )
+    return pd.concat(
+        [new_df.transpose(), kegg_hits.drop("kegg_hit", axis=1)], axis=1, sort=False
+    )
 
 
 def get_uniref_description(uniref_hits, header_dict):
@@ -101,27 +158,35 @@ def get_uniref_description(uniref_hits, header_dict):
     for uniref_hit in uniref_hits.uniref_hit:
         header = header_dict[uniref_hit]
         gene_description.append(header)
-        uniref_list.append(header[header.find('RepID=') + 6:])
-        gene_taxonomy.append(re.search(r'Tax=(.*?) (\S*?)=', header).group(1))
-    new_df = pd.DataFrame([uniref_list, gene_description, gene_taxonomy],
-                          index=['uniref_id', 'uniref_hit', 'uniref_taxonomy'],
-                          columns=uniref_hits.index)
-    return pd.concat([new_df.transpose(), uniref_hits.drop('uniref_hit', axis=1)], axis=1, sort=False)
+        uniref_list.append(header[header.find("RepID=") + 6 :])
+        gene_taxonomy.append(re.search(r"Tax=(.*?) (\S*?)=", header).group(1))
+    new_df = pd.DataFrame(
+        [uniref_list, gene_description, gene_taxonomy],
+        index=["uniref_id", "uniref_hit", "uniref_taxonomy"],
+        columns=uniref_hits.index,
+    )
+    return pd.concat(
+        [new_df.transpose(), uniref_hits.drop("uniref_hit", axis=1)], axis=1, sort=False
+    )
 
 
-def get_basic_description(hits, header_dict, db_name='viral'):
+def get_basic_description(hits, header_dict, db_name="viral"):
     """Get viral gene full descriptions based on headers (text before first space)"""
     hit_list = list()
     description = list()
-    [hit for hit in hits['%s_hit' % db_name]]
-    for hit in hits['%s_hit' % db_name]:
+    [hit for hit in hits["%s_hit" % db_name]]
+    for hit in hits["%s_hit" % db_name]:
         header = header_dict[hit]
         hit_list.append(hit)
         description.append(header)
-    new_df = pd.DataFrame([hit_list, description],
-                          index=['%s_id' % db_name, '%s_hit' % db_name],
-                          columns=hits.index)
-    return pd.concat([new_df.transpose(), hits.drop('%s_hit' % db_name, axis=1)], axis=1, sort=False)
+    new_df = pd.DataFrame(
+        [hit_list, description],
+        index=["%s_id" % db_name, "%s_hit" % db_name],
+        columns=hits.index,
+    )
+    return pd.concat(
+        [new_df.transpose(), hits.drop("%s_hit" % db_name, axis=1)], axis=1, sort=False
+    )
 
 
 def get_peptidase_description(peptidase_hits, header_dict):
@@ -131,51 +196,94 @@ def get_peptidase_description(peptidase_hits, header_dict):
     for peptidase_hit in peptidase_hits.peptidase_hit:
         header = header_dict[peptidase_hit]
         peptidase_list.append(peptidase_hit)
-        peptidase_family.append(re.search(r'#\w*.#', header).group()[1:-1])
+        peptidase_family.append(re.search(r"#\w*.#", header).group()[1:-1])
         peptidase_descirption.append(header)
-    new_df = pd.DataFrame([peptidase_list, peptidase_family, peptidase_descirption],
-                          index=['peptidase_id', 'peptidase_family', 'peptidase_hit'], columns=peptidase_hits.index)
-    return pd.concat([new_df.transpose(), peptidase_hits.drop('peptidase_hit', axis=1)], axis=1, sort=False)
+    new_df = pd.DataFrame(
+        [peptidase_list, peptidase_family, peptidase_descirption],
+        index=["peptidase_id", "peptidase_family", "peptidase_hit"],
+        columns=peptidase_hits.index,
+    )
+    return pd.concat(
+        [new_df.transpose(), peptidase_hits.drop("peptidase_hit", axis=1)],
+        axis=1,
+        sort=False,
+    )
 
 
-def run_mmseqs_profile_search(query_db, pfam_profile, output_loc, logger, output_prefix='mmpro_results', db_handler=None,
-                              threads=10, verbose=False):
+def run_mmseqs_profile_search(
+    query_db,
+    pfam_profile,
+    output_loc,
+    logger,
+    output_prefix="mmpro_results",
+    db_handler=None,
+    threads=10,
+    verbose=False,
+):
     """Use mmseqs to run a search against pfam, currently keeping all hits and not doing any extra filtering"""
-    tmp_dir = path.join(output_loc, 'tmp')
-    output_db = path.join(output_loc, '%s.mmsdb' % output_prefix)
-    run_process(['mmseqs', 'search', query_db, pfam_profile, output_db, tmp_dir, '-k', '5', '-s', '7', '--threads',
-                 str(threads)], logger, verbose=verbose)
-    output_loc = path.join(output_loc, '%s_output.b6' % output_prefix)
-    run_process(['mmseqs', 'convertalis', query_db, pfam_profile, output_db, output_loc], logger, verbose=verbose)
-    pfam_results = pd.read_csv(output_loc, sep='\t', header=None, names=BOUTFMT6_COLUMNS)
+    tmp_dir = path.join(output_loc, "tmp")
+    output_db = path.join(output_loc, "%s.mmsdb" % output_prefix)
+    run_process(
+        [
+            "mmseqs",
+            "search",
+            query_db,
+            pfam_profile,
+            output_db,
+            tmp_dir,
+            "-k",
+            "5",
+            "-s",
+            "7",
+            "--threads",
+            str(threads),
+        ],
+        logger,
+        verbose=verbose,
+    )
+    output_loc = path.join(output_loc, "%s_output.b6" % output_prefix)
+    run_process(
+        ["mmseqs", "convertalis", query_db, pfam_profile, output_db, output_loc],
+        logger,
+        verbose=verbose,
+    )
+    pfam_results = pd.read_csv(
+        output_loc, sep="\t", header=None, names=BOUTFMT6_COLUMNS
+    )
     if pfam_results.shape[0] > 0:
         pfam_dict = dict()
         if db_handler is not None:
-            pfam_descriptions = db_handler.get_descriptions(set(pfam_results.tId), '%s_description' % output_prefix)
+            pfam_descriptions = db_handler.get_descriptions(
+                set(pfam_results.tId), "%s_description" % output_prefix
+            )
         else:
             pfam_descriptions = {}
-        for gene, pfam_frame in pfam_results.groupby('qId'):
+        for gene, pfam_frame in pfam_results.groupby("qId"):
             if len(pfam_descriptions) < 1:
-                pfam_dict[gene] = '; '.join(pfam_frame.tId)
+                pfam_dict[gene] = "; ".join(pfam_frame.tId)
             else:
-                pfam_dict[gene] = '; '.join(['%s [%s]' % (pfam_descriptions[ascession], ascession)
-                                             for ascession in pfam_frame.tId])
+                pfam_dict[gene] = "; ".join(
+                    [
+                        "%s [%s]" % (pfam_descriptions[ascession], ascession)
+                        for ascession in pfam_frame.tId
+                    ]
+                )
         return pd.DataFrame(pfam_dict, index=[f"{output_prefix}_hits"]).T
     else:
         return pd.DataFrame(columns=[f"{output_prefix}_hits"])
 
 
-def find_best_dbcan_hit(genome:str, group:pd.DataFrame):
-    group['perc_cov'] = group.apply(
-        lambda x: (x['target_end'] - x['target_start']) / x['target_length'],
-        axis=1)
-    group.sort_values('perc_cov', inplace=True)
+def find_best_dbcan_hit(genome: str, group: pd.DataFrame):
+    group["perc_cov"] = group.apply(
+        lambda x: (x["target_end"] - x["target_start"]) / x["target_length"], axis=1
+    )
+    group.sort_values("perc_cov", inplace=True)
     group.columns
-    group.sort_values('full_evalue', inplace=True)
-    return group.iloc[0]['target_id']
+    group.sort_values("full_evalue", inplace=True)
+    return group.iloc[0]["target_id"]
 
 
-def dbcan_hmmscan_formater(hits:pd.DataFrame,  db_name:str, db_handler=None):
+def dbcan_hmmscan_formater(hits: pd.DataFrame, db_name: str, db_handler=None):
     """
     format the ouput of the dbcan database.
 
@@ -192,36 +300,44 @@ def dbcan_hmmscan_formater(hits:pd.DataFrame,  db_name:str, db_handler=None):
     if len(hits_sig) == 0:
         # if nothing significant then return nothing, don't get descriptions
         return pd.DataFrame()
-    hit_groups = hits_sig.groupby('query_id')
+    hit_groups = hits_sig.groupby("query_id")
     all_hits = hit_groups.apply(
-        lambda x: '; '.join(x['target_id'].apply(lambda y:y[:-4]).unique())
+        lambda x: "; ".join(x["target_id"].apply(lambda y: y[:-4]).unique())
     )
     hits_df = pd.DataFrame(all_hits)
     hits_df.columns = [f"{db_name}_ids"]
-    def description_pull(x:str):
-         id_list = [re.findall('^[A-Z]*[0-9]*', str(x))[0] for x in  x.split('; ')],
-         id_list = [y for x in  id_list for y in x if len(x) > 0]
-         description_list = db_handler.get_descriptions(id_list, 'dbcan_description').values()
-         description_str = '; '.join(description_list)
-         return description_str
+
+    def description_pull(x: str):
+        id_list = ([re.findall("^[A-Z]*[0-9]*", str(x))[0] for x in x.split("; ")],)
+        id_list = [y for x in id_list for y in x if len(x) > 0]
+        description_list = db_handler.get_descriptions(
+            id_list, "dbcan_description"
+        ).values()
+        description_str = "; ".join(description_list)
+        return description_str
+
     if db_handler is not None:
         hits_df[f"{db_name}_hits"] = hits_df[f"{db_name}_ids"].apply(description_pull)
         hits_df[f"{db_name}_subfam_ec"] = hits_df[f"{db_name}_ids"].apply(
-            lambda x:'; '.join(
+            lambda x: "; ".join(
                 db_handler.get_descriptions(
-                    x.split('; '),
-                    'dbcan_description',
-                    description_name='ec'
-                ).values()))
-    hits_df[f'{db_name}_best_hit'] = [find_best_dbcan_hit(*i) for i in hit_groups]
+                    x.split("; "), "dbcan_description", description_name="ec"
+                ).values()
+            )
+        )
+    hits_df[f"{db_name}_best_hit"] = [find_best_dbcan_hit(*i) for i in hit_groups]
     hits_df.rename_axis(None, inplace=True)
     hits_df.columns
     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)
+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:
@@ -230,22 +346,24 @@ def kofam_hmmscan_formater(hits:pd.DataFrame, hmm_info_path:str=None, use_dbcan2
     if len(hits_sig) == 0:
         return pd.DataFrame()
     kegg_dict = dict()
-    for gene, frame in hits_sig.groupby('query_id'):
+    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']]
+            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()
-
+            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 vogdb_hmmscan_formater(hits:pd.DataFrame,  db_name:str, logger:logging.Logger,
-                           db_handler=None):
+def vogdb_hmmscan_formater(
+    hits: pd.DataFrame, db_name: str, logger: logging.Logger, db_handler=None
+):
     categories_col = f"{db_name}_categories"
     id_col = f"{db_name}_id"
     description_col = f"{db_name}_description"
@@ -255,23 +373,29 @@ def vogdb_hmmscan_formater(hits:pd.DataFrame,  db_name:str, logger:logging.Logge
         # if nothing significant then return nothing, don't get descriptions
         return pd.DataFrame(columns=[categories_col, id_col, description_col])
     # Get the best hits
-    hits_best = hits_sig.sort_values('full_evalue').drop_duplicates(subset=["query_id"])
+    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()
+        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(), description_col),
-            index=[desc_col]).T
-        categories = descriptions[desc_col].apply(lambda x: x.split('; ')[-1])
+            db_handler.get_descriptions(
+                hits_best["target_id"].unique(), description_col
+            ),
+            index=[desc_col],
+        ).T
+        categories = descriptions[desc_col].apply(lambda x: x.split("; ")[-1])
         descriptions[categories_col] = 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)
+            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': id_col}, inplace=True)
+    hits_df.rename(columns={"target_id": id_col}, inplace=True)
     return hits_df
 
 
@@ -280,177 +404,230 @@ def get_gene_data(fasta_loc):
     Based on idba_ud 'scaffold_#' scaffold names with gene name after
     """
     df_dict = dict()
-    for seq in read_sequence(fasta_loc, format='fasta'):
-        split_label = seq.metadata['id'].split('_')
-        scaffold = '_'.join(split_label[:-1])
+    for seq in read_sequence(fasta_loc, format="fasta"):
+        split_label = seq.metadata["id"].split("_")
+        scaffold = "_".join(split_label[:-1])
         gene_position = split_label[-1]
-        start_position, end_position, strandedness = seq.metadata['description'].split('#')[1:4]
-        df_dict[seq.metadata['id']] = [scaffold, int(gene_position), int(start_position), int(end_position),
-                                       int(strandedness)]
-    return pd.DataFrame.from_dict(df_dict, orient='index', columns=['scaffold', 'gene_position', 'start_position',
-                                                                    'end_position', 'strandedness'])
+        start_position, end_position, strandedness = seq.metadata["description"].split(
+            "#"
+        )[1:4]
+        df_dict[seq.metadata["id"]] = [
+            scaffold,
+            int(gene_position),
+            int(start_position),
+            int(end_position),
+            int(strandedness),
+        ]
+    return pd.DataFrame.from_dict(
+        df_dict,
+        orient="index",
+        columns=[
+            "scaffold",
+            "gene_position",
+            "start_position",
+            "end_position",
+            "strandedness",
+        ],
+    )
 
 
 def get_unannotated(fasta_loc, annotations):
     """Get the genes from the fasta which did not get any annotations"""
-    return [seq.metadata['id'] for seq in read_sequence(fasta_loc, format='fasta')
-            if seq.metadata['id'] not in annotations]
+    return [
+        seq.metadata["id"]
+        for seq in read_sequence(fasta_loc, format="fasta")
+        if seq.metadata["id"] not in annotations
+    ]
 
 
 def assign_grades(annotations):
     """Grade genes based on reverse best hits to KEGG, UniRef and Pfam"""
     grades = dict()
     for gene, row in annotations.iterrows():
-        if row.get('kegg_RBH') is True:
-            rank = 'A'
-        elif row.get('uniref_RBH') is True:
-            rank = 'B'
-        elif not pd.isna(row.get('kegg_hit')) or not pd.isna(row.get('uniref_hit')):
-            rank = 'C'
-        elif not pd.isna(row.get('pfam_hits')) or not pd.isna(row.get('cazy_hits'))\
-                or not pd.isna(row.get('peptidase_hit')):
-            rank = 'D'
+        if row.get("kegg_RBH") is True:
+            rank = "A"
+        elif row.get("uniref_RBH") is True:
+            rank = "B"
+        elif not pd.isna(row.get("kegg_hit")) or not pd.isna(row.get("uniref_hit")):
+            rank = "C"
+        elif (
+            not pd.isna(row.get("pfam_hits"))
+            or not pd.isna(row.get("cazy_hits"))
+            or not pd.isna(row.get("peptidase_hit"))
+        ):
+            rank = "D"
         else:
-            rank = 'E'
+            rank = "E"
         grades[gene] = rank
-    return pd.DataFrame(grades, index=['rank']).T
+    return pd.DataFrame(grades, index=["rank"]).T
 
 
-def generate_annotated_fasta(input_fasta, annotations, verbosity='short', name=None):
+def generate_annotated_fasta(input_fasta, annotations, verbosity="short", name=None):
     """Generates fasta entries with added annotation information to the header of a fasta
     either add best annotation (based on grade) (verbosity = short) or all annotations (verbosity = long)
     """
-    for seq in read_sequence(input_fasta, format='fasta'):
-        annotation = annotations.loc[seq.metadata['id']]
-        if 'rank' in annotations.columns and verbosity == 'short':
-            annotation_str = 'rank: %s' % annotation['rank']
-            if annotation['rank'] == 'A':
-                annotation_str += '; %s (db=%s)' % (annotation.kegg_hit, 'kegg')
-            elif annotation['rank'] == 'B':
-                annotation_str += '; %s (db=%s)' % (annotation.uniref_hit, 'uniref')
-            elif annotation['rank'] == 'C':
-                if 'kegg_hit' in annotation:
-                    annotation_str += '; %s (db=%s)' % (annotation.kegg_hit, 'kegg')
-                if 'uniref_hit' in annotation:
-                    annotation_str += '; %s (db=%s)' % (annotation.uniref_hit, 'uniref')
-            elif annotation['rank'] == 'D':
-                annotation_str += '; %s (db=%s)' % (annotation.pfam_hits, 'pfam')
+    for seq in read_sequence(input_fasta, format="fasta"):
+        annotation = annotations.loc[seq.metadata["id"]]
+        if "rank" in annotations.columns and verbosity == "short":
+            annotation_str = "rank: %s" % annotation["rank"]
+            if annotation["rank"] == "A":
+                annotation_str += "; %s (db=%s)" % (annotation.kegg_hit, "kegg")
+            elif annotation["rank"] == "B":
+                annotation_str += "; %s (db=%s)" % (annotation.uniref_hit, "uniref")
+            elif annotation["rank"] == "C":
+                if "kegg_hit" in annotation:
+                    annotation_str += "; %s (db=%s)" % (annotation.kegg_hit, "kegg")
+                if "uniref_hit" in annotation:
+                    annotation_str += "; %s (db=%s)" % (annotation.uniref_hit, "uniref")
+            elif annotation["rank"] == "D":
+                annotation_str += "; %s (db=%s)" % (annotation.pfam_hits, "pfam")
             else:
                 pass
         else:
             annotation_list = []
-            if 'rank' in annotations.columns:
-                annotation_list += ['rank: %s' % annotation['rank']]
-            if 'kegg_hit' in annotations.columns:
+            if "rank" in annotations.columns:
+                annotation_list += ["rank: %s" % annotation["rank"]]
+            if "kegg_hit" in annotations.columns:
                 if not pd.isna(annotation.kegg_hit):
-                    annotation_list += ['%s (db=%s)' % (annotation.kegg_hit, 'kegg')]
-            if 'uniref_hit' in annotations.columns:
+                    annotation_list += ["%s (db=%s)" % (annotation.kegg_hit, "kegg")]
+            if "uniref_hit" in annotations.columns:
                 if not pd.isna(annotation.uniref_hit):
-                    annotation_list += ['%s (db=%s)' % (annotation.uniref_hit, 'uniref')]
-            if 'pfam_hits' in annotations.columns:
+                    annotation_list += [
+                        "%s (db=%s)" % (annotation.uniref_hit, "uniref")
+                    ]
+            if "pfam_hits" in annotations.columns:
                 if not pd.isna(annotation.pfam_hits):
-                    annotation_list += ['%s (db=%s)' % (annotation.pfam_hits, 'pfam')]
-            if 'bin_taxonomy' in annotations.columns:
+                    annotation_list += ["%s (db=%s)" % (annotation.pfam_hits, "pfam")]
+            if "bin_taxonomy" in annotations.columns:
                 if not pd.isna(annotation.bin_taxonomy):
                     annotation_list += [annotation.bin_taxonomy]
-            annotation_str = '; '.join(annotation_list)
+            annotation_str = "; ".join(annotation_list)
         if name is not None:
-            seq.metadata['id'] = '%s_%s' % (name, seq.metadata['id'])
-        seq.metadata['description'] = annotation_str
+            seq.metadata["id"] = "%s_%s" % (name, seq.metadata["id"])
+        seq.metadata["description"] = annotation_str
         yield seq
 
 
-def create_annotated_fasta(input_fasta, annotations, output_fasta, verbosity='short', name=None):
+def create_annotated_fasta(
+    input_fasta, annotations, output_fasta, verbosity="short", name=None
+):
     """For use with genes files, added annotations"""
 
-    write_sequence( generate_annotated_fasta(input_fasta, annotations, verbosity, name), format='fasta', into=output_fasta)
+    write_sequence(
+        generate_annotated_fasta(input_fasta, annotations, verbosity, name),
+        format="fasta",
+        into=output_fasta,
+    )
 
 
 def generate_renamed_fasta(input_fasta, prefix):
     """For use with scaffolds files, merges together bins with fasta name added as a prefix to the file"""
-    for seq in read_sequence(input_fasta, format='fasta'):
-        seq.metadata['id'] = '%s_%s' % (prefix, seq.metadata['id'])
+    for seq in read_sequence(input_fasta, format="fasta"):
+        seq.metadata["id"] = "%s_%s" % (prefix, seq.metadata["id"])
         yield seq
 
 
 def rename_fasta(input_fasta, output_fasta, prefix):
     """See above"""
-    write_sequence(generate_renamed_fasta(input_fasta, prefix), format='fasta', into=output_fasta)
+    write_sequence(
+        generate_renamed_fasta(input_fasta, prefix), format="fasta", into=output_fasta
+    )
 
 
 # TODO: add annotations that don't end with '_id'
 def annotate_gff(input_gff, output_gff, annotations, prefix=None):
     """Go through a gff and add a prefix to the scaffold and gene number for all ID's"""
     f = open(input_gff)
-    o = open(output_gff, 'w')
+    o = open(output_gff, "w")
     for line in f:
         line = line.strip()
-        if not line.startswith('#'):
+        if not line.startswith("#"):
             # replace id with new name
-            old_scaffold = line.split('\t')[0]
-            match = re.search(r'ID=\d*_\d*;', line)
-            gene_number = match.group().split('_')[-1][:-1]
-            old_gene_name = '%s_%s' % (old_scaffold, gene_number)
+            old_scaffold = line.split("\t")[0]
+            match = re.search(r"ID=\d*_\d*;", line)
+            gene_number = match.group().split("_")[-1][:-1]
+            old_gene_name = "%s_%s" % (old_scaffold, gene_number)
             if prefix is not None:  # add prefix to line name and gene name if given
-                line = '%s_%s' % (prefix, line)
-                gene_name = '%s_%s' % (prefix, old_gene_name)
+                line = "%s_%s" % (prefix, line)
+                gene_name = "%s_%s" % (prefix, old_gene_name)
             else:
                 gene_name = old_gene_name
-            line = re.sub(r'ID=\d*_\d*;', 'ID=%s;' % gene_name, line)
+            line = re.sub(r"ID=\d*_\d*;", "ID=%s;" % gene_name, line)
             # 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().replace(':', '_'))
-                                    for key, values in annotations_to_add.items()
-                                    if not pd.isna(values)
-                                    for value in values.split('; ')
-                                    if value != '']
+            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().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)
+                line += "%s;" % ";".join(database_information)
+        o.write("%s\n" % line)
 
 
-def make_gbk_from_gff_and_fasta(gff_loc='genes.gff', fasta_loc='scaffolds.fna', faa_loc='genes.faa', output_gbk=None):
+def make_gbk_from_gff_and_fasta(
+    gff_loc="genes.gff", fasta_loc="scaffolds.fna", faa_loc="genes.faa", output_gbk=None
+):
     # filter scaffolds out of fasta which are not in genes.gff
-    gff_frame = pd.read_csv(gff_loc, sep='\t', comment='#', header=None)
+    gff_frame = pd.read_csv(gff_loc, sep="\t", comment="#", header=None)
     gff_scaffolds = set(gff_frame[0])
     capture_fasta = io.StringIO()
-    write_sequence((i for i in read_sequence(fasta_loc, format='fasta') if i.metadata['id'] in gff_scaffolds),
-                   into=capture_fasta, format='fasta')
+    write_sequence(
+        (
+            i
+            for i in read_sequence(fasta_loc, format="fasta")
+            if i.metadata["id"] in gff_scaffolds
+        ),
+        into=capture_fasta,
+        format="fasta",
+    )
     # scikit-bio can make a genbank for a concatenated gff and fasta file so make this first
     gff = open(gff_loc)
-    aa_records = {i.metadata['id']: i for i in read_sequence(faa_loc, format='fasta')}
+    aa_records = {i.metadata["id"]: i for i in read_sequence(faa_loc, format="fasta")}
     # the gff with fasta format is the gff then ##FASTA then the fasta
-    concat_gff = '%s\n##FASTA\n%s' % (gff.read(), capture_fasta.getvalue())
+    concat_gff = "%s\n##FASTA\n%s" % (gff.read(), capture_fasta.getvalue())
 
     # now we can only ready one record from the gff/fasta hybrid at a time
     # read a record at a time till end of the fasta
-    genbank_records = ''
+    genbank_records = ""
     # +1 because there is no \n> for the first line in the file and +1 because we are indexing starting at 1
-    for i in range(1, capture_fasta.getvalue().count('\n>') + 1 + 1):
-        seq = read_sequence(io.StringIO(concat_gff), format='gff3', into=Sequence, seq_num=i)
-        seq.metadata['LOCUS'] = {'locus_name': seq.metadata['id'],
-                                 'size': len(seq),
-                                 'unit': 'bp',
-                                 'mol_type': 'DNA',
-                                 'shape': 'linear',
-                                 'division': 'ENV',
-                                 'date': time.strftime('%d-%b-%Y').upper()}
-        seq_bounds = (seq.interval_metadata.lower_bound, seq.interval_metadata.upper_bound)
+    for i in range(1, capture_fasta.getvalue().count("\n>") + 1 + 1):
+        seq = read_sequence(
+            io.StringIO(concat_gff), format="gff3", into=Sequence, seq_num=i
+        )
+        seq.metadata["LOCUS"] = {
+            "locus_name": seq.metadata["id"],
+            "size": len(seq),
+            "unit": "bp",
+            "mol_type": "DNA",
+            "shape": "linear",
+            "division": "ENV",
+            "date": time.strftime("%d-%b-%Y").upper(),
+        }
+        seq_bounds = (
+            seq.interval_metadata.lower_bound,
+            seq.interval_metadata.upper_bound,
+        )
         for interval in seq.interval_metadata.query([seq_bounds]):
-            interval.metadata['gene'] = interval.metadata['ID']
-            if interval.metadata['ID'] in aa_records:
-                interval.metadata['translation'] = aa_records[interval.metadata['ID']]
+            interval.metadata["gene"] = interval.metadata["ID"]
+            if interval.metadata["ID"] in aa_records:
+                interval.metadata["translation"] = aa_records[interval.metadata["ID"]]
         # need to capture genbank output so we can combine into a multi-genbank file
         capture_print = io.StringIO()
-        seq.write(capture_print, format='genbank')
+        seq.write(capture_print, format="genbank")
         genbank_records += capture_print.getvalue()
 
     # if no output then return as string, otherwise write write to output location
     if output_gbk is None:
         return genbank_records
     else:
-        open(output_gbk, 'w').write(genbank_records)
+        open(output_gbk, "w").write(genbank_records)
 
 
 def get_dups(columns):
@@ -467,75 +644,130 @@ def get_dups(columns):
 
 def run_trna_scan(fasta, tmp_dir, fasta_name, logger, threads=10, verbose=True):
     """Run tRNAscan-SE on scaffolds and create a table of tRNAs as a separate output"""
-    raw_trnas = path.join(tmp_dir, 'raw_trnas.txt')
-    run_process(['tRNAscan-SE', '-G', '-o', raw_trnas, '--thread', str(threads), fasta], logger, verbose=verbose)
+    raw_trnas = path.join(tmp_dir, "raw_trnas.txt")
+    run_process(
+        ["tRNAscan-SE", "-G", "-o", raw_trnas, "--thread", str(threads), fasta],
+        logger,
+        verbose=verbose,
+    )
     if path.isfile(raw_trnas) and stat(raw_trnas).st_size > 0:
-        trna_frame = pd.read_csv(raw_trnas, sep='\t', skiprows=[0, 2])
+        trna_frame = pd.read_csv(raw_trnas, sep="\t", skiprows=[0, 2])
         trna_frame.columns = [i.strip() for i in trna_frame.columns]
         # if begin.1 or end.1 are in trnas then drop, else drop the second begin or end
-        if 'Begin.1' in trna_frame.columns:
-            trna_frame = trna_frame.drop(['Begin.1'], axis=1)
-        if 'End.1' in trna_frame.columns:
-            trna_frame = trna_frame.drop(['End.1'], axis=1)
+        if "Begin.1" in trna_frame.columns:
+            trna_frame = trna_frame.drop(["Begin.1"], axis=1)
+        if "End.1" in trna_frame.columns:
+            trna_frame = trna_frame.drop(["End.1"], axis=1)
         trna_frame = trna_frame.loc[:, get_dups(trna_frame.columns)]
-        trna_frame.insert(0, 'fasta', fasta_name)
+        trna_frame.insert(0, "fasta", fasta_name)
         return trna_frame
     else:
-        logger.warning('No tRNAs were detected, no trnas.tsv file will be created.')
+        logger.warning("No tRNAs were detected, no trnas.tsv file will be created.")
         return None
 
 
-RAW_RRNA_COLUMNS = ['scaffold', 'tool_name', 'type', 'begin', 'end', 'e-value', 'strand', 'empty', 'note']
-RRNA_COLUMNS = ['fasta', 'begin', 'end', 'strand', 'type', 'e-value', 'note']
+RAW_RRNA_COLUMNS = [
+    "scaffold",
+    "tool_name",
+    "type",
+    "begin",
+    "end",
+    "e-value",
+    "strand",
+    "empty",
+    "note",
+]
+RRNA_COLUMNS = ["fasta", "begin", "end", "strand", "type", "e-value", "note"]
 
 
 def run_barrnap(fasta, fasta_name, logger, threads=10, verbose=True):
-    raw_rrna_str = run_process(['barrnap', '--threads', str(threads), fasta], logger, capture_stdout=True, check=False,
-                               verbose=verbose)
-    raw_rrna_table = pd.read_csv(io.StringIO(raw_rrna_str), skiprows=1, sep='\t', header=None,
-                                 names=RAW_RRNA_COLUMNS, index_col=0)
+    raw_rrna_str = run_process(
+        ["barrnap", "--threads", str(threads), fasta],
+        logger,
+        capture_stdout=True,
+        check=False,
+        verbose=verbose,
+    )
+    raw_rrna_table = pd.read_csv(
+        io.StringIO(raw_rrna_str),
+        skiprows=1,
+        sep="\t",
+        header=None,
+        names=RAW_RRNA_COLUMNS,
+        index_col=0,
+    )
     rrna_table_rows = list()
     for gene, row in raw_rrna_table.iterrows():
-        rrna_row_dict = {entry.split('=')[0]: entry.split('=')[1] for entry in row['note'].split(';')}
-        rrna_table_rows.append([fasta_name, row.begin, row.end, row.strand, rrna_row_dict['Name'].replace('_', ' '),
-                                row['e-value'], rrna_row_dict.get('note', '')])
+        rrna_row_dict = {
+            entry.split("=")[0]: entry.split("=")[1] for entry in row["note"].split(";")
+        }
+        rrna_table_rows.append(
+            [
+                fasta_name,
+                row.begin,
+                row.end,
+                row.strand,
+                rrna_row_dict["Name"].replace("_", " "),
+                row["e-value"],
+                rrna_row_dict.get("note", ""),
+            ]
+        )
     if len(raw_rrna_table) > 0:
-        return pd.DataFrame(rrna_table_rows, index=raw_rrna_table.index, columns=RRNA_COLUMNS).reset_index()
+        return pd.DataFrame(
+            rrna_table_rows, index=raw_rrna_table.index, columns=RRNA_COLUMNS
+        ).reset_index()
     else:
-        logger.warning('No rRNAs were detected, no rrnas.tsv file will be created.')
+        logger.warning("No rRNAs were detected, no rrnas.tsv file will be created.")
         return None
 
 
 def make_trnas_interval(scaffold, row, i):
-    if row['Begin'] < row['End']:
-        begin = row['Begin']
-        end = row['End']
-        strand = '+'
+    if row["Begin"] < row["End"]:
+        begin = row["Begin"]
+        end = row["End"]
+        strand = "+"
     else:
-        begin = row['End']
-        end = row['Begin']
-        strand = '-'
-    metadata = {'source': 'tRNAscan-SE', 'type': 'tRNA', 'score': row['Score'], 'strand': strand, 'phase': 0,
-                'ID': '%s_tRNA_%s' % (scaffold, i), 'codon': row['Codon'], 'product': 'tRNA-%s' % row['Type']}
-    if not pd.isna(row['Note']):
-        metadata['Note'] = row['Note']
+        begin = row["End"]
+        end = row["Begin"]
+        strand = "-"
+    metadata = {
+        "source": "tRNAscan-SE",
+        "type": "tRNA",
+        "score": row["Score"],
+        "strand": strand,
+        "phase": 0,
+        "ID": "%s_tRNA_%s" % (scaffold, i),
+        "codon": row["Codon"],
+        "product": "tRNA-%s" % row["Type"],
+    }
+    if not pd.isna(row["Note"]):
+        metadata["Note"] = row["Note"]
     return begin, end, metadata
 
 
 def make_rrnas_interval(scaffold, row, i):
-    metadata = {'source': 'barrnap', 'type': 'rRNA', 'score': row['e-value'], 'strand': row['strand'], 'phase': 0,
-                'ID': '%s_rRNA_%s' % (scaffold, i), 'product': '%s ribosomal RNA' % row['type'].split()[0],
-                'gene': row['type']}
-    if not pd.isna(row['note']):
-        metadata['Note'] = row['note']
-    return row['begin'], row['end'], metadata
+    metadata = {
+        "source": "barrnap",
+        "type": "rRNA",
+        "score": row["e-value"],
+        "strand": row["strand"],
+        "phase": 0,
+        "ID": "%s_rRNA_%s" % (scaffold, i),
+        "product": "%s ribosomal RNA" % row["type"].split()[0],
+        "gene": row["type"],
+    }
+    if not pd.isna(row["note"]):
+        metadata["Note"] = row["note"]
+    return row["begin"], row["end"], metadata
 
 
 # TODO: make it take an input and output gff location and not overwrite
 # TODO: for some reason 1 is getting added to intervals when added to gff
-def add_intervals_to_gff(annotations_loc, gff_loc, len_dict, interval_function, groupby_column, logger):
+def add_intervals_to_gff(
+    annotations_loc, gff_loc, len_dict, interval_function, groupby_column, logger
+):
     # get fasta length dict so we can merge, I'd love to be able to get this from somewhere else
-    annotation_frame = pd.read_csv(annotations_loc, sep='\t')
+    annotation_frame = pd.read_csv(annotations_loc, sep="\t")
     # process trnas to intervals
     annotation_dict = dict()
     for scaffold, frame in annotation_frame.groupby(groupby_column):
@@ -550,58 +782,99 @@ def add_intervals_to_gff(annotations_loc, gff_loc, len_dict, interval_function,
             begin -= 1
             if begin < 0:
                 begin = 0
-                logger.warning('Interval added to gff started less than zero, set to zero')
+                logger.warning(
+                    "Interval added to gff started less than zero, set to zero"
+                )
             im.add(bounds=[(begin, end)], metadata=metadata)
         annotation_dict[scaffold] = im
     # add trna intervals to gff
-    gff = list(read_sequence(gff_loc, format='gff3'))
-    with open(gff_loc, 'w') as f:
+    gff = list(read_sequence(gff_loc, format="gff3"))
+    with open(gff_loc, "w") as f:
         for scaffold, gff_intervals in gff:
             gff_intervals = IntervalMetadata(len_dict[scaffold], gff_intervals)
             if scaffold in annotation_dict:
                 gff_intervals.merge(annotation_dict[scaffold])
                 gff_intervals.sort()
-            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, formater, logger,
-                          db_name='database', bit_score_threshold=60, rbh_bit_score_threshold=350, threads=10,
-                          verbose=False):
+            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,
+    formater,
+    logger,
+    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"""
     # Get kegg hits
-    logger.info('Getting forward best hits from %s' % db_name)
-    forward_hits = get_best_hits(query_db, target_db, logger, working_dir, 'gene', db_name, bit_score_threshold,
-                                 threads, verbose=verbose)
+    logger.info("Getting forward best hits from %s" % db_name)
+    forward_hits = get_best_hits(
+        query_db,
+        target_db,
+        logger,
+        working_dir,
+        "gene",
+        db_name,
+        bit_score_threshold,
+        threads,
+        verbose=verbose,
+    )
     if stat(forward_hits).st_size == 0:
         return pd.DataFrame()
-    logger.info('Getting reverse best hits from %s' % db_name)
-    reverse_hits = get_reciprocal_best_hits(query_db, target_db, logger, working_dir, 'gene', db_name,
-                                            bit_score_threshold, rbh_bit_score_threshold, threads, verbose=verbose)
+    logger.info("Getting reverse best hits from %s" % db_name)
+    reverse_hits = get_reciprocal_best_hits(
+        query_db,
+        target_db,
+        logger,
+        working_dir,
+        "gene",
+        db_name,
+        bit_score_threshold,
+        rbh_bit_score_threshold,
+        threads,
+        verbose=verbose,
+    )
     hits = process_reciprocal_best_hits(forward_hits, reverse_hits, db_name)
-    logger.info('Getting descriptions of hits from %s' % (db_name))
-    if '%s_description' % db_name in db_handler.get_database_names():
-        header_dict = db_handler.get_descriptions(hits['%s_hit' % db_name], '%s_description' % db_name)
+    logger.info("Getting descriptions of hits from %s" % (db_name))
+    if "%s_description" % db_name in db_handler.get_database_names():
+        header_dict = db_handler.get_descriptions(
+            hits["%s_hit" % db_name], "%s_description" % db_name
+        )
     else:
-        header_dict = multigrep(hits[f'{db_name}_hit'], f'{target_db}_h', logger, '\x00', working_dir)
+        header_dict = multigrep(
+            hits[f"{db_name}_hit"], f"{target_db}_h", logger, "\x00", working_dir
+        )
     hits = formater(hits, header_dict)
     return hits
 
 
-def count_motifs(gene_faa, motif='(C..CH)'):
+def count_motifs(gene_faa, motif="(C..CH)"):
     motif_count_dict = dict()
-    for seq in read_sequence(gene_faa, format='fasta'):
-        motif_count_dict[seq.metadata['id']] = len(list(seq.find_with_regex(motif)))
+    for seq in read_sequence(gene_faa, format="fasta"):
+        motif_count_dict[seq.metadata["id"]] = len(list(seq.find_with_regex(motif)))
     return motif_count_dict
 
 
 def strip_endings(text, suffixes: list):
     for suffix in suffixes:
         if text.endswith(suffix):
-            text = text[:len(text) - len(suffix)]
+            text = text[: len(text) - len(suffix)]
     return text
 
 
-def process_custom_dbs(custom_fasta_loc, custom_db_name, output_dir, logger, threads=1, verbose=False):
+def process_custom_dbs(
+    custom_fasta_loc, custom_db_name, output_dir, logger, threads=1, verbose=False
+):
     # if none is passed from argparse then set to tuple of len 0
     mkdir(output_dir)
 
@@ -610,11 +883,15 @@ def process_custom_dbs(custom_fasta_loc, custom_db_name, output_dir, logger, thr
     if custom_db_name is None:
         custom_db_name = ()
     if len(custom_fasta_loc) != len(custom_db_name):
-        raise ValueError('Lengths of custom db fasta list and custom db name list must be the same.')
-    custom_dbs = {custom_db_name[i]: custom_fasta_loc[i] for i in range(len(custom_db_name))}
+        raise ValueError(
+            "Lengths of custom db fasta list and custom db name list must be the same."
+        )
+    custom_dbs = {
+        custom_db_name[i]: custom_fasta_loc[i] for i in range(len(custom_db_name))
+    }
     custom_db_locs = dict()
     for db_name, db_loc in custom_dbs.items():
-        custom_db_loc = path.join(output_dir, '%s.custom.mmsdb' % db_name)
+        custom_db_loc = path.join(output_dir, "%s.custom.mmsdb" % db_name)
         make_mmseqs_db(db_loc, custom_db_loc, logger, threads=threads, verbose=verbose)
         custom_db_locs[db_name] = custom_db_loc
     return custom_db_locs
@@ -626,27 +903,40 @@ def process_custom_hmms(custom_hmm_loc, custom_hmm_name, logger, verbose=False):
     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.')
+        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]], logger, verbose=verbose)  # all are pressed just in case
+        run_process(
+            ["hmmpress", "-f", custom_hmm_loc[i]], logger, 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, logger:logging.Logger,  verbose=False):
+
+def process_custom_hmm_cutoffs(
+    custom_hmm_cutoffs_loc, custom_hmm_name, logger: logging.Logger, 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.")
+        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):
-        logger.warning(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)}
+        logger.warning(
+            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):
+    def __init__(
+        self, name, scaffolds, genes_faa, genes_fna, gff, gbk, annotations, trnas, rrnas
+    ):
         # TODO: get abspath for every input file/dir
         # TODO: check that files exist
         self.name = name
@@ -660,135 +950,238 @@ def __init__(self, name, scaffolds, genes_faa, genes_fna, gff, gbk, annotations,
         self.rrnas_loc = rrnas
 
     def get_annotations(self):
-        return pd.read_csv(self.annotations_loc, index_col=0, sep='\t')
+        return pd.read_csv(self.annotations_loc, index_col=0, sep="\t")
 
     def get_trnas(self):
-        return pd.read_csv(self.trnas_loc, sep='\t')
+        return pd.read_csv(self.trnas_loc, sep="\t")
 
     def get_rrnas(self):
-        return pd.read_csv(self.rrnas_loc, sep='\t')
-
-def annotate_orfs(gene_faa, db_handler, tmp_dir, logger, 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):
+        return pd.read_csv(self.rrnas_loc, sep="\t")
+
+
+def annotate_orfs(
+    gene_faa,
+    db_handler,
+    tmp_dir,
+    logger,
+    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
-    logger.info('Turning genes from prodigal to mmseqs2 db')
-    query_db = path.join(tmp_dir, 'gene.mmsdb')
-    make_mmseqs_db(gene_faa, query_db, logger, create_index=True, threads=threads, verbose=verbose)
+    logger.info("Turning genes from prodigal to mmseqs2 db")
+    query_db = path.join(tmp_dir, "gene.mmsdb")
+    make_mmseqs_db(
+        gene_faa, query_db, logger, create_index=True, threads=threads, verbose=verbose
+    )
 
     annotation_list = list()
 
-    if db_handler.config['search_databases'].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.config['search_databases']['kegg'], tmp_dir,
-                                                     db_handler, get_kegg_description, logger,
-                                                     'kegg', bit_score_threshold, rbh_bit_score_threshold, threads,
-                                                     verbose))
-    elif db_handler.config['search_databases'].get('kofam_hmm') is not None and db_handler.config['search_databases'].get('kofam_ko_list') is not None:
-        logger.info('Getting hits from kofam')
-        annotation_list.append(run_hmmscan(genes_faa=gene_faa,
-                                           db_loc=db_handler.config['search_databases']['kofam_hmm'],
-                                           db_name='kofam_hmm',
-                                           output_loc=tmp_dir, #check_impliments
-                                           threads=threads, #check_impliments
-                                           verbose=verbose,
-                                           formater=partial(
-                                               kofam_hmmscan_formater,
-                                               hmm_info_path=db_handler.config['search_databases']['kofam_ko_list'],
-                                               top_hit=True,
-                                               use_dbcan2_thresholds=kofam_use_dbcan2_thresholds
-                                           ),
-                                           logger=logger))
+    if db_handler.config["search_databases"].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.config["search_databases"]["kegg"],
+                tmp_dir,
+                db_handler,
+                get_kegg_description,
+                logger,
+                "kegg",
+                bit_score_threshold,
+                rbh_bit_score_threshold,
+                threads,
+                verbose,
+            )
+        )
+    elif (
+        db_handler.config["search_databases"].get("kofam_hmm") is not None
+        and db_handler.config["search_databases"].get("kofam_ko_list") is not None
+    ):
+        logger.info("Getting hits from kofam")
+        annotation_list.append(
+            run_hmmscan(
+                genes_faa=gene_faa,
+                db_loc=db_handler.config["search_databases"]["kofam_hmm"],
+                db_name="kofam_hmm",
+                output_loc=tmp_dir,  # check_impliments
+                threads=threads,  # check_impliments
+                verbose=verbose,
+                formater=partial(
+                    kofam_hmmscan_formater,
+                    hmm_info_path=db_handler.config["search_databases"][
+                        "kofam_ko_list"
+                    ],
+                    top_hit=True,
+                    use_dbcan2_thresholds=kofam_use_dbcan2_thresholds,
+                ),
+                logger=logger,
+            )
+        )
     else:
-        logger.warning('No KEGG source provided so distillation will be of limited use.')
+        logger.warning(
+            "No KEGG source provided so distillation will be of limited use."
+        )
 
     # Get uniref hits
-    if db_handler.config['search_databases'].get('uniref') is not None:
-        annotation_list.append(do_blast_style_search(query_db, db_handler.config['search_databases']['uniref'], tmp_dir,
-                                                     db_handler, get_uniref_description,
-                                                     logger, 'uniref', bit_score_threshold,
-                                                     rbh_bit_score_threshold, threads, verbose))
+    if db_handler.config["search_databases"].get("uniref") is not None:
+        annotation_list.append(
+            do_blast_style_search(
+                query_db,
+                db_handler.config["search_databases"]["uniref"],
+                tmp_dir,
+                db_handler,
+                get_uniref_description,
+                logger,
+                "uniref",
+                bit_score_threshold,
+                rbh_bit_score_threshold,
+                threads,
+                verbose,
+            )
+        )
 
     # Get viral hits
-    if db_handler.config['search_databases'].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_handler.config['search_databases']['viral'], tmp_dir,
-                                                     db_handler, get_viral_description,
-                                                     logger, 'viral', bit_score_threshold,
-                                                     rbh_bit_score_threshold, threads, verbose))
+    if db_handler.config["search_databases"].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_handler.config["search_databases"]["viral"],
+                tmp_dir,
+                db_handler,
+                get_viral_description,
+                logger,
+                "viral",
+                bit_score_threshold,
+                rbh_bit_score_threshold,
+                threads,
+                verbose,
+            )
+        )
 
     # Get peptidase hits
-    if db_handler.config['search_databases'].get('peptidase') is not None:
-        annotation_list.append(do_blast_style_search(query_db, db_handler.config['search_databases']['peptidase'], tmp_dir,
-                                                     db_handler, get_peptidase_description,
-                                                     logger, 'peptidase', bit_score_threshold,
-                                                     rbh_bit_score_threshold, threads, verbose))
+    if db_handler.config["search_databases"].get("peptidase") is not None:
+        annotation_list.append(
+            do_blast_style_search(
+                query_db,
+                db_handler.config["search_databases"]["peptidase"],
+                tmp_dir,
+                db_handler,
+                get_peptidase_description,
+                logger,
+                "peptidase",
+                bit_score_threshold,
+                rbh_bit_score_threshold,
+                threads,
+                verbose,
+            )
+        )
 
     # Get pfam hits
-    if db_handler.config['search_databases'].get('pfam') is not None:
-        logger.info('Getting hits from pfam')
-        annotation_list.append(run_mmseqs_profile_search(query_db, db_handler.config['search_databases']['pfam'], tmp_dir,
-                                                         logger, output_prefix='pfam', db_handler=db_handler, threads=threads,
-                                                         verbose=verbose))
+    if db_handler.config["search_databases"].get("pfam") is not None:
+        logger.info("Getting hits from pfam")
+        annotation_list.append(
+            run_mmseqs_profile_search(
+                query_db,
+                db_handler.config["search_databases"]["pfam"],
+                tmp_dir,
+                logger,
+                output_prefix="pfam",
+                db_handler=db_handler,
+                threads=threads,
+                verbose=verbose,
+            )
+        )
 
     # use hmmer to detect cazy ids using dbCAN
-    if db_handler.config['search_databases'].get('dbcan') is not None:
-        logger.info('Getting hits from dbCAN')
-        annotation_list.append(run_hmmscan(genes_faa=gene_faa,
-                                           db_loc=db_handler.config['search_databases']['dbcan'],
-                                           db_name='cazy',
-                                           output_loc=tmp_dir,
-                                           threads=threads,
-                                           formater=partial(
-                                               dbcan_hmmscan_formater,
-                                               db_name='cazy',
-                                               db_handler=db_handler
-                                           ),
-                                           logger=logger))
+    if db_handler.config["search_databases"].get("dbcan") is not None:
+        logger.info("Getting hits from dbCAN")
+        annotation_list.append(
+            run_hmmscan(
+                genes_faa=gene_faa,
+                db_loc=db_handler.config["search_databases"]["dbcan"],
+                db_name="cazy",
+                output_loc=tmp_dir,
+                threads=threads,
+                formater=partial(
+                    dbcan_hmmscan_formater, db_name="cazy", db_handler=db_handler
+                ),
+                logger=logger,
+            )
+        )
 
     # use hmmer to detect vogdbs
-    if db_handler.config['search_databases'].get('vogdb') is not None:
-        logger.info('Getting hits from VOGDB')
-        annotation_list.append(run_hmmscan(genes_faa=gene_faa,
-                                           db_loc=db_handler.config['search_databases']['vogdb'],
-                                           db_name='vogdb',
-                                           threads=threads,
-                                           output_loc=tmp_dir,
-                                           formater=partial(
-                                               vogdb_hmmscan_formater,
-                                               db_name='vogdb',
-                                               logger=logger,
-                                               db_handler=db_handler
-                                           ),
-                                           logger=logger))
+    if db_handler.config["search_databases"].get("vogdb") is not None:
+        logger.info("Getting hits from VOGDB")
+        annotation_list.append(
+            run_hmmscan(
+                genes_faa=gene_faa,
+                db_loc=db_handler.config["search_databases"]["vogdb"],
+                db_name="vogdb",
+                threads=threads,
+                output_loc=tmp_dir,
+                formater=partial(
+                    vogdb_hmmscan_formater,
+                    db_name="vogdb",
+                    logger=logger,
+                    db_handler=db_handler,
+                ),
+                logger=logger,
+            )
+        )
     for db_name, db_loc in custom_db_locs.items():
-        logger.info('Getting hits from %s' % db_name)
+        logger.info("Getting hits from %s" % db_name)
         get_custom_description = partial(get_basic_description, db_name=db_name)
-        annotation_list.append(do_blast_style_search(query_db, db_loc, tmp_dir, db_handler,
-                                                     get_custom_description, logger, db_name,
-                                                     bit_score_threshold, rbh_bit_score_threshold, threads,
-                                                     verbose))
+        annotation_list.append(
+            do_blast_style_search(
+                query_db,
+                db_loc,
+                tmp_dir,
+                db_handler,
+                get_custom_description,
+                logger,
+                db_name,
+                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
-                                           ),
-                                           logger=logger))
-
-    annotation_list.append(pd.DataFrame(count_motifs(gene_faa, '(C..CH)'), index=['heme_regulatory_motif_count']).T)
+        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,
+                ),
+                logger=logger,
+            )
+        )
+
+    annotation_list.append(
+        pd.DataFrame(
+            count_motifs(gene_faa, "(C..CH)"), index=["heme_regulatory_motif_count"]
+        ).T
+    )
 
     # merge dataframes
-    logger.info('Merging ORF annotations')
+    logger.info("Merging ORF annotations")
     annotations = pd.concat(annotation_list, axis=1, sort=False)
 
     # get scaffold data and assign grades
@@ -798,35 +1191,68 @@ def annotate_orfs(gene_faa, db_handler, tmp_dir, logger, custom_db_locs=(), cust
     return annotations
 
 
-def annotate_fasta(fasta_loc, fasta_name, output_dir, db_handler, logger, 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, threads=1, rename_bins=True, keep_tmp_dir=False,
-                   verbose=False):
+def annotate_fasta(
+    fasta_loc,
+    fasta_name,
+    output_dir,
+    db_handler,
+    logger,
+    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,
+    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')
+    tmp_dir = path.join(output_dir, "tmp")
     mkdir(tmp_dir)
 
     # filter input fasta
-    filtered_fasta = path.join(tmp_dir, 'filtered_fasta.fa')
+    filtered_fasta = path.join(tmp_dir, "filtered_fasta.fa")
     filter_fasta(fasta_loc, min_contig_size, filtered_fasta)
 
     if stat(filtered_fasta).st_size == 0:
-        logger.warning('No sequences were longer than min_contig_size')
+        logger.warning("No sequences were longer than min_contig_size")
         return None
     # predict ORFs with prodigal
     # TODO: handle when prodigal returns no genes
-    gene_gff, gene_fna, gene_faa = run_prodigal(filtered_fasta, tmp_dir, logger, mode=prodigal_mode,
-                                                trans_table=trans_table, verbose=verbose)
+    gene_gff, gene_fna, gene_faa = run_prodigal(
+        filtered_fasta,
+        tmp_dir,
+        logger,
+        mode=prodigal_mode,
+        trans_table=trans_table,
+        verbose=verbose,
+    )
 
     # annotate ORFs
-    annotations = annotate_orfs(gene_faa, db_handler, tmp_dir, logger, custom_db_locs, custom_hmm_locs,
-                                custom_hmm_cutoffs_locs, bit_score_threshold, rbh_bit_score_threshold,
-                                kofam_use_dbcan2_thresholds, threads, verbose)
+    annotations = annotate_orfs(
+        gene_faa,
+        db_handler,
+        tmp_dir,
+        logger,
+        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')
+    renamed_scaffolds = path.join(output_dir, "scaffolds.annotated.fa")
     if rename_bins:
         # rename scaffolds to match prodigal names
         prefix = fasta_name
@@ -836,88 +1262,158 @@ def annotate_fasta(fasta_loc, fasta_name, output_dir, db_handler, logger, min_co
         copy2(filtered_fasta, renamed_scaffolds)
 
     # generate fna and faa output files with annotations
-    annotated_fna = path.join(output_dir, 'genes.annotated.fna')
+    annotated_fna = path.join(output_dir, "genes.annotated.fna")
     create_annotated_fasta(gene_fna, annotations, annotated_fna, name=prefix)
-    annotated_faa = path.join(output_dir, 'genes.annotated.faa')
+    annotated_faa = path.join(output_dir, "genes.annotated.faa")
     create_annotated_fasta(gene_faa, annotations, annotated_faa, name=prefix)
 
     # rename gff entries to match prodigal names
-    renamed_gffs = path.join(output_dir, 'genes.annotated.gff')
+    renamed_gffs = path.join(output_dir, "genes.annotated.gff")
     annotate_gff(gene_gff, renamed_gffs, annotations, prefix=prefix)
 
     # add fasta name to frame and index, append to list
-    annotations.insert(0, 'fasta', fasta_name)
+    annotations.insert(0, "fasta", fasta_name)
     if rename_bins:
-        annotations.index = annotations.fasta + '_' + annotations.index
-    annotations_loc = path.join(output_dir, 'annotations.tsv')
-    annotations.to_csv(annotations_loc, sep='\t')
+        annotations.index = annotations.fasta + "_" + annotations.index
+    annotations_loc = path.join(output_dir, "annotations.tsv")
+    annotations.to_csv(annotations_loc, sep="\t")
 
     # get tRNAs and rRNAs
-    len_dict = {i.metadata['id']: len(i) for i in read_sequence(renamed_scaffolds, format='fasta')}
+    len_dict = {
+        i.metadata["id"]: len(i)
+        for i in read_sequence(renamed_scaffolds, format="fasta")
+    }
     if not skip_trnascan:
-        trna_table = run_trna_scan(renamed_scaffolds, tmp_dir, fasta_name, logger, threads=threads, verbose=verbose)
+        trna_table = run_trna_scan(
+            renamed_scaffolds,
+            tmp_dir,
+            fasta_name,
+            logger,
+            threads=threads,
+            verbose=verbose,
+        )
         if trna_table is not None:
-            trna_loc = path.join(output_dir, 'trnas.tsv')
-            trna_table.to_csv(trna_loc, sep='\t', index=False)
-            add_intervals_to_gff(trna_loc, renamed_gffs, len_dict, make_trnas_interval, 'Name', logger)
+            trna_loc = path.join(output_dir, "trnas.tsv")
+            trna_table.to_csv(trna_loc, sep="\t", index=False)
+            add_intervals_to_gff(
+                trna_loc, renamed_gffs, len_dict, make_trnas_interval, "Name", logger
+            )
         else:
             trna_loc = None
     else:
         trna_loc = None
 
-    rrna_table = run_barrnap(renamed_scaffolds, fasta_name, logger,  threads=threads, verbose=verbose)
+    rrna_table = run_barrnap(
+        renamed_scaffolds, fasta_name, logger, threads=threads, verbose=verbose
+    )
     if rrna_table is not None:
-        rrna_loc = path.join(output_dir, 'rrnas.tsv')
-        rrna_table.to_csv(rrna_loc, sep='\t', index=False)
-        add_intervals_to_gff(rrna_loc, renamed_gffs, len_dict, make_rrnas_interval, 'scaffold', logger)
+        rrna_loc = path.join(output_dir, "rrnas.tsv")
+        rrna_table.to_csv(rrna_loc, sep="\t", index=False)
+        add_intervals_to_gff(
+            rrna_loc, renamed_gffs, len_dict, make_rrnas_interval, "scaffold", logger
+        )
     else:
         rrna_loc = None
 
     # make genbank file
-    current_gbk = path.join(output_dir, '%s.gbk' % fasta_name)
-    make_gbk_from_gff_and_fasta(renamed_gffs, renamed_scaffolds, annotated_faa, current_gbk)
+    current_gbk = path.join(output_dir, "%s.gbk" % fasta_name)
+    make_gbk_from_gff_and_fasta(
+        renamed_gffs, renamed_scaffolds, annotated_faa, current_gbk
+    )
 
     if not keep_tmp_dir:
         rmtree(tmp_dir)
 
-    return Annotation(name=fasta_name, scaffolds=renamed_scaffolds, genes_faa=annotated_faa, genes_fna=annotated_fna,
-                      gff=renamed_gffs, gbk=current_gbk, annotations=annotations_loc, trnas=trna_loc, rrnas=rrna_loc)
+    return Annotation(
+        name=fasta_name,
+        scaffolds=renamed_scaffolds,
+        genes_faa=annotated_faa,
+        genes_fna=annotated_fna,
+        gff=renamed_gffs,
+        gbk=current_gbk,
+        annotations=annotations_loc,
+        trnas=trna_loc,
+        rrnas=rrna_loc,
+    )
 
 
 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_handler, logger, 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,
-                    threads=10, verbose=True):
+    return path.splitext(path.basename(remove_suffix(fasta_loc, ".gz")))[0]
+
+
+def annotate_fastas(
+    fasta_locs,
+    output_dir,
+    db_handler,
+    logger,
+    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,
+    threads=10,
+    verbose=True,
+):
     # check for no conflicting options/configurations
-    tmp_dir = path.join(output_dir, 'working_dir')
+    tmp_dir = path.join(output_dir, "working_dir")
     mkdir(tmp_dir)
 
     # setup custom databases to be searched
-    custom_db_locs = process_custom_dbs(custom_fasta_loc, custom_db_name, path.join(tmp_dir, 'custom_dbs'),
-                                        logger, threads, verbose)
+    custom_db_locs = process_custom_dbs(
+        custom_fasta_loc,
+        custom_db_name,
+        path.join(tmp_dir, "custom_dbs"),
+        logger,
+        threads,
+        verbose,
+    )
     custom_hmm_locs = process_custom_hmms(custom_hmm_loc, custom_hmm_name, logger)
-    custom_hmm_cutoffs_locs= process_custom_hmm_cutoffs(custom_hmm_cutoffs_loc, custom_hmm_name, logger)
-    logger.info('Retrieved database locations and descriptions')
+    custom_hmm_cutoffs_locs = process_custom_hmm_cutoffs(
+        custom_hmm_cutoffs_loc, custom_hmm_name, logger
+    )
+    logger.info("Retrieved database locations and descriptions")
 
     # iterate over list of fastas and annotate each individually
     annotations_list = list()
     for fasta_loc in fasta_locs:
         # get name of file e.g. /home/shaffemi/my_genome.fa -> my_genome
         fasta_name = get_fasta_name(fasta_loc)
-        logger.info('Annotating %s' % fasta_name)
+        logger.info("Annotating %s" % 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_handler, logger, 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, threads, rename_bins, keep_tmp_dir, verbose))
-    logger.info('Annotations complete, processing annotations')
+            annotate_fasta(
+                fasta_loc,
+                fasta_name,
+                fasta_dir,
+                db_handler,
+                logger,
+                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,
+                threads,
+                rename_bins,
+                keep_tmp_dir,
+                verbose,
+            )
+        )
+    logger.info("Annotations complete, processing annotations")
 
     all_annotations = merge_annotations(annotations_list, output_dir)
 
@@ -929,82 +1425,140 @@ def annotate_fastas(fasta_locs, output_dir, db_handler, logger, min_contig_size=
 
 def make_fasta_namses_df(fasta_loc):
     fasta_name = get_fasta_name(fasta_loc)
-    names = [{'fasta': fasta_name, 'seq': seq.metadata['id']} for seq in read_sequence(fasta_loc,format='fasta')]
+    names = [
+        {"fasta": fasta_name, "seq": seq.metadata["id"]}
+        for seq in read_sequence(fasta_loc, format="fasta")
+    ]
     return pd.DataFrame(names)
 
 
 # TODO: Add force flag to remove output dir if it already exists
 # TODO: Add continute flag to continue if output directory already exists
 # TODO: make fasta loc either a string or list to remove annotate_bins_cmd and annotate_called_genes_cmd?
-def annotate_bins(input_fasta:list, 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=(),
-                  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,
-                  log_file_path:str=None, join_fastas:bool=False, config_loc:str=None):
+def annotate_bins(
+    input_fasta: list,
+    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=(),
+    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,
+    log_file_path: str = None,
+    join_fastas: bool = False,
+    config_loc: str = None,
+):
     rename_bins = True
     fasta_locs = [j for i in input_fasta for j in glob(i)]
     mkdir(output_dir)
     if log_file_path is None:
         log_file_path = path.join(output_dir, "annotate.log")
-    logger = logging.getLogger('annotation_log')
+    logger = logging.getLogger("annotation_log")
     setup_logger(logger, log_file_path)
     logger.info(f"The log file is created at {log_file_path}.")
 
     if len(fasta_locs) == 0:
-        raise ValueError('Given fasta locations return no paths: %s' % input_fasta)
+        raise ValueError("Given fasta locations return no paths: %s" % input_fasta)
     fasta_names = [get_fasta_name(i) for i in fasta_locs]
     if len(fasta_names) != len(set(fasta_names)):
-        raise ValueError('Genome file names must be unique. At least one name appears twice in this search.')
-    logger.info('%s FASTAs found' % len(fasta_locs))
+        raise ValueError(
+            "Genome file names must be unique. At least one name appears twice in this search."
+        )
+    logger.info("%s FASTAs found" % len(fasta_locs))
     # set up
     db_handler = DatabaseHandler(logger, config_loc)
-    db_handler.filter_db_locs(low_mem_mode, use_uniref,
-                              use_vogdb, master_list=MAG_DBS_TO_ANNOTATE,)
+    db_handler.filter_db_locs(
+        low_mem_mode,
+        use_uniref,
+        use_vogdb,
+        master_list=MAG_DBS_TO_ANNOTATE,
+    )
     db_conf = db_handler.get_settings_str()
-    logger.info(f"Starting the Annotation of Bins with database configuration: \n {db_conf}")
+    logger.info(
+        f"Starting the Annotation of Bins with database configuration: \n {db_conf}"
+    )
 
     # check inputs
-    prodigal_modes = ['train', 'meta', 'single']
+    prodigal_modes = ["train", "meta", "single"]
     if prodigal_mode not in prodigal_modes:
-        raise ValueError('Prodigal mode must be one of %s.' % ', '.join(prodigal_modes))
-    elif prodigal_mode in ['normal', 'single']:
-        logger.warning('When running prodigal in single mode your bins must have long contigs (average length >3 Kbp), '
-                      'be long enough (total length > 500 Kbp) and have very low contamination in order for prodigal '
-                      'training to work well.')
+        raise ValueError("Prodigal mode must be one of %s." % ", ".join(prodigal_modes))
+    elif prodigal_mode in ["normal", "single"]:
+        logger.warning(
+            "When running prodigal in single mode your bins must have long contigs (average length >3 Kbp), "
+            "be long enough (total length > 500 Kbp) and have very low contamination in order for prodigal "
+            "training to work well."
+        )
 
     # prodigal_trans_tables = ['auto'] + [str(i) for i in range(1, 26)]
     prodigal_trans_tables = [str(i) for i in range(1, 26)]
     if trans_table not in prodigal_trans_tables:
         # raise ValueError('Prodigal translation table must be 1-25 or auto')
-        raise ValueError('Prodigal translation table must be 1-25')
-
-
-    all_annotations = annotate_fastas(fasta_locs, output_dir, db_handler, logger, 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,
-                                      threads, verbose)
+        raise ValueError("Prodigal translation table must be 1-25")
+
+    all_annotations = annotate_fastas(
+        fasta_locs,
+        output_dir,
+        db_handler,
+        logger,
+        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,
+        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])
+        gtdb_taxonomy = pd.concat(
+            [pd.read_csv(i, sep="\t", index_col=0) for i in gtdb_taxonomy]
+        )
         taxonomy = list()
         taxonomy_missing_bins = list()
         for i in all_annotations.fasta:
             # add taxonomy
             if i in gtdb_taxonomy.index:
-                taxonomy.append(gtdb_taxonomy.loc[i, 'classification'])
+                taxonomy.append(gtdb_taxonomy.loc[i, "classification"])
             else:
                 taxonomy.append(i)
                 taxonomy_missing_bins.append(i)
         for i in set(taxonomy_missing_bins):
-            logger.warning('Bin %s was not found in taxonomy file, replaced with bin name.' % i)
-        all_annotations['bin_taxonomy'] = taxonomy
+            logger.warning(
+                "Bin %s was not found in taxonomy file, replaced with bin name." % i
+            )
+        all_annotations["bin_taxonomy"] = taxonomy
     # if given add quality information
     if len(checkm_quality) > 0:
-        checkm_quality = pd.concat([pd.read_csv(i, sep='\t', index_col=0) for i in checkm_quality])
-        checkm_quality.index = [strip_endings(i, ['.fa', '.fasta', '.fna']) for i in checkm_quality.index]
+        checkm_quality = pd.concat(
+            [pd.read_csv(i, sep="\t", index_col=0) for i in checkm_quality]
+        )
+        checkm_quality.index = [
+            strip_endings(i, [".fa", ".fasta", ".fna"]) for i in checkm_quality.index
+        ]
 
         completeness = list()
         contamination = list()
@@ -1012,90 +1566,155 @@ def annotate_bins(input_fasta:list, output_dir='.', min_contig_size=2500, prodig
         for i in all_annotations.fasta:
             # add completeness and contamination
             if i in checkm_quality.index:
-                completeness.append(checkm_quality.loc[i, 'Completeness'])
-                contamination.append(checkm_quality.loc[i, 'Contamination'])
+                completeness.append(checkm_quality.loc[i, "Completeness"])
+                contamination.append(checkm_quality.loc[i, "Contamination"])
             else:
                 completeness.append(0)
                 contamination.append(100)
                 quality_missing_bins.append(i)
             for j in set(quality_missing_bins):
-                logger.warning('Bin %s was not found in quality file, '
-                              'replaced with completeness 0 and contamination 100.' % j)
-        all_annotations['bin_completeness'] = completeness
-        all_annotations['bin_contamination'] = contamination
-    all_annotations.to_csv(path.join(output_dir, 'annotations.tsv'), sep='\t')
+                logger.warning(
+                    "Bin %s was not found in quality file, "
+                    "replaced with completeness 0 and contamination 100." % j
+                )
+        all_annotations["bin_completeness"] = completeness
+        all_annotations["bin_contamination"] = contamination
+    all_annotations.to_csv(path.join(output_dir, "annotations.tsv"), sep="\t")
 
     logger.info("Completed annotations")
 
 
-def annotate_called_genes_cmd(input_faa, output_dir='.', bit_score_threshold=60,
-                              rbh_bit_score_threshold=350,
-                              custom_db_name=(), custom_fasta_loc=(), custom_hmm_loc=(), custom_hmm_name=(),
-                              custom_hmm_cutoffs_loc=(), use_uniref=False, log_file_path:str=None,
-                              use_vogdb=False, kofam_use_dbcan2_thresholds=False, rename_genes=True,
-                              keep_tmp_dir=True, low_mem_mode=False, threads=10, verbose=True,
-                              config_loc:str=None):
+def annotate_called_genes_cmd(
+    input_faa,
+    output_dir=".",
+    bit_score_threshold=60,
+    rbh_bit_score_threshold=350,
+    custom_db_name=(),
+    custom_fasta_loc=(),
+    custom_hmm_loc=(),
+    custom_hmm_name=(),
+    custom_hmm_cutoffs_loc=(),
+    use_uniref=False,
+    log_file_path: str = None,
+    use_vogdb=False,
+    kofam_use_dbcan2_thresholds=False,
+    rename_genes=True,
+    keep_tmp_dir=True,
+    low_mem_mode=False,
+    threads=10,
+    verbose=True,
+    config_loc: str = None,
+):
     fasta_locs = glob(input_faa)
-    annotate_called_genes(fasta_locs, output_dir, bit_score_threshold, rbh_bit_score_threshold,
-                          custom_db_name, 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, None,
-                          config_loc)
+    annotate_called_genes(
+        fasta_locs,
+        output_dir,
+        bit_score_threshold,
+        rbh_bit_score_threshold,
+        custom_db_name,
+        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,
+        None,
+        config_loc,
+    )
+
 
 def perform_fasta_checks(fasta_locs, logger):
     """Perform all checks related to integrity of fastas"""
     if len(fasta_locs) == 0:
-        raise ValueError('Given fasta locations returns no paths: %s' % input_faa)
-    logger.info('%s fastas found' % len(fasta_locs))
+        raise ValueError("Given fasta locations returns no paths: %s" % input_faa)
+    logger.info("%s fastas found" % len(fasta_locs))
     logger.info("Checking for duplicate names")
     names = [get_fasta_name(i) for i in fasta_locs]
     try:
         if len(names) != len(set(names)):
-            raise ValueError('Genome file names must be unique. At least one name appears twice in this search.')
-        fastas_dup_check(fasta_locs, '>')
+            raise ValueError(
+                "Genome file names must be unique. At least one name appears twice in this search."
+            )
+        fastas_dup_check(fasta_locs, ">")
     except ValueError as error:
         logger.critical(error)
         raise error
     logger.info("No duplicate names found")
 
-def annotate_called_genes(fasta_locs, output_dir='.', bit_score_threshold=60, rbh_bit_score_threshold=350,
-                          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, log_file_path:str=None, config_loc:str=None):
+
+def annotate_called_genes(
+    fasta_locs,
+    output_dir=".",
+    bit_score_threshold=60,
+    rbh_bit_score_threshold=350,
+    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,
+    log_file_path: str = None,
+    config_loc: str = None,
+):
     mkdir(output_dir)
 
     # Get a logger
     if log_file_path is None:
         log_file_path = path.join(output_dir, "annotate.log")
-    logger = logging.getLogger('annotation_log')
+    logger = logging.getLogger("annotation_log")
     setup_logger(logger, log_file_path)
     logger.info(f"The log file is created at {log_file_path}")
     # get database locations
     db_handler = DatabaseHandler(logger, config_loc)
-    db_handler.filter_db_locs(low_mem_mode, use_uniref, use_vogdb, master_list=MAG_DBS_TO_ANNOTATE)
+    db_handler.filter_db_locs(
+        low_mem_mode, use_uniref, use_vogdb, master_list=MAG_DBS_TO_ANNOTATE
+    )
 
     # Check fastas
     perform_fasta_checks(fasta_locs, logger)
 
-    logger.info(f"Starting the Annotation of Genes with database configuration: \n {db_handler.get_settings_str()}")
+    logger.info(
+        f"Starting the Annotation of Genes with database configuration: \n {db_handler.get_settings_str()}"
+    )
 
     logger.info("Checking for duplicate names")
     fasta_names = [get_fasta_name(i) for i in fasta_locs]
     if len(fasta_names) != len(set(fasta_names)):
-        raise ValueError('Genome file names must be unique. At least one name appears twice in this search.')
+        raise ValueError(
+            "Genome file names must be unique. At least one name appears twice in this search."
+        )
     logger.info("No duplicate names found")
 
-
-    tmp_dir = path.join(output_dir, 'working_dir')
+    tmp_dir = path.join(output_dir, "working_dir")
     mkdir(tmp_dir)
 
     # setup custom databases to be searched
-    custom_db_locs = process_custom_dbs(custom_fasta_loc, custom_db_name, path.join(tmp_dir, 'custom_dbs'),
-                                        logger, threads, verbose)
+    custom_db_locs = process_custom_dbs(
+        custom_fasta_loc,
+        custom_db_name,
+        path.join(tmp_dir, "custom_dbs"),
+        logger,
+        threads,
+        verbose,
+    )
     custom_hmm_locs = process_custom_hmms(custom_hmm_loc, custom_hmm_name, logger)
-    custom_hmm_cutoffs_locs= process_custom_hmm_cutoffs(custom_hmm_cutoffs_loc, custom_hmm_name, logger)
-    logger.info('Retrieved database locations and descriptions')
+    custom_hmm_cutoffs_locs = process_custom_hmm_cutoffs(
+        custom_hmm_cutoffs_loc, custom_hmm_name, logger
+    )
+    logger.info("Retrieved database locations and descriptions")
 
     # annotate
     annotation_locs = list()
@@ -1107,28 +1726,41 @@ def annotate_called_genes(fasta_locs, output_dir='.', bit_score_threshold=60, rb
         mkdir(fasta_dir)
 
         # annotate
-        annotations = annotate_orfs(fasta_loc, db_handler, fasta_dir, logger, 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')
+        annotations = annotate_orfs(
+            fasta_loc,
+            db_handler,
+            fasta_dir,
+            logger,
+            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)
         faa_locs.append(annotated_faa)
 
         # add fasta name to frame and index, write file
-        annotations.insert(0, 'fasta', fasta_name)
+        annotations.insert(0, "fasta", fasta_name)
         if rename_genes:
-            annotations.index = annotations.fasta + '_' + annotations.index
-        annotation_loc = path.join(fasta_dir, 'annotations.tsv')
-        annotations.to_csv(annotation_loc, sep='\t')
+            annotations.index = annotations.fasta + "_" + annotations.index
+        annotation_loc = path.join(fasta_dir, "annotations.tsv")
+        annotations.to_csv(annotation_loc, sep="\t")
         annotation_locs.append(annotation_loc)
 
     # merge
-    all_annotations = pd.concat([pd.read_csv(i, sep='\t', index_col=0) for i in annotation_locs], sort=False)
-    all_annotations = all_annotations.sort_values('fasta')
-    all_annotations.to_csv(path.join(output_dir, 'annotations.tsv'), sep='\t')
-    merge_files(faa_locs, path.join(output_dir, 'genes.faa'))
+    all_annotations = pd.concat(
+        [pd.read_csv(i, sep="\t", index_col=0) for i in annotation_locs], sort=False
+    )
+    all_annotations = all_annotations.sort_values("fasta")
+    all_annotations.to_csv(path.join(output_dir, "annotations.tsv"), sep="\t")
+    merge_files(faa_locs, path.join(output_dir, "genes.faa"))
 
     # clean up
     if not keep_tmp_dir:
@@ -1140,34 +1772,61 @@ def annotate_called_genes(fasta_locs, output_dir='.', bit_score_threshold=60, rb
 def merge_annotations(annotations_list, output_dir, write_annotations=False):
 
     # merge annotation dicts
-    all_annotations = pd.concat([i.get_annotations() for i in annotations_list if i is not None], sort=False)
-    all_annotations = all_annotations.sort_values(['fasta', 'scaffold', 'gene_position'])
+    all_annotations = pd.concat(
+        [i.get_annotations() for i in annotations_list if i is not None], sort=False
+    )
+    all_annotations = all_annotations.sort_values(
+        ["fasta", "scaffold", "gene_position"]
+    )
 
     # merge gene files
-    merge_files([i.genes_fna_loc for i in annotations_list if i is not None], path.join(output_dir, 'genes.fna'))
-    merge_files([i.genes_faa_loc for i in annotations_list if i is not None], path.join(output_dir, 'genes.faa'))
-    merge_files([i.scaffolds_loc for i in annotations_list if i is not None], path.join(output_dir, 'scaffolds.fna'))
-    merge_files([i.gff_loc for i in annotations_list if i is not None], path.join(output_dir, 'genes.gff'), True)
-    trnas_locs = [i.trnas_loc for i in annotations_list if i is not None if i.trnas_loc is not None]
+    merge_files(
+        [i.genes_fna_loc for i in annotations_list if i is not None],
+        path.join(output_dir, "genes.fna"),
+    )
+    merge_files(
+        [i.genes_faa_loc for i in annotations_list if i is not None],
+        path.join(output_dir, "genes.faa"),
+    )
+    merge_files(
+        [i.scaffolds_loc for i in annotations_list if i is not None],
+        path.join(output_dir, "scaffolds.fna"),
+    )
+    merge_files(
+        [i.gff_loc for i in annotations_list if i is not None],
+        path.join(output_dir, "genes.gff"),
+        True,
+    )
+    trnas_locs = [
+        i.trnas_loc
+        for i in annotations_list
+        if i is not None
+        if i.trnas_loc is not None
+    ]
     if len(trnas_locs) > 0:
-        merge_files(trnas_locs, path.join(output_dir, 'trnas.tsv'), True)
-    rrnas_locs = [i.rrnas_loc for i in annotations_list if i is not None if i.rrnas_loc is not None]
+        merge_files(trnas_locs, path.join(output_dir, "trnas.tsv"), True)
+    rrnas_locs = [
+        i.rrnas_loc
+        for i in annotations_list
+        if i is not None
+        if i.rrnas_loc is not None
+    ]
     if len(rrnas_locs) > 0:
-        merge_files(rrnas_locs, path.join(output_dir, 'rrnas.tsv'), True)
+        merge_files(rrnas_locs, path.join(output_dir, "rrnas.tsv"), True)
 
     # make output gbk dir
-    gbk_dir = path.join(output_dir, 'genbank')
+    gbk_dir = path.join(output_dir, "genbank")
     mkdir(gbk_dir)
     for anno in annotations_list:
         if anno is not None:
             # TODO: make annotate_fasta generate a genbank dir and then copy it's contents, get rid of Annotation.name
             if path.isfile(anno.gbk_loc):
-                copy2(anno.gbk_loc, path.join(gbk_dir, '%s.gbk' % anno.name))
+                copy2(anno.gbk_loc, path.join(gbk_dir, "%s.gbk" % anno.name))
             else:
-                for gbk_loc in glob(path.join(anno.gbk_loc, '*.gbk')):
+                for gbk_loc in glob(path.join(anno.gbk_loc, "*.gbk")):
                     copy2(gbk_loc, gbk_dir)
     if write_annotations:
-        all_annotations.to_csv(path.join(output_dir, 'annotations.tsv'), sep='\t')
+        all_annotations.to_csv(path.join(output_dir, "annotations.tsv"), sep="\t")
     else:
         return all_annotations
 
@@ -1178,33 +1837,45 @@ def merge_annotations_cmd(input_dirs, output_dir):
     # Get a logger
     annotations_list = list()
     log_file_path = path.join(output_dir, "annotate.log")
-    logger = logging.getLogger('annotate.log')
+    logger = logging.getLogger("annotate.log")
     setup_logger(logger, log_file_path)
     logger.info(f"The log file is created at {log_file_path}")
 
     # make Annotation objects per directory
     for i, annotation_dir in enumerate(glob(input_dirs)):
         if not path.exists(path.join(annotation_dir, "annotations.tsv")):
-            logger.warning("Skipping the fallowing directory because"
-                           f" no 'annotations.tsv' file was found: {annotation_dir}" )
+            logger.warning(
+                "Skipping the fallowing directory because"
+                f" no 'annotations.tsv' file was found: {annotation_dir}"
+            )
             continue
-        name = 'annotation_%s' % i
-        scaffolds = path.join(annotation_dir, 'scaffolds.fna')
-        genes_faa = path.join(annotation_dir, 'genes.faa')
-        genes_fna = path.join(annotation_dir, 'genes.fna')
-        gff = path.join(annotation_dir, 'genes.gff')
-        gbk = path.join(annotation_dir, 'genbank')
-        annotations = path.join(annotation_dir, 'annotations.tsv')
-        trnas = path.join(annotation_dir, 'trnas.tsv')
+        name = "annotation_%s" % i
+        scaffolds = path.join(annotation_dir, "scaffolds.fna")
+        genes_faa = path.join(annotation_dir, "genes.faa")
+        genes_fna = path.join(annotation_dir, "genes.fna")
+        gff = path.join(annotation_dir, "genes.gff")
+        gbk = path.join(annotation_dir, "genbank")
+        annotations = path.join(annotation_dir, "annotations.tsv")
+        trnas = path.join(annotation_dir, "trnas.tsv")
         if not path.isfile(trnas):
             logger.warning("No trnas.tsv file found in directory %s" % annotation_dir)
             trnas = None
-        rrnas = path.join(annotation_dir, 'rrnas.tsv')
+        rrnas = path.join(annotation_dir, "rrnas.tsv")
         if not path.isfile(rrnas):
             logger.warning("No rrnas.tsv file found in directory %s" % annotation_dir)
             rrnas = None
-        annotations_list.append(Annotation(name=name, scaffolds=scaffolds, genes_faa=genes_faa, genes_fna=genes_fna,
-                                           gff=gff, gbk=gbk, annotations=annotations, trnas=trnas, rrnas=rrnas))
+        annotations_list.append(
+            Annotation(
+                name=name,
+                scaffolds=scaffolds,
+                genes_faa=genes_faa,
+                genes_fna=genes_fna,
+                gff=gff,
+                gbk=gbk,
+                annotations=annotations,
+                trnas=trnas,
+                rrnas=rrnas,
+            )
+        )
     # run merge_annotations
     merge_annotations(annotations_list, output_dir, write_annotations=True)
-
diff --git a/mag_annotator/database_handler.py b/mag_annotator/database_handler.py
index 4bbafd59..0e9bf560 100644
--- a/mag_annotator/database_handler.py
+++ b/mag_annotator/database_handler.py
@@ -17,11 +17,25 @@
 from mag_annotator.database_setup import TABLE_NAME_TO_CLASS_DICT, create_description_db
 from mag_annotator.utils import divide_chunks, setup_logger
 
-SEARCH_DATABASES = {'kegg', 'kofam_hmm', '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', 'dbcan_fam_activities', 'vog_annotations')
+SEARCH_DATABASES = {
+    "kegg",
+    "kofam_hmm",
+    "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", "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
@@ -31,11 +45,11 @@
 
 
 def get_config_loc():
-    loc = getenv('DRAM_CONFIG_LOCATION')
+    loc = getenv("DRAM_CONFIG_LOCATION")
     if loc:
         return loc
     else:
-        return path.abspath(resource_filename('mag_annotator', 'CONFIG'))
+        return path.abspath(resource_filename("mag_annotator", "CONFIG"))
 
 
 def clear_dict(val):
@@ -46,7 +60,6 @@ def clear_dict(val):
 
 
 class DatabaseHandler:
-
     def __init__(self, logger, config_loc=None):
         # read in new configuration
         # TODO: validate config file after reading it in
@@ -63,47 +76,55 @@ def __init__(self, logger, config_loc=None):
 
         self.config_loc = config_loc
 
-
-
     def load_config(self, config_file):
         conf = json.loads(open(config_file).read())
         if len(conf) == 0:
-            self.logger.warn('There is no config information in the provided file')
+            self.logger.warn("There is no config information in the provided file")
             self.clear_config(write_config=False)
-        if 'dram_version' not in conf:
-            warnings.warn("The DRAM version in your config is empty."
-                          " This may not be a problem, but if this"
-                          " import fails then you should check that"
-                          " the origin of the file is valid.")
+        if "dram_version" not in conf:
+            warnings.warn(
+                "The DRAM version in your config is empty."
+                " This may not be a problem, but if this"
+                " import fails then you should check that"
+                " the origin of the file is valid."
+            )
             self.__construct_from_dram_pre_1_4_0(conf)
         else:
-            conf_version = conf.get('dram_version')
+            conf_version = conf.get("dram_version")
             if conf_version is None:
                 self.__construct_from_dram_pre_1_4_0(conf)
-            elif conf_version not in {current_dram_version, "1.4.0",  "1.4.0rc1", "1.4.0rc2", "1.4.0rc3", "1.4.0rc4"}: # Known suported versions
-                warnings.warn("The DRAM version in your config is not listed in the versions "
-                              "that are known to work. This may not be a problem, but if this "
-                              "import fails then you should contact suport.")
+            elif conf_version not in {
+                current_dram_version,
+                "1.4.0",
+                "1.4.0rc1",
+                "1.4.0rc2",
+                "1.4.0rc3",
+                "1.4.0rc4",
+            }:  # Known suported versions
+                warnings.warn(
+                    "The DRAM version in your config is not listed in the versions "
+                    "that are known to work. This may not be a problem, but if this "
+                    "import fails then you should contact suport."
+                )
             self.__construct_default(conf)
 
     def get_log_path(self):
-        path = self.config.get('log_path')
+        path = self.config.get("log_path")
         if path is None:
-            path = path.join(self.config_loc, 'database_processing.log')
+            path = path.join(self.config_loc, "database_processing.log")
         return path
 
-
-    def __construct_default(self, conf:dict):
+    def __construct_default(self, conf: dict):
         self.config = conf
 
         # set up description database connection
-        description_loc = self.config.get('description_db')
+        description_loc = self.config.get("description_db")
         if description_loc is None:
             self.session = None
-            warnings.warn('Database does not exist at path %s' % description_loc)
+            warnings.warn("Database does not exist at path %s" % description_loc)
         elif not path.exists(description_loc):
             self.session = None
-            warnings.warn('Database does not exist at path %s' % description_loc)
+            warnings.warn("Database does not exist at path %s" % description_loc)
         else:
             self.start_db_session()
 
@@ -119,59 +140,77 @@ def __construct_from_dram_pre_1_4_0(self, config_old):
         self.config = clear_dict(json.loads(open(system_config_loc).read()))
         self.config_loc = system_config_loc
 
-
         # read in configuration # TODO: validate config file after reading it in
-        self.config['search_databases'] = {
-            key: value for key, value in config_old.items() if key in SEARCH_DATABASES}
-        self.config['database_descriptions']  = {
-            key: value for key, value in config_old.items() if key in DATABASE_DESCRIPTIONS}
-        self.config['dram_sheets']  = {
-            key: value for key, value in config_old.items() if key in DRAM_SHEETS}
+        self.config["search_databases"] = {
+            key: value for key, value in config_old.items() if key in SEARCH_DATABASES
+        }
+        self.config["database_descriptions"] = {
+            key: value
+            for key, value in config_old.items()
+            if key in DATABASE_DESCRIPTIONS
+        }
+        self.config["dram_sheets"] = {
+            key: value for key, value in config_old.items() if key in DRAM_SHEETS
+        }
         self.config["dram_version"] = current_dram_version
 
         # set up description database connection
-        self.config['description_db'] = config_old.get('description_db')
-        if self.config.get('description_db') is None:
+        self.config["description_db"] = config_old.get("description_db")
+        if self.config.get("description_db") is None:
             self.session = None
-            warnings.warn('Database does not exist at path %s' % self.config.get('description_db'))
-        elif not path.exists(self.config.get('description_db')):
+            warnings.warn(
+                "Database does not exist at path %s" % self.config.get("description_db")
+            )
+        elif not path.exists(self.config.get("description_db")):
             self.session = None
-            warnings.warn('Database does not exist at path %s' % self.config.get('description_db'))
+            warnings.warn(
+                "Database does not exist at path %s" % self.config.get("description_db")
+            )
         else:
             self.start_db_session()
 
     def start_db_session(self):
-        engine = create_engine('sqlite:///%s' % self.config.get('description_db'))
+        engine = create_engine("sqlite:///%s" % self.config.get("description_db"))
         db_session = sessionmaker(bind=engine)
         self.session = db_session()
 
-
     # functions for adding descriptions to tables
     def add_descriptions_to_database(self, description_list, db_name, clear_table=True):
         description_class = TABLE_NAME_TO_CLASS_DICT[db_name]
         if clear_table:
             self.session.query(description_class).delete()
         # TODO: try batching
-        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]
+        )
         self.session.commit()
         self.session.expunge_all()
 
     # functions for getting descriptions from tables
     def get_description(self, annotation_id, db_name, return_ob=False):
-        return self.session.query(TABLE_NAME_TO_CLASS_DICT[db_name]).filter_by(id=annotation_id).one().description
-
-    def get_descriptions(self, ids, db_name, description_name='description'):
+        return (
+            self.session.query(TABLE_NAME_TO_CLASS_DICT[db_name])
+            .filter_by(id=annotation_id)
+            .one()
+            .description
+        )
+
+    def get_descriptions(self, ids, db_name, description_name="description"):
         description_class = TABLE_NAME_TO_CLASS_DICT[db_name]
         descriptions = [
             des
             for chunk in divide_chunks(list(ids), 499)
-            for des in self.session.query(description_class).filter(description_class.id.in_(chunk)).all()
+            for des in self.session.query(description_class)
+            .filter(description_class.id.in_(chunk))
+            .all()
         ]
         # [des for des in self.session.query(description_class).filter(description_class.id.in_(list(ids))).all() ]
         # [i.id for i in self.session.query(TABLE_NAME_TO_CLASS_DICT['dbcan_description']).all()]
         if len(descriptions) == 0:
-            warnings.warn("No descriptions were found for your id's. Does this %s look like an id from %s" % (list(ids)[0],
-                                                                                                     db_name))
+            warnings.warn(
+                "No descriptions were found for your id's. Does this %s look like an id from %s"
+                % (list(ids)[0], db_name)
+            )
         return {i.id: i.__dict__[description_name] for i in descriptions}
 
     @staticmethod
@@ -180,29 +219,49 @@ def get_database_names():
 
     def get_settings_str(self):
         out_str = ""
-        settings = self.config.get('setup_info')
+        settings = self.config.get("setup_info")
         if settings is None:
-            warnings.warn('there are no settings, the config is corrupted or too old.', DeprecationWarning)
-            return 'there are no settings, the config is corrupted or too old.'
+            warnings.warn(
+                "there are no settings, the config is corrupted or too old.",
+                DeprecationWarning,
+            )
+            return "there are no settings, the config is corrupted or too old."
         for i in ["search_databases", "database_descriptions", "dram_sheets"]:
             out_str += "\n"
             for k in self.config.get(i):
                 if settings.get(k) is not None:
                     out_str += f"\n{settings[k]['name']}:"
                     for l, w in settings[k].items():
-                        if l =='name':
+                        if l == "name":
                             continue
                         out_str += f"\n    {l.title()}: {w}"
         return out_str
 
-
-    def set_database_paths(self, kegg_loc=None, kofam_hmm_loc=None, kofam_ko_list_loc=None, uniref_loc=None,
-                           pfam_loc=None, pfam_hmm_loc=None, dbcan_loc=None, dbcan_fam_activities_loc=None,
-                           dbcan_subfam_ec_loc=None, viral_loc=None, peptidase_loc=None, vogdb_loc=None,
-                           vog_annotations_loc=None, description_db_loc=None, log_path_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 set_database_paths(
+        self,
+        kegg_loc=None,
+        kofam_hmm_loc=None,
+        kofam_ko_list_loc=None,
+        uniref_loc=None,
+        pfam_loc=None,
+        pfam_hmm_loc=None,
+        dbcan_loc=None,
+        dbcan_fam_activities_loc=None,
+        dbcan_subfam_ec_loc=None,
+        viral_refseq_loc=None,
+        peptidase_loc=None,
+        vogdb_loc=None,
+        vog_annotations_loc=None,
+        description_db_loc=None,
+        log_path_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,
+        update_description_db=False
+    ):
         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
@@ -210,143 +269,184 @@ def check_exists_and_add_to_location_dict(loc, old_value):
                 return path.realpath(loc)
             else:  # if the location doesn't exist then raise error
                 raise ValueError("Database location does not exist: %s" % loc)
+
         locs = {
-                "search_databases": {
-                  'kegg': kegg_loc,
-                  'kofam_hmm': kofam_hmm_loc,
-                  'kofam_ko_list': kofam_ko_list_loc,
-                  'uniref': uniref_loc,
-                  'pfam': pfam_loc,
-                  'dbcan': dbcan_loc,
-                  'viral': viral_loc,
-                  'peptidase': peptidase_loc,
-                  'vogdb': vogdb_loc,
-                },
-                "database_descriptions": {
-                  'pfam_hmm': pfam_hmm_loc,
-                  'dbcan_fam_activities': dbcan_fam_activities_loc,
-                  'dbcan_subfam_ec': dbcan_subfam_ec_loc,
-                  'vog_annotations': vog_annotations_loc,
-                },
-                "dram_sheets": {
-                  'genome_summary_form': genome_summary_form_loc,
-                  'module_step_form': module_step_form_loc,
-                  'etc_module_database': etc_module_database_loc,
-                  'function_heatmap_form': function_heatmap_form_loc,
-                  'amg_database': amg_database_loc,
-                },
+            "search_databases": {
+                "kegg": kegg_loc,
+                "kofam_hmm": kofam_hmm_loc,
+                "kofam_ko_list": kofam_ko_list_loc,
+                "uniref": uniref_loc,
+                "pfam": pfam_loc,
+                "dbcan": dbcan_loc,
+                "viral": viral_refseq_loc,
+                "peptidase": peptidase_loc,
+                "vogdb": vogdb_loc,
+            },
+            "database_descriptions": {
+                "pfam_hmm": pfam_hmm_loc,
+                "dbcan_fam_activities": dbcan_fam_activities_loc,
+                "dbcan_subfam_ec": dbcan_subfam_ec_loc,
+                "vog_annotations": vog_annotations_loc,
+            },
+            "dram_sheets": {
+                "genome_summary_form": genome_summary_form_loc,
+                "module_step_form": module_step_form_loc,
+                "etc_module_database": etc_module_database_loc,
+                "function_heatmap_form": function_heatmap_form_loc,
+                "amg_database": amg_database_loc,
+            },
         }
 
-        self.config.update({i:{
-            k:check_exists_and_add_to_location_dict(locs[i][k], self.config.get(i).get(k))
-            for k in locs[i]} for i in locs})
-
-        self.config['description_db'] = check_exists_and_add_to_location_dict(description_db_loc, self.config.get('description_db'))
-        self.config['log_path'] = check_exists_and_add_to_location_dict(log_path_loc, self.config.get('log_path_db'))
+        self.config.update(
+            {
+                i: {
+                    k: check_exists_and_add_to_location_dict(
+                        locs[i][k], self.config.get(i).get(k)
+                    )
+                    for k in locs[i]
+                }
+                for i in locs
+            }
+        )
+        if update_description_db:
+            self.populate_description_db(output_loc=description_db_loc)
+        else:
+            self.config["description_db"] = check_exists_and_add_to_location_dict(
+                description_db_loc, self.config.get("description_db")
+            )
+        self.config["log_path"] = check_exists_and_add_to_location_dict(
+            log_path_loc, self.config.get("log_path_db")
+        )
         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:
+        with open(config_loc, "w") as f:
             f.write(json.dumps(self.config, indent=2))
 
     @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]
+        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
+        ]
         mmseqs_headers_split = []
         mmseqs_ids_unique = set()
         mmseqs_ids_not_unique = set()
         # TODO this could be faster with numpy
         for i in mmseqs_headers:
-            header = {'id': i.split(' ')[0], 'description': i}
-            if header['id'] not in mmseqs_ids_unique:
+            header = {"id": i.split(" ")[0], "description": i}
+            if header["id"] not in mmseqs_ids_unique:
                 mmseqs_headers_split += [header]
-                mmseqs_ids_unique.add(header['id'])
+                mmseqs_ids_unique.add(header["id"])
             else:
-                mmseqs_ids_not_unique.add(header['id'])
+                mmseqs_ids_not_unique.add(header["id"])
         if len(mmseqs_ids_not_unique) > 0:
-            warnings.warn(f'There are {len(mmseqs_ids_not_unique)} non unique headers '
-                          f'in {mmseqs_db}! You should definitly investigate this!')
+            warnings.warn(
+                f"There are {len(mmseqs_ids_not_unique)} non unique headers "
+                f"in {mmseqs_db}! You should definitly investigate this!"
+            )
         return mmseqs_headers_split
 
     @staticmethod
     def process_pfam_descriptions(pfam_hmm):
-        if pfam_hmm.endswith('.gz'):
-            f = gzip.open(pfam_hmm, 'r').read().decode('utf-8')
+        if pfam_hmm.endswith(".gz"):
+            f = gzip.open(pfam_hmm, "r").read().decode("utf-8")
         else:
             f = open(pfam_hmm).read()
-        entries = f.strip().split('//')
+        entries = f.strip().split("//")
         description_list = list()
         for i, entry in enumerate(entries):
             if len(entry) > 0:
-                entry = entry.split('\n')
+                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})
+                    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, dbcan_subfam_ec):
         def line_reader(line):
-            if not line.startswith('#') and len(line.strip()) != 0:
+            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:])
+                    description = " ".join(line[1:])
                 else:
-                    description = ' '.join(line)
-                return pd.DataFrame({'id': line[0], 'description': description.replace('\n', ' ')}, index=[0])
+                    description = " ".join(line)
+                return pd.DataFrame(
+                    {"id": line[0], "description": description.replace("\n", " ")},
+                    index=[0],
+                )
+
         with open(dbcan_fam_activities) as f:
             description_data = pd.concat([line_reader(line) for line in f.readlines()])
 
-        ec_data = (pd.read_csv(dbcan_subfam_ec, sep='\t',names=['id', 'id2','ec'], comment='#')[['id', 'ec']]
-                           .drop_duplicates())
-        ec_data = (pd.concat([ec_data['id'],
-                                     ec_data['ec'].str.split('|', expand=True)]
-                                    ,axis=1)
-                           .melt(id_vars='id',value_name='ec')
-                           .dropna(subset=['ec'])[['id', 'ec']]
-                           .groupby('id')
-                           .apply(lambda x: ','.join(x['ec'].unique()))
-                           )
-        ec_data = pd.DataFrame(ec_data, columns=['ec']).reset_index()
-        data = pd.merge(description_data, ec_data, how='outer', on='id').fillna('')
+        ec_data = pd.read_csv(
+            dbcan_subfam_ec, sep="\t", names=["id", "id2", "ec"], comment="#"
+        )[["id", "ec"]].drop_duplicates()
+        ec_data = (
+            pd.concat(
+                [ec_data["id"], ec_data["ec"].str.split("|", expand=True)], axis=1
+            )
+            .melt(id_vars="id", value_name="ec")
+            .dropna(subset=["ec"])[["id", "ec"]]
+            .groupby("id")
+            .apply(lambda x: ",".join(x["ec"].unique()))
+        )
+        ec_data = pd.DataFrame(ec_data, columns=["ec"]).reset_index()
+        data = pd.merge(description_data, ec_data, how="outer", on="id").fillna("")
         return [i.to_dict() for _, i in data.iterrows()]
 
-
     @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()]
+        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, select_db=None, update_config=True,  erase_old_db=False):
-        if self.config.get('description_db') is None and output_loc is None:  # description db location must be set somewhere
-            self.logger.critical('Must provide output location if description db location is not set in configuration')
-            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.config['description_db']= output_loc
+    def populate_description_db(
+        self, output_loc=None, select_db=None, update_config=True, erase_old_db=False
+    ):
+        if (
+            self.config.get("description_db") is None and output_loc is None
+        ):  # description db location must be set somewhere
+            self.logger.critical(
+                "Must provide output location if description db location is not set in configuration"
+            )
+            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.config["description_db"] = output_loc
             self.start_db_session()
         # I don't think this is needed
-        if path.exists(self.config.get('description_db')) and erase_old_db:
-            remove(self.config.get('description_db'))
-        create_description_db(self.config.get('description_db'))
+        if path.exists(self.config.get("description_db")) and erase_old_db:
+            remove(self.config.get("description_db"))
+        create_description_db(self.config.get("description_db"))
+
         def check_db(db_name, db_function):
             # TODO add these sorts of checks to a separate function
             # if self.config.get('search_databases').get(db_name) is None:
@@ -357,30 +457,44 @@ def check_db(db_name, db_function):
             #                 f"{self.config['search_databases'][db_name]}")
             #     return
             self.add_descriptions_to_database(
-                db_function(),
-                f'{db_name}_description',
-                clear_table=True)
-            self.config['setup_info'][db_name]['description_db_updated'] = \
-                datetime.now().strftime("%m/%d/%Y, %H:%M:%S")
-            self.logger.info(f'Description updated for the {db_name} database')
+                db_function(), f"{db_name}_description", clear_table=True
+            )
+            self.config["setup_info"][db_name][
+                "description_db_updated"
+            ] = datetime.now().strftime("%m/%d/%Y, %H:%M:%S")
+            self.logger.info(f"Description updated for the {db_name} database")
+
         # fill database
-        mmseqs_database = ['kegg', 'uniref',  'viral', 'peptidase']
-        process_functions = {i:partial(self.make_header_dict_from_mmseqs_db,
-                                       self.config['search_databases'][i])
-                             for i in mmseqs_database
-                             if self.config['search_databases'][i] is not None}
+        mmseqs_database = ["kegg", "uniref", "viral", "peptidase"]
+        process_functions = {
+            i: partial(
+                self.make_header_dict_from_mmseqs_db, self.config["search_databases"][i]
+            )
+            for i in mmseqs_database
+            if self.config["search_databases"][i] is not None
+        }
         # Use table names
-        process_functions.update({
-            'pfam': partial(self.process_pfam_descriptions,
-                                self.config.get('database_descriptions')['pfam_hmm']),
-            'dbcan': partial(self.process_dbcan_descriptions,
-                                            self.config.get('database_descriptions')['dbcan_fam_activities'],
-                                            self.config.get('database_descriptions')['dbcan_subfam_ec']),
-            'vogdb': partial(self.process_vogdb_descriptions,
-                                         self.config.get('database_descriptions')['vog_annotations'])
-        })
+        process_functions.update(
+            {
+                "pfam": partial(
+                    self.process_pfam_descriptions,
+                    self.config.get("database_descriptions")["pfam_hmm"],
+                ),
+                "dbcan": partial(
+                    self.process_dbcan_descriptions,
+                    self.config.get("database_descriptions")["dbcan_fam_activities"],
+                    self.config.get("database_descriptions")["dbcan_subfam_ec"],
+                ),
+                "vogdb": partial(
+                    self.process_vogdb_descriptions,
+                    self.config.get("database_descriptions")["vog_annotations"],
+                ),
+            }
+        )
         if select_db is not None:
-            process_functions = {i:k for i, k in process_functions.items() if i in select_db}
+            process_functions = {
+                i: k for i, k in process_functions.items() if i in select_db
+            }
 
         for i, k in process_functions.items():
             check_db(i, k)
@@ -388,83 +502,116 @@ def check_db(db_name, db_function):
         if update_config:  # if new description db is set then save it
             self.write_config()
 
-    def filter_db_locs(self, low_mem_mode=False, use_uniref=True, use_vogdb=True, master_list=None):
+    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.config['search_databases'].keys()
+            dbs_to_use = self.config["search_databases"].keys()
         else:
             dbs_to_use = master_list
         # filter out dbs for low mem mode
         if low_mem_mode:
-            if ('kofam_hmm' not in self.config.get('search_databases')) or ('kofam_ko_list' not in self.config.get('search_databases')):
-                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')]
+            if ("kofam_hmm" not in self.config.get("search_databases")) or (
+                "kofam_ko_list" not in self.config.get("search_databases")
+            ):
+                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.config.get('search_databases'):
-                warnings.warn('Sequences will not be annoated against uniref as it is not configured for use in DRAM')
+            if "uniref" not in self.config.get("search_databases"):
+                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']
+            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.config.get('search_databases'):
-                warnings.warn('Sequences will not be annoated against VOGDB as it is not configured for use in DRAM')
+            if "vogdb" not in self.config.get("search_databases"):
+                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.config['search_databases'] = {key: value for key, value in self.config.get('search_databases').items() if key in dbs_to_use}
+            dbs_to_use = [i for i in dbs_to_use if i != "vogdb"]
+        self.config["search_databases"] = {
+            key: value
+            for key, value in self.config.get("search_databases").items()
+            if key in dbs_to_use
+        }
 
     def clear_config(self, write_config=False):
         self.config = {
-                        "search_databases": {},
-                        "database_descriptions": {},
-                        "dram_sheets": {},
-                        "dram_version": current_dram_version,
-                        "description_db": None,
-                        "setup_info": {},
-                        "log_path": None
-                      }
+            "search_databases": {},
+            "database_descriptions": {},
+            "dram_sheets": {},
+            "dram_version": current_dram_version,
+            "description_db": None,
+            "setup_info": {},
+            "log_path": None,
+        }
         if write_config:
             self.write_config()
 
 
-def set_database_paths(clear_config=False, update_description_db=False, **kargs):
-    #TODO Add tests
+def set_database_paths(clear_config=False, **kargs):
+    # TODO Add tests
     db_handler = DatabaseHandler(None)
     if clear_config:
         db_handler.clear_config(write_config=True)
     db_handler.set_database_paths(**kargs, write_config=True)
-    if update_description_db:
-        db_handler.populate_description_db()
 
 
 def print_database_locations(config_loc=None):
     conf = DatabaseHandler(None, config_loc)
     # search databases
-    print('Processed search databases')
-    print('KEGG db: %s' % conf.config.get('search_databases').get('kegg'))
-    print('KOfam db: %s' % conf.config.get('search_databases').get('kofam_hmm'))
-    print('KOfam KO list: %s' % conf.config.get('search_databases').get('kofam_ko_list'))
-    print('UniRef db: %s' % conf.config.get('search_databases').get('uniref'))
-    print('Pfam db: %s' % conf.config.get('search_databases').get('pfam'))
-    print('dbCAN db: %s' % conf.config.get('search_databases').get('dbcan'))
-    print('RefSeq Viral db: %s' % conf.config.get('search_databases').get('viral'))
-    print('MEROPS peptidase db: %s' % conf.config.get('search_databases').get('peptidase'))
-    print('VOGDB db: %s' % conf.config.get('search_databases').get('vogdb'))
+    print("Processed search databases")
+    print("KEGG db: %s" % conf.config.get("search_databases").get("kegg"))
+    print("KOfam db: %s" % conf.config.get("search_databases").get("kofam_hmm"))
+    print(
+        "KOfam KO list: %s" % conf.config.get("search_databases").get("kofam_ko_list")
+    )
+    print("UniRef db: %s" % conf.config.get("search_databases").get("uniref"))
+    print("Pfam db: %s" % conf.config.get("search_databases").get("pfam"))
+    print("dbCAN db: %s" % conf.config.get("search_databases").get("dbcan"))
+    print("RefSeq Viral db: %s" % conf.config.get("search_databases").get("viral"))
+    print(
+        "MEROPS peptidase db: %s" % conf.config.get("search_databases").get("peptidase")
+    )
+    print("VOGDB db: %s" % conf.config.get("search_databases").get("vogdb"))
     # database descriptions used during description db population
-    print('Descriptions of search database entries')
-    print('Pfam hmm dat: %s' % conf.config.get('database_descriptions').get('pfam_hmm'))
-    print('dbCAN family activities: %s' % conf.config.get('database_descriptions').get('dbcan_fam_activities'))
-    print('VOG annotations: %s' % conf.config.get('database_descriptions').get('vog_annotations'))
+    print("Descriptions of search database entries")
+    print("Pfam hmm dat: %s" % conf.config.get("database_descriptions").get("pfam_hmm"))
+    print(
+        "dbCAN family activities: %s"
+        % conf.config.get("database_descriptions").get("dbcan_fam_activities")
+    )
+    print(
+        "VOG annotations: %s"
+        % conf.config.get("database_descriptions").get("vog_annotations")
+    )
     print()
     # description database
-    print('Description db: %s' % conf.config.get('description_db'))
+    print("Description db: %s" % conf.config.get("description_db"))
     print()
     # DRAM sheets
-    print('DRAM distillation sheets')
-    print('Genome summary form: %s' % conf.config.get('dram_sheets').get('genome_summary_form'))
-    print('Module step form: %s' % conf.config.get('dram_sheets').get('module_step_form'))
-    print('ETC module database: %s' % conf.config.get('dram_sheets').get('etc_module_database'))
-    print('Function heatmap form: %s' % conf.config.get('dram_sheets').get('function_heatmap_form'))
-    print('AMG database: %s' % conf.config.get('dram_sheets').get('amg_database'))
+    print("DRAM distillation sheets")
+    print(
+        "Genome summary form: %s"
+        % conf.config.get("dram_sheets").get("genome_summary_form")
+    )
+    print(
+        "Module step form: %s" % conf.config.get("dram_sheets").get("module_step_form")
+    )
+    print(
+        "ETC module database: %s"
+        % conf.config.get("dram_sheets").get("etc_module_database")
+    )
+    print(
+        "Function heatmap form: %s"
+        % conf.config.get("dram_sheets").get("function_heatmap_form")
+    )
+    print("AMG database: %s" % conf.config.get("dram_sheets").get("amg_database"))
 
 
 def print_database_settings(config_loc=None):
@@ -472,8 +619,7 @@ def print_database_settings(config_loc=None):
     print(conf.get_settings_str())
 
 
-
-def populate_description_db(output_loc=None, select_db=None,  config_loc=None):
+def populate_description_db(output_loc=None, select_db=None, config_loc=None):
     db_handler = DatabaseHandler(None, config_loc)
     db_handler.populate_description_db(output_loc, select_db)
 
@@ -491,31 +637,33 @@ def import_config(config_loc):
     db_handler = DatabaseHandler(None, config_loc)
     with open(system_config, "w") as outfile:
         json.dump(db_handler.config, outfile, indent=2)
-    print('Import, appears to be successfull.')
+    print("Import, appears to be successfull.")
 
 
-def mv_db_folder(new_location:str='./', old_config_file:str=None):
+def mv_db_folder(new_location: str = "./", old_config_file: str = None):
     new_location = path.abspath(new_location)
     old_config_file = path.abspath(old_config_file)
     db_handler = DatabaseHandler(None)
     if old_config_file is not None:
         db_handler.load_config(old_config_file)
     paths = ["search_databases", "dram_sheets", "database_descriptions"]
-    def auto_move_path(k:str, v:str):
+
+    def auto_move_path(k: str, v: str):
         if v is None:
             db_handler.logger.warn(f"The path for {k} was not set, so can't update.")
             return
         new_path = path.join(new_location, path.basename(v))
         if not path.exists(new_path):
-            db_handler.logger.warn(f"There is no file at path {new_path},"
-                                   f" so no new location will be set for {k}.")
+            db_handler.logger.warn(
+                f"There is no file at path {new_path},"
+                f" so no new location will be set for {k}."
+            )
             return
         db_handler.logger.info(f"Moving {k} to {new_path}")
         db_handler.set_database_paths(**{f"{k}_loc": new_path}, write_config=True)
-    auto_move_path('log_path', db_handler.config.get('log_path'))
-    auto_move_path('description_db', db_handler.config.get('description_db'))
+
+    auto_move_path("log_path", db_handler.config.get("log_path"))
+    auto_move_path("description_db", db_handler.config.get("description_db"))
     for i in paths:
         for k, v in db_handler.config.get(i).items():
             auto_move_path(k, v)
-
-
diff --git a/scripts/DRAM-setup.py b/scripts/DRAM-setup.py
index 0ecc17a1..430d2bc9 100644
--- a/scripts/DRAM-setup.py
+++ b/scripts/DRAM-setup.py
@@ -113,11 +113,11 @@ def print_version():
     set_db_locs_parser.add_argument('--pfam_hmm_loc', default=None, help='pfam hmm .dat file to get PF descriptions')
     set_db_locs_parser.add_argument('--dbcan_loc', default=None,
                                     help='hmm file for dbcan, already processed with hmmpress')
-    set_db_locs_parser.add_argument('--dbcan_fam_activities', default=None, help='CAZY family activities file')
-    set_db_locs_parser.add_argument('--dbcan_subfam_ec', default=None, help='CAZY sub-family ECs file')
+    set_db_locs_parser.add_argument('--dbcan_fam_activities_loc', default=None, help='CAZY family activities file')
+    set_db_locs_parser.add_argument('--dbcan_subfam_ec_loc', default=None, help='CAZY sub-family ECs file')
     set_db_locs_parser.add_argument('--vogdb_loc', default=None,
                                     help='hmm file for vogdb, already processed with hmmpress')
-    set_db_locs_parser.add_argument('--vog_annotations', default=None,
+    set_db_locs_parser.add_argument('--vog_annotations_loc', default=None,
                                     help='vog annotations file') # add loc to vog_annotations to match the rest
 
     set_db_locs_parser.add_argument('--viral_refseq_loc', default=None,


Test patch:
diff --git a/tests/test_annotate_bins.py b/tests/test_annotate_bins.py
index c2681d83..2f335dad 100644
--- a/tests/test_annotate_bins.py
+++ b/tests/test_annotate_bins.py
@@ -1,4 +1,3 @@
-from re import L
 import pytest
 
 import os
@@ -14,42 +13,68 @@
 from skbio.io import read as read_sequence
 
 from mag_annotator.utils import make_mmseqs_db, parse_hmmsearch_domtblout, setup_logger
-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_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, \
-    annotate_gff, make_gbk_from_gff_and_fasta, make_trnas_interval, make_rrnas_interval,\
-    add_intervals_to_gff, vogdb_hmmscan_formater, dbcan_hmmscan_formater, \
-    kofam_hmmscan_formater
-
+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_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,
+    annotate_gff,
+    make_gbk_from_gff_and_fasta,
+    make_trnas_interval,
+    make_rrnas_interval,
+    add_intervals_to_gff,
+    vogdb_hmmscan_formater,
+    dbcan_hmmscan_formater,
+    kofam_hmmscan_formater,
+)
 
 
 @pytest.fixture()
 def fasta_loc():
-    return os.path.join('tests', 'data', 'NC_001422.fasta')
+    return os.path.join("tests", "data", "NC_001422.fasta")
 
 
 @pytest.fixture()
 def logger(tmpdir):
-    logger = logging.getLogger('test_log')
+    logger = logging.getLogger("test_log")
     setup_logger(logger)
     return logger
 
+
 def test_filter_fasta(fasta_loc, tmpdir):
     filtered_seq_5000 = filter_fasta(fasta_loc, min_len=5000)
     assert len(filtered_seq_5000) == 1
     assert len(filtered_seq_5000[0]) == 5386
     filtered_seq_6000 = filter_fasta(fasta_loc, min_len=6000)
     assert len(filtered_seq_6000) == 0
-    filt_fasta = tmpdir.mkdir('test_filt').join('filtered_fasta.fasta')
+    filt_fasta = tmpdir.mkdir("test_filt").join("filtered_fasta.fasta")
     filter_fasta(fasta_loc, min_len=5000, output_loc=str(filt_fasta))
     assert os.path.isfile(filt_fasta)
 
 
 @pytest.fixture()
 def prodigal_dir(fasta_loc, tmpdir, logger):
-    prodigal_output = tmpdir.mkdir('prodigal_output')
+    prodigal_output = tmpdir.mkdir("prodigal_output")
     gff, fna, faa = run_prodigal(fasta_loc, str(prodigal_output), logger)
     return gff, fna, faa
 
@@ -77,32 +102,34 @@ def test_run_prodigal(prodigal_gff, prodigal_fna, prodigal_faa):
 
 @pytest.fixture()
 def mmseqs_db_dir(tmpdir):
-    output_loc = tmpdir.mkdir('make_mmseqs_db_test')
+    output_loc = tmpdir.mkdir("make_mmseqs_db_test")
     return output_loc
 
 
 @pytest.fixture()
 def mmseqs_db(prodigal_faa, mmseqs_db_dir, logger):
-    output_file = str(mmseqs_db_dir.join('mmseqs_db.mmsdb'))
+    output_file = str(mmseqs_db_dir.join("mmseqs_db.mmsdb"))
     make_mmseqs_db(prodigal_faa, output_file, logger, True, 1)
     return output_file
 
 
 @pytest.fixture()
 def phix_proteins():
-    return os.path.join('tests', 'data', 'NC_001422.faa')
+    return os.path.join("tests", "data", "NC_001422.faa")
 
 
 @pytest.fixture()
 def target_mmseqs_db(mmseqs_db_dir, phix_proteins, logger):
-    output_file = str(mmseqs_db_dir.join('target.mmsdb'))
+    output_file = str(mmseqs_db_dir.join("target.mmsdb"))
     make_mmseqs_db(phix_proteins, output_file, logger, True, 1)
     return output_file
 
 
 @pytest.fixture()
 def best_hits_loc(mmseqs_db, target_mmseqs_db, mmseqs_db_dir, logger):
-    best_hits_loc = get_best_hits(mmseqs_db, target_mmseqs_db, logger, mmseqs_db_dir,threads=1, verbose=False)
+    best_hits_loc = get_best_hits(
+        mmseqs_db, target_mmseqs_db, logger, mmseqs_db_dir, threads=1, verbose=False
+    )
     return best_hits_loc
 
 
@@ -111,10 +138,12 @@ def test_get_best_hits(best_hits_loc):
 
 
 @pytest.fixture()
-def reverse_best_hits_loc(best_hits_loc, mmseqs_db, target_mmseqs_db, mmseqs_db_dir, logger):
+def reverse_best_hits_loc(
+    best_hits_loc, mmseqs_db, target_mmseqs_db, mmseqs_db_dir, logger
+):
     reverse_best_hits_loc = get_reciprocal_best_hits(
-        mmseqs_db, target_mmseqs_db, logger, mmseqs_db_dir, threads=1,
-        verbose=False)
+        mmseqs_db, target_mmseqs_db, logger, mmseqs_db_dir, threads=1, verbose=False
+    )
     return reverse_best_hits_loc
 
 
@@ -124,144 +153,180 @@ def test_get_reciprocal_best_hits(reverse_best_hits_loc):
 
 @pytest.fixture()
 def processed_hits():
-    forward = os.path.join('tests', 'data', 'query_target_hits.b6')
-    reverse = os.path.join('tests', 'data', 'target_query_hits.b6')
+    forward = os.path.join("tests", "data", "query_target_hits.b6")
+    reverse = os.path.join("tests", "data", "target_query_hits.b6")
     processed_hits = process_reciprocal_best_hits(forward, reverse)
     return processed_hits
 
+
 def test_dbcan_hmmscan_formater():
     # TODO can we test with db-handler?
-    output_expt = pd.DataFrame({"bin_1.scaffold_1": ["GT4; GT5", "GT5.hmm"],
-                           "bin_1.scaffold_2": ["GH1", "GH1.hmm"]}, index=["cazy_ids", "cazy_best_hit"]).T
-    input_b6 = os.path.join('tests', 'data', 'unformatted_cazy.b6')
+    output_expt = pd.DataFrame(
+        {
+            "bin_1.scaffold_1": ["GT4; GT5", "GT5.hmm"],
+            "bin_1.scaffold_2": ["GH1", "GH1.hmm"],
+        },
+        index=["cazy_ids", "cazy_best_hit"],
+    ).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 = 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"
 
 
-
 def test_process_reciprocal_best_hits(processed_hits):
     assert processed_hits.shape == (7, 5)
-    assert set(processed_hits.loc[processed_hits.target_RBH].index) == {'NC_001422.1_5', 'NC_001422.1_4',
-                                                                        'NC_001422.1_7', 'NC_001422.1_6'}
+    assert set(processed_hits.loc[processed_hits.target_RBH].index) == {
+        "NC_001422.1_5",
+        "NC_001422.1_4",
+        "NC_001422.1_7",
+        "NC_001422.1_6",
+    }
 
 
 def test_get_kegg_description():
-    header_dict = {'aad:TC41_2367': 'aad:TC41_2367  ABC-type molybdate transport system periplasmic component-like '
-                                    'protein',
-                   'aar:Acear_0854': 'aar:Acear_0854  hypothetical protein; K05810 conserved hypothetical protein',
-                   'aar:Acear_1520': 'aar:Acear_1520  hypothetical protein'}
-    kegg_hits_data = [['aad:TC41_2367', 10e-5],
-                      ['aar:Acear_0854', 10e-6],
-                      ['aar:Acear_1520', 10e-10]]
-    kegg_hits = pd.DataFrame(kegg_hits_data, index=['gene1', 'gene2', 'gene3'], columns=['kegg_hit', 'eVal'])
+    header_dict = {
+        "aad:TC41_2367": "aad:TC41_2367  ABC-type molybdate transport system periplasmic component-like "
+        "protein",
+        "aar:Acear_0854": "aar:Acear_0854  hypothetical protein; K05810 conserved hypothetical protein",
+        "aar:Acear_1520": "aar:Acear_1520  hypothetical protein",
+    }
+    kegg_hits_data = [
+        ["aad:TC41_2367", 10e-5],
+        ["aar:Acear_0854", 10e-6],
+        ["aar:Acear_1520", 10e-10],
+    ]
+    kegg_hits = pd.DataFrame(
+        kegg_hits_data, index=["gene1", "gene2", "gene3"], columns=["kegg_hit", "eVal"]
+    )
     kegg_hits_add_description = get_kegg_description(kegg_hits, header_dict)
     print(kegg_hits_add_description.head())
     assert kegg_hits_add_description.shape == (3, 4)
-    assert kegg_hits_add_description.loc['gene2', 'ko_id'] == 'K05810'
-    assert kegg_hits_add_description.loc['gene2', 'kegg_genes_id'] == 'aar:Acear_0854'
-    assert kegg_hits_add_description.loc['gene1', 'ko_id'] == ''
-    assert kegg_hits_add_description.loc['gene1', 'kegg_genes_id'] == 'aad:TC41_2367'
+    assert kegg_hits_add_description.loc["gene2", "ko_id"] == "K05810"
+    assert kegg_hits_add_description.loc["gene2", "kegg_genes_id"] == "aar:Acear_0854"
+    assert kegg_hits_add_description.loc["gene1", "ko_id"] == ""
+    assert kegg_hits_add_description.loc["gene1", "kegg_genes_id"] == "aad:TC41_2367"
 
 
 def test_get_uniref_description():
-    header_dict = {'UniRef90_A0A139CGD2': 'UniRef90_A0A139CGD2 Phosphate transport system permease protein PstA n=2 '
-                                          'Tax=Candidatus Frackibacter TaxID=2017975 RepID=A0A139CGD2_9FIRM',
-                   'UniRef90_2642661139': 'UniRef90_2642661139 Ga0073286_10147 conserved hypothetical protein n=1 '
-                                          'Tax=Fuchsiella alkaliacetigena WG11 RepID=Ga0073286_10147'}
-    uniref_hits_data = [['UniRef90_A0A139CGD2', 10e-5],
-                        ['UniRef90_2642661139', 10e-20]]
-    uniref_hits = pd.DataFrame(uniref_hits_data, index=['gene1', 'gene2'], columns=['uniref_hit', 'eVal'])
+    header_dict = {
+        "UniRef90_A0A139CGD2": "UniRef90_A0A139CGD2 Phosphate transport system permease protein PstA n=2 "
+        "Tax=Candidatus Frackibacter TaxID=2017975 RepID=A0A139CGD2_9FIRM",
+        "UniRef90_2642661139": "UniRef90_2642661139 Ga0073286_10147 conserved hypothetical protein n=1 "
+        "Tax=Fuchsiella alkaliacetigena WG11 RepID=Ga0073286_10147",
+    }
+    uniref_hits_data = [["UniRef90_A0A139CGD2", 10e-5], ["UniRef90_2642661139", 10e-20]]
+    uniref_hits = pd.DataFrame(
+        uniref_hits_data, index=["gene1", "gene2"], columns=["uniref_hit", "eVal"]
+    )
     uniref_hits_add_description = get_uniref_description(uniref_hits, header_dict)
     assert uniref_hits_add_description.shape == (2, 4)
-    assert uniref_hits_add_description.loc['gene1', 'uniref_id'] == 'A0A139CGD2_9FIRM'
-    assert uniref_hits_add_description.loc['gene1', 'uniref_taxonomy'] == 'Candidatus Frackibacter'
+    assert uniref_hits_add_description.loc["gene1", "uniref_id"] == "A0A139CGD2_9FIRM"
+    assert (
+        uniref_hits_add_description.loc["gene1", "uniref_taxonomy"]
+        == "Candidatus Frackibacter"
+    )
 
 
 def test_get_basic_description():
-    header_dict = {'YP_009015653.1': 'YP_009015653.1 gp350 [Bacillus virus G]',
-                   'NP_077550.1': 'NP_077550.1 EsV-1-65 [Ectocarpus siliculosus virus 1]'}
-    viral_hits_data = [['NP_077550.1'],
-                       ['YP_009015653.1']]
-    viral_hits = pd.DataFrame(viral_hits_data, index=['gene1', 'gene2'], columns=['viral_hit'])
+    header_dict = {
+        "YP_009015653.1": "YP_009015653.1 gp350 [Bacillus virus G]",
+        "NP_077550.1": "NP_077550.1 EsV-1-65 [Ectocarpus siliculosus virus 1]",
+    }
+    viral_hits_data = [["NP_077550.1"], ["YP_009015653.1"]]
+    viral_hits = pd.DataFrame(
+        viral_hits_data, index=["gene1", "gene2"], columns=["viral_hit"]
+    )
     viral_hits_add_description = get_basic_description(viral_hits, header_dict)
     assert viral_hits_add_description.shape == (2, 2)
-    assert viral_hits_add_description.loc['gene1', 'viral_id'] == 'NP_077550.1'
+    assert viral_hits_add_description.loc["gene1", "viral_id"] == "NP_077550.1"
 
 
 def test_get_peptidase_description():
-    header_dict = {'MER0025711': 'MER0025711 - family S12 unassigned peptidases (Cytophaga hutchinsonii) [S12.UPW]#S12#'
-                                 '{peptidase unit: 588-980}~source ZP_00119535~',
-                   'MER0068848': 'MER0068848 - family C56 non-peptidase homologues (Silicibacter sp. TM1040) [C56.UNW]#'
-                                 'C56#{peptidase unit: 31-167}~source YP_611907~'}
-    peptidase_hits_data = [['MER0025711'], ['MER0068848']]
-    peptidase_hits = pd.DataFrame(peptidase_hits_data, index=['gene1', 'gene2'], columns=['peptidase_hit'])
-    peptidase_hits_add_description = get_peptidase_description(peptidase_hits, header_dict)
+    header_dict = {
+        "MER0025711": "MER0025711 - family S12 unassigned peptidases (Cytophaga hutchinsonii) [S12.UPW]#S12#"
+        "{peptidase unit: 588-980}~source ZP_00119535~",
+        "MER0068848": "MER0068848 - family C56 non-peptidase homologues (Silicibacter sp. TM1040) [C56.UNW]#"
+        "C56#{peptidase unit: 31-167}~source YP_611907~",
+    }
+    peptidase_hits_data = [["MER0025711"], ["MER0068848"]]
+    peptidase_hits = pd.DataFrame(
+        peptidase_hits_data, index=["gene1", "gene2"], columns=["peptidase_hit"]
+    )
+    peptidase_hits_add_description = get_peptidase_description(
+        peptidase_hits, header_dict
+    )
     assert peptidase_hits_add_description.shape == (2, 3)
-    assert peptidase_hits_add_description.loc['gene1', 'peptidase_id'] == 'MER0025711'
-    assert peptidase_hits_add_description.loc['gene1', 'peptidase_family'] == 'S12'
-    assert peptidase_hits_add_description.loc['gene2', 'peptidase_hit'] == 'MER0068848 - family C56 non-peptidase ' \
-                                                                           'homologues (Silicibacter sp. TM1040) [C56' \
-                                                                           '.UNW]#C56#{peptidase unit: 31-167}~source' \
-                                                                           ' YP_611907~'
+    assert peptidase_hits_add_description.loc["gene1", "peptidase_id"] == "MER0025711"
+    assert peptidase_hits_add_description.loc["gene1", "peptidase_family"] == "S12"
+    assert (
+        peptidase_hits_add_description.loc["gene2", "peptidase_hit"]
+        == "MER0068848 - family C56 non-peptidase "
+        "homologues (Silicibacter sp. TM1040) [C56"
+        ".UNW]#C56#{peptidase unit: 31-167}~source"
+        " YP_611907~"
+    )
 
 
 def test_get_sig_row():
-    names = ['target_start', 'target_end', 'target_length', 'full_evalue']
+    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():
-    phix_seq = ">NC_001422.1_1 # 51 # 221 # 1 # ID=1_1;partial=00;start_type=ATG;rbs_motif=GGA/GAG/AGG;rbs_spacer=5-" \
-               "10bp;gc_cont=0.404\n" \
-               "MSRKIILIKQELLLLVYELNRSGLLAENEKIRPILAQLEKLLLCDLSPSTNDSVKN*\n" \
-               ">NC_001422.1_2 # 390 # 848 # 1 # ID=1_2;partial=00;start_type=ATG;rbs_motif=GGA/GAG/AGG;rbs_spacer=" \
-               "5-10bp;gc_cont=0.468\n" \
-               "MSQVTEQSVRFQTALASIKLIQASAVLDLTEDDFDFLTSNKVWIATDRSRARRCVEACVY\n" \
-               "GTLDFVGYPRFPAPVEFIAAVIAYYVHPVNIQTACLIMEGAEFTENIINGVERPVKAAEL\n" \
-               "FAFTLRVRAGNTDVLTDAEENVRQKLRAEGVM*\n" \
-               ">NC_001422.1_3 # 848 # 964 # 1 # ID=1_3;partial=00;start_type=ATG;rbs_motif=AGGAG;rbs_spacer=5-10bp;" \
-               "gc_cont=0.496\n" \
-               "MSKGKKRSGARPGRPQPLRGTKGKRKGARLWYVGGQQF*\n" \
-               ">NC_001422.1_4 # 1001 # 2284 # 1 # ID=1_4;partial=00;start_type=ATG;rbs_motif=GGAG/GAGG;rbs_spacer=5" \
-               "-10bp;gc_cont=0.449\n" \
-               "MSNIQTGAERMPHDLSHLGFLAGQIGRLITISTTPVIAGDSFEMDAVGALRLSPLRRGLA\n" \
-               "IDSTVDIFTFYVPHRHVYGEQWIKFMKDGVNATPLPTVNTTGYIDHAAFLGTINPDTNKI\n" \
-               "PKHLFQGYLNIYNNYFKAPWMPDRTEANPNELNQDDARYGFRCCHLKNIWTAPLPPETEL\n" \
-               "SRQMTTSTTSIDIMGLQAAYANLHTDQERDYFMQRYHDVISSFGGKTSYDADNRPLLVMR\n" \
-               "SNLWASGYDVDGTDQTSLGQFSGRVQQTYKHSVPRFFVPEHGTMFTLALVRFPPTATKEI\n" \
-               "QYLNAKGALTYTDIAGDPVLYGNLPPREISMKDVFRSGDSSKKFKIAEGQWYRYAPSYVS\n" \
-               "PAYHLLEGFPFIQEPPSGDLQERVLIRHHDYDQCFQSVQLLQWNSQVKFNVTVYRNLPTT\n" \
-               "RDSIMTS*\n" \
-               ">NC_001422.1_5 # 2395 # 2922 # 1 # ID=1_5;partial=00;start_type=ATG;rbs_motif=AGGAG;rbs_spacer=5-10" \
-               "bp;gc_cont=0.420\n" \
-               "MFQTFISRHNSNFFSDKLVLTSVTPASSAPVLQTPKATSSTLYFDSLTVNAGNGGFLHCI\n" \
-               "QMDTSVNAANQVVSVGADIAFDADPKFFACLVRFESSSVPTTLPTAYDVYPLNGRHDGGY\n" \
-               "YTVKDCVTIDVLPRTPGNNVYVGFMVWSNFTATKCRGLVSLNQVIKEIICLQPLK*\n" \
-               ">NC_001422.1_6 # 2931 # 3917 # 1 # ID=1_6;partial=00;start_type=ATG;rbs_motif=GGAG/GAGG;rbs_spacer=" \
-               "5-10bp;gc_cont=0.451\n" \
-               "MFGAIAGGIASALAGGAMSKLFGGGQKAASGGIQGDVLATDNNTVGMGDAGIKSAIQGSN\n" \
-               "VPNPDEAAPSFVSGAMAKAGKGLLEGTLQAGTSAVSDKLLDLVGLGGKSAADKGKDTRDY\n" \
-               "LAAAFPELNAWERAGADASSAGMVDAGFENQKELTKMQLDNQKEIAEMQNETQKEIAGIQ\n" \
-               "SATSRQNTKDQVYAQNEMLAYQQKESTARVASIMENTNLSKQQQVSEIMRQMLTQAQTAG\n" \
-               "QYFTNDQIKEMTRKVSAEVDLVHQQTQNQRYGSSHIGATAKDISNVVTDAASGVVDIFHG\n" \
-               "IDKAVADTWNNFWKDGKADGIGSNLSRK*\n" \
-               ">NC_001422.1_7 # 3981 # 5384 # 1 # ID=1_7;partial=01;start_type=ATG;rbs_motif=GGAGG;rbs_spacer=5-10" \
-               "bp;gc_cont=0.455\n" \
-               "MVRSYYPSECHADYFDFERIEALKPAIEACGISTLSQSPMLGFHKQMDNRIKLLEEILSF\n" \
-               "RMQGVEFDNGDMYVDGHKAASDVRDEFVSVTEKLMDELAQCYNVLPQLDINNTIDHRPEG\n" \
-               "DEKWFLENEKTVTQFCRKLAAERPLKDIRDEYNYPKKKGIKDECSRLLEASTMKSRRGFA\n" \
-               "IQRLMNAMRQAHADGWFIVFDTLTLADDRLEAFYDNPNALRDYFRDIGRMVLAAEGRKAN\n" \
-               "DSHADCYQYFCVPEYGTANGRLHFHAVHFMRTLPTGSVDPNFGRRVRNRRQLNSLQNTWP\n" \
-               "YGYSMPIAVRYTQDAFSRSGWLWPVDAKGEPLKATSYMAVGFYVAKYVNKKSDMDLAAKG\n" \
-               "LGAKEWNNSLKTKLSLLPKKLFRIRMSRNFGMKMLTMTNLSTECLIQLTKLGYDATPFNQ\n" \
-               "ILKQNAKREMRLRLGKVTVADVLAAQPVTTNLLKFMRASIKMIGVSNL\n"
+    phix_seq = (
+        ">NC_001422.1_1 # 51 # 221 # 1 # ID=1_1;partial=00;start_type=ATG;rbs_motif=GGA/GAG/AGG;rbs_spacer=5-"
+        "10bp;gc_cont=0.404\n"
+        "MSRKIILIKQELLLLVYELNRSGLLAENEKIRPILAQLEKLLLCDLSPSTNDSVKN*\n"
+        ">NC_001422.1_2 # 390 # 848 # 1 # ID=1_2;partial=00;start_type=ATG;rbs_motif=GGA/GAG/AGG;rbs_spacer="
+        "5-10bp;gc_cont=0.468\n"
+        "MSQVTEQSVRFQTALASIKLIQASAVLDLTEDDFDFLTSNKVWIATDRSRARRCVEACVY\n"
+        "GTLDFVGYPRFPAPVEFIAAVIAYYVHPVNIQTACLIMEGAEFTENIINGVERPVKAAEL\n"
+        "FAFTLRVRAGNTDVLTDAEENVRQKLRAEGVM*\n"
+        ">NC_001422.1_3 # 848 # 964 # 1 # ID=1_3;partial=00;start_type=ATG;rbs_motif=AGGAG;rbs_spacer=5-10bp;"
+        "gc_cont=0.496\n"
+        "MSKGKKRSGARPGRPQPLRGTKGKRKGARLWYVGGQQF*\n"
+        ">NC_001422.1_4 # 1001 # 2284 # 1 # ID=1_4;partial=00;start_type=ATG;rbs_motif=GGAG/GAGG;rbs_spacer=5"
+        "-10bp;gc_cont=0.449\n"
+        "MSNIQTGAERMPHDLSHLGFLAGQIGRLITISTTPVIAGDSFEMDAVGALRLSPLRRGLA\n"
+        "IDSTVDIFTFYVPHRHVYGEQWIKFMKDGVNATPLPTVNTTGYIDHAAFLGTINPDTNKI\n"
+        "PKHLFQGYLNIYNNYFKAPWMPDRTEANPNELNQDDARYGFRCCHLKNIWTAPLPPETEL\n"
+        "SRQMTTSTTSIDIMGLQAAYANLHTDQERDYFMQRYHDVISSFGGKTSYDADNRPLLVMR\n"
+        "SNLWASGYDVDGTDQTSLGQFSGRVQQTYKHSVPRFFVPEHGTMFTLALVRFPPTATKEI\n"
+        "QYLNAKGALTYTDIAGDPVLYGNLPPREISMKDVFRSGDSSKKFKIAEGQWYRYAPSYVS\n"
+        "PAYHLLEGFPFIQEPPSGDLQERVLIRHHDYDQCFQSVQLLQWNSQVKFNVTVYRNLPTT\n"
+        "RDSIMTS*\n"
+        ">NC_001422.1_5 # 2395 # 2922 # 1 # ID=1_5;partial=00;start_type=ATG;rbs_motif=AGGAG;rbs_spacer=5-10"
+        "bp;gc_cont=0.420\n"
+        "MFQTFISRHNSNFFSDKLVLTSVTPASSAPVLQTPKATSSTLYFDSLTVNAGNGGFLHCI\n"
+        "QMDTSVNAANQVVSVGADIAFDADPKFFACLVRFESSSVPTTLPTAYDVYPLNGRHDGGY\n"
+        "YTVKDCVTIDVLPRTPGNNVYVGFMVWSNFTATKCRGLVSLNQVIKEIICLQPLK*\n"
+        ">NC_001422.1_6 # 2931 # 3917 # 1 # ID=1_6;partial=00;start_type=ATG;rbs_motif=GGAG/GAGG;rbs_spacer="
+        "5-10bp;gc_cont=0.451\n"
+        "MFGAIAGGIASALAGGAMSKLFGGGQKAASGGIQGDVLATDNNTVGMGDAGIKSAIQGSN\n"
+        "VPNPDEAAPSFVSGAMAKAGKGLLEGTLQAGTSAVSDKLLDLVGLGGKSAADKGKDTRDY\n"
+        "LAAAFPELNAWERAGADASSAGMVDAGFENQKELTKMQLDNQKEIAEMQNETQKEIAGIQ\n"
+        "SATSRQNTKDQVYAQNEMLAYQQKESTARVASIMENTNLSKQQQVSEIMRQMLTQAQTAG\n"
+        "QYFTNDQIKEMTRKVSAEVDLVHQQTQNQRYGSSHIGATAKDISNVVTDAASGVVDIFHG\n"
+        "IDKAVADTWNNFWKDGKADGIGSNLSRK*\n"
+        ">NC_001422.1_7 # 3981 # 5384 # 1 # ID=1_7;partial=01;start_type=ATG;rbs_motif=GGAGG;rbs_spacer=5-10"
+        "bp;gc_cont=0.455\n"
+        "MVRSYYPSECHADYFDFERIEALKPAIEACGISTLSQSPMLGFHKQMDNRIKLLEEILSF\n"
+        "RMQGVEFDNGDMYVDGHKAASDVRDEFVSVTEKLMDELAQCYNVLPQLDINNTIDHRPEG\n"
+        "DEKWFLENEKTVTQFCRKLAAERPLKDIRDEYNYPKKKGIKDECSRLLEASTMKSRRGFA\n"
+        "IQRLMNAMRQAHADGWFIVFDTLTLADDRLEAFYDNPNALRDYFRDIGRMVLAAEGRKAN\n"
+        "DSHADCYQYFCVPEYGTANGRLHFHAVHFMRTLPTGSVDPNFGRRVRNRRQLNSLQNTWP\n"
+        "YGYSMPIAVRYTQDAFSRSGWLWPVDAKGEPLKATSYMAVGFYVAKYVNKKSDMDLAAKG\n"
+        "LGAKEWNNSLKTKLSLLPKKLFRIRMSRNFGMKMLTMTNLSTECLIQLTKLGYDATPFNQ\n"
+        "ILKQNAKREMRLRLGKVTVADVLAAQPVTTNLLKFMRASIKMIGVSNL\n"
+    )
     return StringIO(phix_seq)
 
 
@@ -272,142 +337,223 @@ def test_get_gene_data(phix_prodigal_genes):
 
 
 def test_get_unannotated(phix_proteins):
-    annotated_genes = ['NP_040704.1', 'NP_040703.1', 'NP_040713.1', 'NP_040712.1', 'NP_040711.1', 'NP_040710.1',
-                       'NP_040709.1', 'NP_040707.1']
-    unannotated_genes = ['NP_040705.1', 'NP_040706.1', 'NP_040708.1']
+    annotated_genes = [
+        "NP_040704.1",
+        "NP_040703.1",
+        "NP_040713.1",
+        "NP_040712.1",
+        "NP_040711.1",
+        "NP_040710.1",
+        "NP_040709.1",
+        "NP_040707.1",
+    ]
+    unannotated_genes = ["NP_040705.1", "NP_040706.1", "NP_040708.1"]
     test_unannotated_genes = get_unannotated(phix_proteins, annotated_genes)
     assert set(unannotated_genes) == set(test_unannotated_genes)
 
 
 def test_assign_grades():
-    annotations_data = [[True, 'K00001', False, 'KOER09234OK', ['PF00001']],
-                        [False, 'K00002', True, 'KLODKJFSO234KL', ['PF01234']],
-                        [False, 'K00003', False, 'EIORWU234KLKDS', np.NaN],
-                        [False, np.NaN, False, np.NaN, np.NaN],
-                        [False, np.NaN, False, np.NaN, ['PF01235']]]
-    annotations = pd.DataFrame(annotations_data, index=['gene1', 'gene2', 'gene3', 'gene4', 'gene5'],
-                               columns=['kegg_RBH', 'kegg_hit', 'uniref_RBH', 'uniref_hit', 'pfam_hits'])
+    annotations_data = [
+        [True, "K00001", False, "KOER09234OK", ["PF00001"]],
+        [False, "K00002", True, "KLODKJFSO234KL", ["PF01234"]],
+        [False, "K00003", False, "EIORWU234KLKDS", np.NaN],
+        [False, np.NaN, False, np.NaN, np.NaN],
+        [False, np.NaN, False, np.NaN, ["PF01235"]],
+    ]
+    annotations = pd.DataFrame(
+        annotations_data,
+        index=["gene1", "gene2", "gene3", "gene4", "gene5"],
+        columns=["kegg_RBH", "kegg_hit", "uniref_RBH", "uniref_hit", "pfam_hits"],
+    )
     test_grades = assign_grades(annotations)
-    assert test_grades.loc['gene1', 'rank'] == 'A'
-    assert test_grades.loc['gene2', 'rank'] == 'B'
-    assert test_grades.loc['gene3', 'rank'] == 'C'
-    assert test_grades.loc['gene4', 'rank'] == 'E'
-    assert test_grades.loc['gene5', 'rank'] == 'D'
+    assert test_grades.loc["gene1", "rank"] == "A"
+    assert test_grades.loc["gene2", "rank"] == "B"
+    assert test_grades.loc["gene3", "rank"] == "C"
+    assert test_grades.loc["gene4", "rank"] == "E"
+    assert test_grades.loc["gene5", "rank"] == "D"
     # test no uniref
-    annotations_data2 = [[True, 'K00001', ['PF00001']],
-                         [False, 'K00002', ['PF01234']],
-                         [False, 'K00003', np.NaN],
-                         [False, np.NaN, np.NaN],
-                         [False, np.NaN, ['PF01235']]]
-    annotations2 = pd.DataFrame(annotations_data2, index=['gene1', 'gene2', 'gene3', 'gene4', 'gene5'],
-                                columns=['kegg_RBH', 'kegg_hit', 'pfam_hits'])
+    annotations_data2 = [
+        [True, "K00001", ["PF00001"]],
+        [False, "K00002", ["PF01234"]],
+        [False, "K00003", np.NaN],
+        [False, np.NaN, np.NaN],
+        [False, np.NaN, ["PF01235"]],
+    ]
+    annotations2 = pd.DataFrame(
+        annotations_data2,
+        index=["gene1", "gene2", "gene3", "gene4", "gene5"],
+        columns=["kegg_RBH", "kegg_hit", "pfam_hits"],
+    )
     test_grades2 = assign_grades(annotations2)
-    assert test_grades2.loc['gene1', 'rank'] == 'A'
-    assert test_grades2.loc['gene2', 'rank'] == 'C'
-    assert test_grades2.loc['gene3', 'rank'] == 'C'
-    assert test_grades2.loc['gene4', 'rank'] == 'E'
-    assert test_grades2.loc['gene5', 'rank'] == 'D'
+    assert test_grades2.loc["gene1", "rank"] == "A"
+    assert test_grades2.loc["gene2", "rank"] == "C"
+    assert test_grades2.loc["gene3", "rank"] == "C"
+    assert test_grades2.loc["gene4", "rank"] == "E"
+    assert test_grades2.loc["gene5", "rank"] == "D"
 
 
 @pytest.fixture()
 def phix_annotations():
-    return pd.DataFrame([['A', 'K1', 'U1', None, 'a_bug1'],
-                         ['B', 'K2', 'U2', 'P2', 'a_bug2'],
-                         ['C', None, 'U3', 'P3', 'a_bug3'],
-                         ['D', 'K4', 'U4', 'P4', 'a_bug4'],
-                         ['A', 'K5', 'U5', 'P5', 'a_bug5'],
-                         ['E', 'K6', 'U6', 'P6', 'a_bug6'],
-                         ['C', 'K7', 'U7', 'P7', 'a_bug7']],
-                        index=['NC_001422.1_1', 'NC_001422.1_2', 'NC_001422.1_3', 'NC_001422.1_4', 'NC_001422.1_5',
-                               'NC_001422.1_6', 'NC_001422.1_7'],
-                        columns=['rank', 'kegg_hit', 'uniref_hit', 'pfam_hits', 'bin_taxonomy'])
+    return pd.DataFrame(
+        [
+            ["A", "K1", "U1", None, "a_bug1"],
+            ["B", "K2", "U2", "P2", "a_bug2"],
+            ["C", None, "U3", "P3", "a_bug3"],
+            ["D", "K4", "U4", "P4", "a_bug4"],
+            ["A", "K5", "U5", "P5", "a_bug5"],
+            ["E", "K6", "U6", "P6", "a_bug6"],
+            ["C", "K7", "U7", "P7", "a_bug7"],
+        ],
+        index=[
+            "NC_001422.1_1",
+            "NC_001422.1_2",
+            "NC_001422.1_3",
+            "NC_001422.1_4",
+            "NC_001422.1_5",
+            "NC_001422.1_6",
+            "NC_001422.1_7",
+        ],
+        columns=["rank", "kegg_hit", "uniref_hit", "pfam_hits", "bin_taxonomy"],
+    )
 
 
 def test_generate_annotated_fasta_short(phix_prodigal_genes, phix_annotations):
-    fasta_generator_short = generate_annotated_fasta(phix_prodigal_genes, phix_annotations, verbosity='short',
-                                                     name='phiX')
-    short_fasta_header_dict = {seq.metadata['id']: seq.metadata['description'] for seq in fasta_generator_short}
+    fasta_generator_short = generate_annotated_fasta(
+        phix_prodigal_genes, phix_annotations, verbosity="short", name="phiX"
+    )
+    short_fasta_header_dict = {
+        seq.metadata["id"]: seq.metadata["description"] for seq in fasta_generator_short
+    }
 
-    assert short_fasta_header_dict['phiX_NC_001422.1_1'] == 'rank: A; K1 (db=kegg)'
-    assert short_fasta_header_dict['phiX_NC_001422.1_2'] == 'rank: B; U2 (db=uniref)'
-    assert short_fasta_header_dict['phiX_NC_001422.1_4'] == 'rank: D; P4 (db=pfam)'
+    assert short_fasta_header_dict["phiX_NC_001422.1_1"] == "rank: A; K1 (db=kegg)"
+    assert short_fasta_header_dict["phiX_NC_001422.1_2"] == "rank: B; U2 (db=uniref)"
+    assert short_fasta_header_dict["phiX_NC_001422.1_4"] == "rank: D; P4 (db=pfam)"
 
 
 @pytest.fixture()
 def phix_annotations_no_kegg():
-    return pd.DataFrame([['B', 'U1', None, 'a_bug1'],
-                         ['B', 'U2', 'P2', 'a_bug2'],
-                         ['C', 'U3', 'P3', 'a_bug3'],
-                         ['D', 'U4', 'P4', 'a_bug4'],
-                         ['B', 'U5', 'P5', 'a_bug5'],
-                         ['E', 'U6', 'P6', 'a_bug6'],
-                         ['C', 'U7', 'P7', 'a_bug7']],
-                        index=['NC_001422.1_1', 'NC_001422.1_2', 'NC_001422.1_3', 'NC_001422.1_4', 'NC_001422.1_5',
-                               'NC_001422.1_6', 'NC_001422.1_7'],
-                        columns=['rank', 'uniref_hit', 'pfam_hits', 'bin_taxonomy'])
-
-
-def test_generate_annotated_fasta_short_no_kegg(phix_prodigal_genes, phix_annotations_no_kegg):
+    return pd.DataFrame(
+        [
+            ["B", "U1", None, "a_bug1"],
+            ["B", "U2", "P2", "a_bug2"],
+            ["C", "U3", "P3", "a_bug3"],
+            ["D", "U4", "P4", "a_bug4"],
+            ["B", "U5", "P5", "a_bug5"],
+            ["E", "U6", "P6", "a_bug6"],
+            ["C", "U7", "P7", "a_bug7"],
+        ],
+        index=[
+            "NC_001422.1_1",
+            "NC_001422.1_2",
+            "NC_001422.1_3",
+            "NC_001422.1_4",
+            "NC_001422.1_5",
+            "NC_001422.1_6",
+            "NC_001422.1_7",
+        ],
+        columns=["rank", "uniref_hit", "pfam_hits", "bin_taxonomy"],
+    )
+
+
+def test_generate_annotated_fasta_short_no_kegg(
+    phix_prodigal_genes, phix_annotations_no_kegg
+):
     # drop kegg hit column
-    fasta_generator_short_no_kegg = generate_annotated_fasta(phix_prodigal_genes, phix_annotations_no_kegg,
-                                                             verbosity='short', name='phiX')
-    short_fasta_header_dict_no_kegg = {seq.metadata['id']: seq.metadata['description']
-                                       for seq in fasta_generator_short_no_kegg}
+    fasta_generator_short_no_kegg = generate_annotated_fasta(
+        phix_prodigal_genes, phix_annotations_no_kegg, verbosity="short", name="phiX"
+    )
+    short_fasta_header_dict_no_kegg = {
+        seq.metadata["id"]: seq.metadata["description"]
+        for seq in fasta_generator_short_no_kegg
+    }
     print(short_fasta_header_dict_no_kegg)
 
-    assert short_fasta_header_dict_no_kegg['phiX_NC_001422.1_1'] == 'rank: B; U1 (db=uniref)'
-    assert short_fasta_header_dict_no_kegg['phiX_NC_001422.1_2'] == 'rank: B; U2 (db=uniref)'
-    assert short_fasta_header_dict_no_kegg['phiX_NC_001422.1_4'] == 'rank: D; P4 (db=pfam)'
+    assert (
+        short_fasta_header_dict_no_kegg["phiX_NC_001422.1_1"]
+        == "rank: B; U1 (db=uniref)"
+    )
+    assert (
+        short_fasta_header_dict_no_kegg["phiX_NC_001422.1_2"]
+        == "rank: B; U2 (db=uniref)"
+    )
+    assert (
+        short_fasta_header_dict_no_kegg["phiX_NC_001422.1_4"] == "rank: D; P4 (db=pfam)"
+    )
 
 
 def test_generate_annotated_fasta_long(phix_prodigal_genes, phix_annotations):
-    fasta_generator_long = generate_annotated_fasta(phix_prodigal_genes, phix_annotations, verbosity='long',
-                                                    name='phiX')
-    long_fasta_header_dict = {seq.metadata['id']: seq.metadata['description'] for seq in fasta_generator_long}
-    assert long_fasta_header_dict['phiX_NC_001422.1_1'] == 'rank: A; K1 (db=kegg); U1 (db=uniref); a_bug1'
+    fasta_generator_long = generate_annotated_fasta(
+        phix_prodigal_genes, phix_annotations, verbosity="long", name="phiX"
+    )
+    long_fasta_header_dict = {
+        seq.metadata["id"]: seq.metadata["description"] for seq in fasta_generator_long
+    }
+    assert (
+        long_fasta_header_dict["phiX_NC_001422.1_1"]
+        == "rank: A; K1 (db=kegg); U1 (db=uniref); a_bug1"
+    )
 
 
 def test_create_annotated_fasta(phix_prodigal_genes, phix_annotations, tmpdir):
-    filt_fasta = tmpdir.mkdir('test_annotate_fasta').join('annotated_fasta.faa')
+    filt_fasta = tmpdir.mkdir("test_annotate_fasta").join("annotated_fasta.faa")
     create_annotated_fasta(phix_prodigal_genes, phix_annotations, str(filt_fasta))
     assert os.path.isfile(filt_fasta)
 
 
 def test_generate_renamed_fasta(fasta_loc):
-    genome_fasta_header_dict = {seq.metadata['id']: seq.metadata['description']
-                                for seq in generate_renamed_fasta(fasta_loc, 'phiX')}
+    genome_fasta_header_dict = {
+        seq.metadata["id"]: seq.metadata["description"]
+        for seq in generate_renamed_fasta(fasta_loc, "phiX")
+    }
     assert len(genome_fasta_header_dict) == 1
-    assert genome_fasta_header_dict['phiX_NC_001422.1'] == 'Coliphage phi-X174, complete genome'
+    assert (
+        genome_fasta_header_dict["phiX_NC_001422.1"]
+        == "Coliphage phi-X174, complete genome"
+    )
 
 
 def test_rename_fasta(fasta_loc, tmpdir):
-    filt_fasta = tmpdir.mkdir('test_annotate_scaffolds').join('annotated_fasta.fasta')
-    rename_fasta(fasta_loc, str(filt_fasta), 'phiX')
+    filt_fasta = tmpdir.mkdir("test_annotate_scaffolds").join("annotated_fasta.fasta")
+    rename_fasta(fasta_loc, str(filt_fasta), "phiX")
     assert os.path.isfile(filt_fasta)
 
 
 def test_run_trna_scan(tmpdir, logger):
-    filt_fasta = tmpdir.mkdir('test_trnascan1')
-    no_trna = run_trna_scan(os.path.join('tests', 'data', 'e_coli_16S.fasta'),
-                            str(filt_fasta), 'no_trnas', logger, threads=1, verbose=False)
+    filt_fasta = tmpdir.mkdir("test_trnascan1")
+    no_trna = run_trna_scan(
+        os.path.join("tests", "data", "e_coli_16S.fasta"),
+        str(filt_fasta),
+        "no_trnas",
+        logger,
+        threads=1,
+        verbose=False,
+    )
     assert no_trna is None
 
-    filt_fasta = tmpdir.mkdir('test_trnascan2')
-    trnas_loc = os.path.join('tests', 'data', 'trnas.fa')
-    trnas = run_trna_scan(trnas_loc, str(filt_fasta), 'phiX', logger, threads=1, verbose=False)
+    filt_fasta = tmpdir.mkdir("test_trnascan2")
+    trnas_loc = os.path.join("tests", "data", "trnas.fa")
+    trnas = run_trna_scan(
+        trnas_loc, str(filt_fasta), "phiX", logger, threads=1, verbose=False
+    )
     assert trnas.shape == (6, 9)
 
 
 def test_run_barrnap(fasta_loc, logger):
-    no_rrnas = run_barrnap(fasta_loc, 'phiX', logger, threads=1, verbose=False)
+    no_rrnas = run_barrnap(fasta_loc, "phiX", logger, threads=1, verbose=False)
     assert no_rrnas is None
-    rrna_table = run_barrnap(os.path.join('tests', 'data', 'e_coli_16S.fasta'), 'coli', logger, threads=1, verbose=False)
+    rrna_table = run_barrnap(
+        os.path.join("tests", "data", "e_coli_16S.fasta"),
+        "coli",
+        logger,
+        threads=1,
+        verbose=False,
+    )
     assert rrna_table.shape == (1, 8)
-    assert rrna_table.loc[0, 'type'] == '16S rRNA'
-    assert rrna_table.loc[0, 'fasta'] == 'coli'
-    assert rrna_table.loc[0, 'begin'] == 100
-    assert rrna_table.loc[0, 'end'] == 1637
-    assert rrna_table.loc[0, 'strand'] == '-'
+    assert rrna_table.loc[0, "type"] == "16S rRNA"
+    assert rrna_table.loc[0, "fasta"] == "coli"
+    assert rrna_table.loc[0, "begin"] == 100
+    assert rrna_table.loc[0, "end"] == 1637
+    assert rrna_table.loc[0, "strand"] == "-"
 
 
 class FakeDatabaseHandler:
@@ -418,126 +564,158 @@ def get_database_names():
 
 def test_do_blast_style_search(mmseqs_db, target_mmseqs_db, tmpdir, logger):
     database_handler = FakeDatabaseHandler()
-    working_dir = tmpdir.mkdir('test_blast_search')
-    get_fake_description = partial(get_basic_description, db_name='fake')
-    do_blast_style_search(mmseqs_db, target_mmseqs_db, str(working_dir),
-                          database_handler, get_fake_description, logger,
-                          db_name='fake', threads=1)
+    working_dir = tmpdir.mkdir("test_blast_search")
+    get_fake_description = partial(get_basic_description, db_name="fake")
+    do_blast_style_search(
+        mmseqs_db,
+        target_mmseqs_db,
+        str(working_dir),
+        database_handler,
+        get_fake_description,
+        logger,
+        db_name="fake",
+        threads=1,
+    )
 
 
 def test_count_motifs(phix_proteins):
-    motif_counts = count_motifs(phix_proteins, motif='(A.A)')
+    motif_counts = count_motifs(phix_proteins, motif="(A.A)")
     assert len(motif_counts) == 11
-    assert motif_counts['NP_040713.1'] == 7
-    assert motif_counts['NP_040709.1'] == 0
+    assert motif_counts["NP_040713.1"] == 7
+    assert motif_counts["NP_040709.1"] == 0
 
 
 def test_strip_endings():
-    assert strip_endings('abc.efg', ['.efg', '.jkl']) == 'abc'
-    assert strip_endings('abc.jkl', ['.efg', '.jkl']) == 'abc'
-    assert strip_endings('123456', ['.efg', '.jkl']) == '123456'
+    assert strip_endings("abc.efg", [".efg", ".jkl"]) == "abc"
+    assert strip_endings("abc.jkl", [".efg", ".jkl"]) == "abc"
+    assert strip_endings("123456", [".efg", ".jkl"]) == "123456"
 
 
 def test_process_custom_dbs(phix_proteins, tmpdir, logger):
-    custom_db_dir = tmpdir.mkdir('custom_dbs')
-    process_custom_dbs([phix_proteins], ['phix'], os.path.join(custom_db_dir, 'custom_dbs0'), logger)
-    assert os.path.isfile(os.path.join(custom_db_dir, 'custom_dbs0', 'phix.custom.mmsdb'))
-    process_custom_dbs(None, None, os.path.join(custom_db_dir, 'custom_dbs1'), logger)
-    assert len(os.listdir(os.path.join(custom_db_dir, 'custom_dbs1'))) == 0
+    custom_db_dir = tmpdir.mkdir("custom_dbs")
+    process_custom_dbs(
+        [phix_proteins], ["phix"], os.path.join(custom_db_dir, "custom_dbs0"), logger
+    )
+    assert os.path.isfile(
+        os.path.join(custom_db_dir, "custom_dbs0", "phix.custom.mmsdb")
+    )
+    process_custom_dbs(None, None, os.path.join(custom_db_dir, "custom_dbs1"), logger)
+    assert len(os.listdir(os.path.join(custom_db_dir, "custom_dbs1"))) == 0
     with pytest.raises(ValueError):
-        process_custom_dbs(['thing1', 'thing2'], ['thing1'], os.path.join(custom_db_dir, 'custom_dbs2'), logger)
+        process_custom_dbs(
+            ["thing1", "thing2"],
+            ["thing1"],
+            os.path.join(custom_db_dir, "custom_dbs2"),
+            logger,
+        )
 
 
 def test_get_dubs():
     w_dups = [True, False, True, True]
     assert get_dups(w_dups) == [True, True, False, False]
-    no_dups = ['a', 'b', 'c']
+    no_dups = ["a", "b", "c"]
     assert get_dups(no_dups) == [True, True, True]
 
 
 def test_parse_hmmsearch_domtblout():
-    parsed_hit = parse_hmmsearch_domtblout(os.path.join('tests', 'data', 'hmmsearch_hit.txt'))
+    parsed_hit = parse_hmmsearch_domtblout(
+        os.path.join("tests", "data", "hmmsearch_hit.txt")
+    )
     assert parsed_hit.shape == (1, 23)
-    assert parsed_hit.loc[0, 'query_id'] == 'NP_040710.1'
-    assert parsed_hit.loc[0, 'query_length'] == 38
-    assert parsed_hit.loc[0, 'target_id'] == 'Microvir_J'
-    assert parsed_hit.loc[0, 'target_ascession'] == 'PF04726.13'
-    assert parsed_hit.loc[0, 'full_evalue'] == 6.900000e-31
+    assert parsed_hit.loc[0, "query_id"] == "NP_040710.1"
+    assert parsed_hit.loc[0, "query_length"] == 38
+    assert parsed_hit.loc[0, "target_id"] == "Microvir_J"
+    assert parsed_hit.loc[0, "target_ascession"] == "PF04726.13"
+    assert parsed_hit.loc[0, "full_evalue"] == 6.900000e-31
 
 
 @pytest.fixture()
 def annotated_fake_gff_loc():
-    return os.path.join('tests', 'data', 'annotated_fake_gff.gff')
+    return os.path.join("tests", "data", "annotated_fake_gff.gff")
 
 
 @pytest.fixture()
 def fake_phix_annotations():
-    return pd.DataFrame([[np.NaN],
-                         ['GH13'],
-                         [np.NaN],
-                         [np.NaN],
-                         [np.NaN],
-                         [np.NaN],
-                         [np.NaN]],
-                        index=['NC_001422.1_1', 'NC_001422.1_2', 'NC_001422.1_3', 'NC_001422.1_4',
-                               'NC_001422.1_5', 'NC_001422.1_6', 'NC_001422.1_7'],
-                        columns=['cazy_id'])
+    return pd.DataFrame(
+        [[np.NaN], ["GH13"], [np.NaN], [np.NaN], [np.NaN], [np.NaN], [np.NaN]],
+        index=[
+            "NC_001422.1_1",
+            "NC_001422.1_2",
+            "NC_001422.1_3",
+            "NC_001422.1_4",
+            "NC_001422.1_5",
+            "NC_001422.1_6",
+            "NC_001422.1_7",
+        ],
+        columns=["cazy_id"],
+    )
 
 
 @pytest.fixture()
 def fake_gff_loc():
-    return os.path.join('tests', 'data', 'fake_gff.gff')
+    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')
-    annotate_gff(fake_gff_loc, test_annotated_gff_loc, fake_phix_annotations, 'fake')
+# 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")
+    annotate_gff(fake_gff_loc, test_annotated_gff_loc, fake_phix_annotations, "fake")
     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')
+    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'),
+        hmm_info_path=os.path.join("tests", "data", "hmm_thresholds.txt"),
         top_hit=True,
-        use_dbcan2_thresholds=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"
+    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')
+    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'),
+        hmm_info_path=os.path.join("tests", "data", "hmm_thresholds.txt"),
         top_hit=True,
-        use_dbcan2_thresholds=False)
+        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_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')
+    # 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', logger=logger)
+    output_rcvd = vogdb_hmmscan_formater(hits=hits, db_name="vogdb", logger=logger)
     output_rcvd.sort_index(inplace=True)
     output_expt.sort_index(inplace=True)
     assert output_rcvd.equals(output_expt), "Error in vogdb_hmmscan_formater"
@@ -545,19 +723,25 @@ def test_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", "GT5.hmm"],
-                           "bin_1.scaffold_2": ["GH1", "GH1.hmm"]}, index=["cazy_ids", "cazy_best_hit"]).T
-    input_b6 = os.path.join('tests', 'data', 'unformatted_cazy.b6')
+    output_expt = pd.DataFrame(
+        {
+            "bin_1.scaffold_1": ["GT4; GT5", "GT5.hmm"],
+            "bin_1.scaffold_2": ["GH1", "GH1.hmm"],
+        },
+        index=["cazy_ids", "cazy_best_hit"],
+    ).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 = 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
+test_gbk = (
+    """LOCUS       NC_001422.1   5386 bp   DNA   linear   ENV   %s
 FEATURES             Location/Qualifiers
      CDS             51..221
                      /gene=1_1
@@ -686,52 +870,96 @@ def test_dbcan_hmmscan_formater():
      5281 tgggaaaagt tactgtagcc gacgttttgg cggcgcaacc tgtgacgaca aatctgctca
      5341 aatttatgcg cgcttcgata aaaatgattg gcgtatccaa cctgca
 //
-""" % time.strftime('%d-%b-%Y').upper()
+"""
+    % time.strftime("%d-%b-%Y").upper()
+)
 
 
 # TODO: Make better, add features that should carry over with the right names
 def test_make_gbk_from_gff_and_fasta(fake_gff_loc, fasta_loc, phix_proteins, tmpdir):
     gbk = make_gbk_from_gff_and_fasta(fake_gff_loc, fasta_loc, phix_proteins)
     assert gbk == test_gbk
-    gbk_test_file = os.path.join(tmpdir.mkdir('gbk_test'), 'test.gbk')
+    gbk_test_file = os.path.join(tmpdir.mkdir("gbk_test"), "test.gbk")
     make_gbk_from_gff_and_fasta(fake_gff_loc, fasta_loc, phix_proteins, gbk_test_file)
     assert os.path.isfile(gbk_test_file)
 
 
 def test_make_trnas_interval():
-    test_scaffold = 'scaffold_1'
+    test_scaffold = "scaffold_1"
     i = 1
 
-    test_row1 = {'Begin': 1, 'End': 13, 'Score': 1000, 'Codon': 'AUG', 'Type': 'AUG codon coding tRNA',
-                 'Note': np.NaN}
+    test_row1 = {
+        "Begin": 1,
+        "End": 13,
+        "Score": 1000,
+        "Codon": "AUG",
+        "Type": "AUG codon coding tRNA",
+        "Note": np.NaN,
+    }
     begin1, end1, metadata1 = make_trnas_interval(test_scaffold, test_row1, i)
     assert begin1 == 1
     assert end1 == 13
-    assert metadata1 == {'source': 'tRNAscan-SE', 'type': 'tRNA', 'score': 1000, 'strand': '+', 'phase': 0,
-                         'ID': 'scaffold_1_tRNA_1', 'codon': 'AUG', 'product': 'tRNA-AUG codon coding tRNA'}
-
-    test_row2 = {'End': 1, 'Begin': 13, 'Score': 1000, 'Codon': 'AUG', 'Type': 'AUG codon coding tRNA',
-                 'Note': 'a note'}
+    assert metadata1 == {
+        "source": "tRNAscan-SE",
+        "type": "tRNA",
+        "score": 1000,
+        "strand": "+",
+        "phase": 0,
+        "ID": "scaffold_1_tRNA_1",
+        "codon": "AUG",
+        "product": "tRNA-AUG codon coding tRNA",
+    }
+
+    test_row2 = {
+        "End": 1,
+        "Begin": 13,
+        "Score": 1000,
+        "Codon": "AUG",
+        "Type": "AUG codon coding tRNA",
+        "Note": "a note",
+    }
     begin2, end2, metadata2 = make_trnas_interval(test_scaffold, test_row2, i)
     assert begin2 == 1
     assert end2 == 13
-    assert metadata2 == {'source': 'tRNAscan-SE', 'type': 'tRNA', 'score': 1000, 'strand': '-', 'phase': 0,
-                         'ID': 'scaffold_1_tRNA_1', 'codon': 'AUG', 'product': 'tRNA-AUG codon coding tRNA',
-                         'Note': 'a note'}
+    assert metadata2 == {
+        "source": "tRNAscan-SE",
+        "type": "tRNA",
+        "score": 1000,
+        "strand": "-",
+        "phase": 0,
+        "ID": "scaffold_1_tRNA_1",
+        "codon": "AUG",
+        "product": "tRNA-AUG codon coding tRNA",
+        "Note": "a note",
+    }
 
 
 def test_make_rrnas_interval():
-    test_scaffold = 'scaffold_1'
+    test_scaffold = "scaffold_1"
     i = 1
 
-    test_row1 = {'begin': 1, 'end': 13, 'e-value': 1000, 'strand': '+', 'type': '16S ribosomal rRNA gene',
-                 'note': 'its a pseudo'}
+    test_row1 = {
+        "begin": 1,
+        "end": 13,
+        "e-value": 1000,
+        "strand": "+",
+        "type": "16S ribosomal rRNA gene",
+        "note": "its a pseudo",
+    }
     begin1, end1, metadata1 = make_rrnas_interval(test_scaffold, test_row1, i)
     assert begin1 == 1
     assert end1 == 13
-    assert metadata1 == {'source': 'barrnap', 'type': 'rRNA', 'score': 1000, 'strand': '+', 'phase': 0,
-                         'ID': 'scaffold_1_rRNA_1', 'product': '16S ribosomal RNA',
-                         'gene': '16S ribosomal rRNA gene', 'Note': 'its a pseudo'}
+    assert metadata1 == {
+        "source": "barrnap",
+        "type": "rRNA",
+        "score": 1000,
+        "strand": "+",
+        "phase": 0,
+        "ID": "scaffold_1_rRNA_1",
+        "product": "16S ribosomal RNA",
+        "gene": "16S ribosomal rRNA gene",
+        "Note": "its a pseudo",
+    }
 
 
 annotated_fake_gff_w_rna = """##gff-version 3
@@ -747,16 +975,34 @@ def test_make_rrnas_interval():
 
 
 def test_add_intervals_to_gff(annotated_fake_gff_loc, tmpdir):
-    add_intervals_test_loc = tmpdir.mkdir('fake_rrnas_loc')
-    annotate_fake_gff_loc_w_rna = os.path.join(add_intervals_test_loc, 'fake.gff')
+    add_intervals_test_loc = tmpdir.mkdir("fake_rrnas_loc")
+    annotate_fake_gff_loc_w_rna = os.path.join(add_intervals_test_loc, "fake.gff")
     copy(annotated_fake_gff_loc, annotate_fake_gff_loc_w_rna)
-    fake_rrnas_loc = os.path.join(add_intervals_test_loc, 'rrnas.tsv')
-    fake_rrnas = pd.DataFrame([['fake_NC_001422.1', 990, 1000, 101.0, '+', '16S ribosomal rRNA gene', np.NaN]],
-                              columns=['scaffold', 'begin', 'end', 'e-value', 'strand', 'type', 'note'])
-    fake_rrnas.to_csv(fake_rrnas_loc, sep='\t')
-    add_intervals_to_gff(fake_rrnas_loc, annotate_fake_gff_loc_w_rna, {'fake_NC_001422.1': 6000},
-                         make_rrnas_interval, 'scaffold', logger)
+    fake_rrnas_loc = os.path.join(add_intervals_test_loc, "rrnas.tsv")
+    fake_rrnas = pd.DataFrame(
+        [
+            [
+                "fake_NC_001422.1",
+                990,
+                1000,
+                101.0,
+                "+",
+                "16S ribosomal rRNA gene",
+                np.NaN,
+            ]
+        ],
+        columns=["scaffold", "begin", "end", "e-value", "strand", "type", "note"],
+    )
+    fake_rrnas.to_csv(fake_rrnas_loc, sep="\t")
+    add_intervals_to_gff(
+        fake_rrnas_loc,
+        annotate_fake_gff_loc_w_rna,
+        {"fake_NC_001422.1": 6000},
+        make_rrnas_interval,
+        "scaffold",
+        logger,
+    )
     assert os.path.isfile(annotate_fake_gff_loc_w_rna)
-    gff = list(read_sequence(annotate_fake_gff_loc_w_rna, format='gff3'))
+    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