###############################################################################
#                                                                             #
#    This program is free software: you can redistribute it and/or modify     #
#    it under the terms of the GNU General Public License as published by     #
#    the Free Software Foundation, either version 3 of the License, or        #
#    (at your option) any later version.                                      #
#                                                                             #
#    This program is distributed in the hope that it will be useful,          #
#    but WITHOUT ANY WARRANTY; without even the implied warranty of           #
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            #
#    GNU General Public License for more details.                             #
#                                                                             #
#    You should have received a copy of the GNU General Public License        #
#    along with this program. If not, see <http://www.gnu.org/licenses/>.     #
#                                                                             #
###############################################################################

import logging
import os
import re
import sys
from collections import defaultdict, Counter

import psycopg2
from biolib.seq_io import read_fasta
from biolib.taxonomy import Taxonomy
from psycopg2.extensions import AsIs

import Config
import ConfigMetadata
from Exceptions import GenomeDatabaseError


class MetadataManager(object):

    def __init__(self, cur, currentUser):
        """Initialize.

        Parameters
        ----------
        cur : psycopg2.cursor
            Database cursor.
        currentUser : User
            Current user of database.
        """

        self.logger = logging.getLogger()

        self.cur = cur
        self.currentUser = currentUser

    def viewMetadata(self):
        '''
         Function: viewMetadata
        Lists all metadata field available in the database

        Returns:
        Print lists of Metadata on Screen
        '''
        try:
            self.cur.execute("SELECT * FROM view_list_meta_columns")
            print "\t".join(("Table", "Field", "Datatype", "Description"))
            for tup in self.cur.fetchall():
                info_list = list(tup)
                if info_list[2] == "double precision":
                    info_list[2] = "float"
                elif str(info_list[2]).startswith("timestamp"):
                    info_list[2] = "timestamp"
                print "\t".join(info_list)
        except GenomeDatabaseError as e:
            raise e

    def exportLPSNMetadata(self, path, outformat=None):
        '''
        Function: exportLPSNMetadata
        Export LPSN metadata for all LPSN ranks to either a tab-delimited or csv file.
        The metadata is generated by parsing all html files in LPSN

        :param path: Path to the output file
        :param outformat : CSV or Tab-delimited output
        '''

        try:
            query_tmp = 'CREATE TEMPORARY TABLE temp_tb as SELECT * FROM lpsn_metadata ORDER BY index,"Name";'
            query_tmp += "ALTER TABLE temp_tb DROP COLUMN index;"
            self.cur.execute(query_tmp)
            query = "SELECT * FROM temp_tb"
            if outformat == 'csv':
                outputquery = 'copy ({0}) to stdout with csv header'.format(
                    query)
            elif outformat == 'tab':
                outputquery = "copy ({0}) to stdout with DELIMITER '\t' NULL AS 'none' csv header;".format(
                    query)
            with open(path, 'w') as f:
                self.cur.copy_expert(outputquery, f)
            print "Export Successful"
        except GenomeDatabaseError as e:
            raise e

    def exportMetadata(self, path, outformat=None):
        '''
        Function: exportMetadata
        Export metadata for all genomes to either a tab-delimited or csv file

        :param path: Path to the output file
        :param outformat : CSV or Tab-delimited output
        '''

        try:
            query_tmp = "CREATE TEMPORARY TABLE temp_tb as SELECT * FROM metadata_view;"
            query_tmp += "ALTER TABLE temp_tb DROP COLUMN id;"
            self.cur.execute(query_tmp)
            query = "SELECT * FROM temp_tb"
            if outformat == 'csv':
                outputquery = 'copy ({0}) to stdout with csv header'.format(
                    query)
            elif outformat == 'tab':
                outputquery = "copy ({0}) to stdout with DELIMITER '\t' NULL AS 'none' csv header;".format(
                    query)
            with open(path, 'w') as f:
                self.cur.copy_expert(outputquery, f)
            print "Export Successful"
        except GenomeDatabaseError as e:
            raise e

    def ExportGenomePaths(self, path):
        '''
        Function: ExportGenomePaths
        Export the full path for all genomes to a csv file

        :param path: Path to the output file
        '''

        try:
            query_tmp = "SELECT id_at_source,external_id_prefix,fasta_file_location FROM genomes,genome_sources WHERE genomes.genome_source_id=genome_sources.id;"
            self.cur.execute(query_tmp)
            with open(path, "w") as f:
                for (id, prefix, file_location) in self.cur:
                    dir_prefix = None
                    if prefix == 'U':
                        dir_prefix = Config.GTDB_GENOME_USR_DIR
                    elif prefix == 'RS':
                        dir_prefix = Config.GTDB_GENOME_RSQ_DIR
                    elif prefix == 'GB':
                        dir_prefix = Config.GTDB_GENOME_GBK_DIR
                    else:
                        raise GenomeDatabaseError(
                            "Unrecognized database prefix: %s" % prefix)
                    f.write("{0}\t{1}\n".format(
                        prefix + "_" + id, os.path.dirname(os.path.join(dir_prefix, file_location))))
            print "Export Successful"
        except GenomeDatabaseError as e:
            raise e

    def exportTaxonomy(self, taxonomy_src, output_file):
        """Write taxonomy to file.

        Parameters
        ----------
        taxonomy_src : str
          Indicates desired taxonomy ('GTDB' or 'NCBI').
        output_file : str
          Output file.
        """

        try:
            if taxonomy_src == 'NCBI':
                self.cur.execute(
                    "SELECT accession, ncbi_taxonomy FROM metadata_view")
            else:
                self.cur.execute(
                    "SELECT accession, gtdb_taxonomy FROM metadata_view")

            fout = open(output_file, 'w')
            for genome_id, taxonomy in self.cur.fetchall():
                if taxonomy:
                    fout.write('%s\t%s\n' % (genome_id, taxonomy))
                else:
                    fout.write('%s\t%s\n' %
                               (genome_id, ';'.join(Taxonomy.rank_prefixes)))
            fout.close()

            print 'Taxonomy information written to: %s' % output_file
        except GenomeDatabaseError as e:
            raise e

    def exportTaxonomyMapping(self, src, dest, output_file):
        # type: (str, str, str) -> None
        """Summarises the source taxonomy mapping to the destination taxonomy (only for representative species).

        Parameters
        ----------
        :param src: str
            Indicates the source taxonomy ('GTDB' or 'SILVA')
        :param dest: str
            Indicates the destination taxonomy ('SILVA', or 'GTDB')
        :param output_file: str
            Output file.
        :return: None
        """

        valid_db = ['GTDB', 'SILVA']
        assert (src in valid_db and dest in valid_db and src != dest)

        # Get the taxonomy for each representative species
        try:
            self.cur.execute("""SELECT mt.id, mt.gtdb_domain, mt.gtdb_phylum, mt.gtdb_class, mt.gtdb_order,
                                    mt.gtdb_family, mt.gtdb_genus, mt.gtdb_species, mn.ncbi_wgs_master
                                FROM metadata_taxonomy mt
                                LEFT JOIN metadata_ncbi mn ON mn.id = mt.id
                                WHERE gtdb_representative""")
        except GenomeDatabaseError as e:
            raise e

        gtdb_taxonomy = dict()
        for id, gtdb_domain, gtdb_phylum, gtdb_class, gtdb_order, gtdb_family, gtdb_genus, gtdb_species, ncbi_wgs_master in self.cur.fetchall():
            gtdb_taxonomy[id] = {'domain': gtdb_domain, 'phylum': gtdb_phylum, 'class': gtdb_class, 'order': gtdb_order,
                             'family': gtdb_family, 'genus': gtdb_genus, 'species': gtdb_species}

        # Load the SILVA taxonomic rank information
        silva_mapping = self._parse_silva_mapping(ConfigMetadata.GTDB_SSU_SILVA_MAPPING)

        # Load the aggregate SILVA SSU blast hits
        ssu_agg = self._parse_silva_ssu_aggregate(ConfigMetadata.GTDB_SSU_SILVA_AGGREGATE)

        def _get_most_significant_hit(blast_hits, min_align_length):
            # Remove those under the threshold
            trim_set = list()
            for hit in blast_hits:
                if hit['blast_align_len'] > min_align_length:
                    trim_set.append(hit)

            # Annotate with the consensus
            tax_seen = defaultdict(list)
            for hit in trim_set:
                tax_seen[hit['taxonomy']].append(True)

            new_hits = list()
            for hit in blast_hits:
                new_hit = hit
                hit['tax_count'] = len(tax_seen[hit['taxonomy']])
                new_hits.append(new_hit)

            hits_sorted = sorted(new_hits, key=lambda k: (-k['tax_count'], -k['blast_bitscore'], -k['blast_align_len']))
            for hit in hits_sorted:
                return hit
            return None

        def _match_silva_taxonomy(input_taxonomy, silva_taxonomy):
            out = dict()

            # Try and match all of the ranks
            input_ranks = input_taxonomy.split(';')
            for i in range(1, len(input_ranks)):
                cur_tax_string = ';'.join(input_ranks[0:i])
                cur_rank = silva_taxonomy.get(cur_tax_string + ';')
                out[cur_rank] = input_ranks[i - 1]

            out['species'] = input_ranks[-1]

            return out

        def _create_output_file(a_to_b, output_path):
            with open(output_path, 'w') as f:
                f.write('%s\t%s\t%s\t%s\t%s\n' % ('tax_rank', 'tax_from', 'num_genomes', 'num_ranks', 'tax_to'))
                for rank in ['domain', 'phylum', 'class', 'order', 'family', 'genus', 'species']:
                    a_tax_dict = a_to_b[rank]
                    for a_tax, a_tax_list in a_tax_dict.items():
                        counts = Counter(a_tax_list)

                        list_props = list()
                        for cur_rank, cur_count in sorted(counts.items(), key=lambda v: -v[1]):
                            cur_prop = (float(cur_count) / len(a_tax_list)) * 100.0
                            list_props.append('%s (%.2f%%)' % (cur_rank, cur_prop))

                        prop_str = ', '.join(list_props)
                        f.write('%s\t%s\t%d\t%d\t%s\n' % (rank, a_tax, len(a_tax_list), len(counts), prop_str))

        # Variables for reporting.
        removed_due_to_eukaryota_contamination = set()
        had_no_significant_hits = set()

        # Iterate over all BLAST hits for each genome id
        gtdb_to_silva = defaultdict(lambda: defaultdict(list))
        silva_to_gtdb = defaultdict(lambda: defaultdict(list))
        for gid, blast_hits in ssu_agg.items():

            # Take the consensus taxonomy
            most_sig_hit = _get_most_significant_hit(blast_hits, min_align_length=900)

            if most_sig_hit is None:
                had_no_significant_hits.add(gid)
                continue
            if most_sig_hit['taxonomy'][0:9] == 'Eukaryota':
                removed_due_to_eukaryota_contamination.add(gid)
                continue

            # Find the corresponding taxonomy in GTDB
            cur_silva_tax = _match_silva_taxonomy(most_sig_hit['taxonomy'], silva_mapping)
            cur_gtdb_tax = gtdb_taxonomy[gid]

            # Find where the taxonomies agree/disagree.
            identity_thresholds = [('domain', 0), ('phylum', 75.0), ('class', 78.5), ('order', 82.0), ('family', 86.5),
                                   ('genus', 94.5), ('species', 98.7)]
            for cur_domain, cur_thresh in identity_thresholds:
                if most_sig_hit['blast_perc_identity'] >= cur_thresh:
                    if cur_domain not in cur_silva_tax:
                        continue
                    gtdb_to_silva[cur_domain][cur_gtdb_tax[cur_domain]].append(cur_silva_tax[cur_domain])
                    silva_to_gtdb[cur_domain][cur_silva_tax[cur_domain]].append(cur_gtdb_tax[cur_domain])

        # Create the output file
        if src == 'GTDB' and dest == 'SILVA':
            _create_output_file(gtdb_to_silva, output_file)
        else:
            _create_output_file(silva_to_gtdb, output_file)

        print('Number of genomes which were ignored due to poor BLAST hits: %d' % len(had_no_significant_hits))
        print('Number of genomes which were ignored due to Eukaryota contamination in BLAST hits: %d' % len(removed_due_to_eukaryota_contamination))


    def importMetadata(self, table=None, field=None, typemeta=None, metafile=None):
        '''
        Function importMetadata
        import one field of Metadata for a list of Genomes

        :param table: Table where the column is located
        :param field: Name of the Column
        :param typemeta: Data type of the column
        :param metafile: TSV file with the format (Genome_id \t Value)
        '''
        try:
            data_list = []
            with open(metafile, 'r') as metaf:
                for line in metaf:
                    data_list.append(tuple(line.strip().split('\t')))
            data_zip = zip(*data_list)
            genome_id = list(data_zip[0])
            meta_value = list(data_zip[1])
            for n, i in enumerate(genome_id):
                new_i = i.split("_", 1)[1]
                genome_id[n] = new_i
            query = "SELECT upsert('{0}','{1}','{2}',%s,%s)".format(
                table, field, typemeta)
            self.cur.execute(query, (genome_id, meta_value))
        except GenomeDatabaseError as e:
            raise e
        except psycopg2.Error as e:
            raise GenomeDatabaseError(e.pgerror)

    def createMetadata(self, metadatafile):
        '''
        Function createMetadata
        Create or Update metaddata columns in the database

        :param metadatafile: TSV file listing one new field per line
        Format of the TSV file is new_field \t description \t type \t table
        '''
        try:
            data_dict = {}
            with open(metadatafile, 'r') as metaf:
                for line in metaf:
                    array_line = line.strip().split('\t')
                    if not array_line[3].startswith("metadata_"):
                        raise GenomeDatabaseError(
                            "Only Metadata Tables can be modified")
                    data_dict[array_line[0]] = {
                        "table": array_line[3], "type": array_line[2], "desc": array_line[1]}

            query = "SELECT v.field,v.table from view_list_meta_columns as v"
            self.cur.execute(query)
            all_col_dict = dict(self.cur.fetchall())
            for key, value in data_dict.iteritems():
                if key in all_col_dict:
                    if all_col_dict.get(key) == value['table']:
                        query_comment = "COMMENT ON COLUMN {0}.{1} IS '{2}'".format(
                            value['table'], key, value['desc'])
                        self.cur.execute(query_comment)
                    else:
                        logging.warning("Column {0} is already presents in the {1} table .".format(
                            key, all_col_dict.get(key)))
                else:
                    query_add_col = "ALTER TABLE {0} ADD COLUMN {1} {2}".format(
                        value['table'], key, value['type'])
                    self.cur.execute(query_add_col)
                    query_add_comment = "COMMENT ON COLUMN {0}.{1} IS '{2}'".format(
                        value['table'], key, value['desc'])
                    self.cur.execute(query_add_comment)

# ---------- PSQL are not refresh automatically so we need to drop the existing view and recreate it with a new Definition.
            self.cur.execute("SELECT refreshView()")

        except GenomeDatabaseError as e:
            raise e

    def addMetadata(self, db_genome_id, genome_file, gff_file, checkm_results, output_dir):
        """Calculate and add metadata to DB.

        Parameters
        ----------
        db_genome_id : str
            Unique database identifer of genome.
        genome_file : str
            Name of FASTA file containing nucleotide sequences.
        gff_file : str
            Name of generic feature file describing genes.
        checkm_results : dict
            CheckM metadata.
        output_dir : str
            Output directory.
        """

        # create rows for genome in metadata tables
        self.cur.execute(
            "INSERT INTO metadata_nucleotide (id) VALUES ({0})".format(db_genome_id))
        self.cur.execute(
            "INSERT INTO metadata_genes (id) VALUES ({0})".format(db_genome_id))
        self.cur.execute(
            "INSERT INTO metadata_taxonomy (id) VALUES ({0})".format(db_genome_id))
        self.cur.execute(
            "INSERT INTO metadata_rna (id) VALUES ({0})".format(db_genome_id))
        self.cur.execute(
            "INSERT INTO metadata_rrna_sequences (id) VALUES ({0})".format(db_genome_id))

        self._calculateMetadata(genome_file, gff_file, output_dir)
        self._storeMetadata(db_genome_id, output_dir)
        self._storeCheckM(db_genome_id, checkm_results)

        return True

    def _calculateMetadata(self, genome_file, gff_file, output_dir):
        """Calculate metadata for new genome.

        Parameters
        ----------
        genome_file : str
            Name of FASTA file containing nucleotide sequences.
        gff_file : str
            Name of generic feature file describing genes.
        output_dir : str
            Output directory.
        """
        os.system('genometk nucleotide --silent %s %s' %
                  (genome_file, output_dir))
        os.system('genometk gene --silent %s %s %s' %
                  (genome_file, gff_file, output_dir))
        os.system('genometk rna --silent --db %s --taxonomy_file %s %s ssu %s' % (
            ConfigMetadata.GTDB_SSU_GG_DB,
            ConfigMetadata.GTDB_SSU_GG_TAXONOMY,
            genome_file,
            os.path.join(output_dir, ConfigMetadata.GTDB_SSU_GG_OUTPUT_DIR)))
        os.system('genometk rna --silent --db %s --taxonomy_file %s %s ssu %s' % (
            ConfigMetadata.GTDB_SSU_SILVA_DB,
            ConfigMetadata.GTDB_SSU_SILVA_TAXONOMY,
            genome_file,
            os.path.join(output_dir, ConfigMetadata.GTDB_SSU_SILVA_OUTPUT_DIR)))
        os.system('genometk rna --silent --db %s --taxonomy_file %s %s lsu_23S %s' % (
            ConfigMetadata.GTDB_LSU_SILVA_DB,
            ConfigMetadata.GTDB_LSU_SILVA_TAXONOMY,
            genome_file,
            os.path.join(output_dir, ConfigMetadata.GTDB_LSU_SILVA_OUTPUT_DIR)))
        os.system('genometk rna --silent --cpus 1 --min_len 80 %s lsu_5S %s' % (
            genome_file,
            os.path.join(output_dir, ConfigMetadata.GTDB_LSU_5S_OUTPUT_DIR)))
        return True

    def _storeMetadata(self, db_genome_id, genome_dir):
        """Parse metadata files for genome and store in database.

        Parameters
        ----------
        db_genome_id : str
            Unique database identifier of genome.
        genome_dir : str
            Directory containing metadata files to parse.
        """
        try:
            # nucleotide metadata
            metadata_nt_path = os.path.join(
                genome_dir, ConfigMetadata.GTDB_NT_FILE)
            genome_list_nt = [tuple(line.rstrip().split('\t'))
                              for line in open(metadata_nt_path)]
            query_nt = "UPDATE metadata_nucleotide SET %s = %s WHERE id = {0}".format(
                db_genome_id)
            for c, v in genome_list_nt:
                try:
                    v = float(v)
                    self.cur.execute(query_nt, [AsIs(c), v])
                except:
                    self.cur.execute(query_nt, [AsIs(c), v])

            # protein metadata
            metadata_gene_path = os.path.join(
                genome_dir, ConfigMetadata.GTDB_GENE_FILE)
            genome_list_gene = [tuple(line.rstrip().split('\t'))
                                for line in open(metadata_gene_path)]
            query_gene = "UPDATE metadata_genes SET %s = %s WHERE id = {0}".format(
                db_genome_id)
            for c, v in genome_list_gene:
                try:
                    v = float(v)
                    self.cur.execute(query_gene, [AsIs(c), v])
                except:
                    self.cur.execute(query_gene, [AsIs(c), v])

            # Greengenes SSU metadata
            query_taxonomy = "UPDATE metadata_rna SET %s = %s WHERE id = {0}".format(
                db_genome_id)
            query_sequence = "UPDATE metadata_rrna_sequences SET %s = %s WHERE id = {0}".format(
                db_genome_id)
            metadata_ssu_gg_path = os.path.join(
                genome_dir, ConfigMetadata.GTDB_SSU_GG_OUTPUT_DIR, ConfigMetadata.GTDB_SSU_FILE)
            genome_list_taxonomy, _ssu_count, ssu_query_id = self._parse_taxonomy_file(
                metadata_ssu_gg_path, ConfigMetadata.GTDB_SSU_GG_PREFIX)
            if genome_list_taxonomy:
                for c, v in genome_list_taxonomy:
                    try:
                        if "blast_subject_id" not in c:
                            v = float(v)
                        self.cur.execute(query_taxonomy, [AsIs(c), v])
                    except:
                        self.cur.execute(query_taxonomy, [AsIs(c), v])

            # SILVA SSU metadata saved in metadata_ssu table [HACK: eventually
            # information will only be stored in this table]
            query_taxonomy = "UPDATE metadata_rna SET %s = %s WHERE id = {0}".format(
                db_genome_id)
            query_sequence = "UPDATE metadata_rrna_sequences SET %s = %s WHERE id = {0}".format(
                db_genome_id)
            metadata_ssu_silva_path = os.path.join(
                genome_dir, ConfigMetadata.GTDB_SSU_SILVA_OUTPUT_DIR, ConfigMetadata.GTDB_SSU_FILE)
            metadata_ssu_fna_silva_path = os.path.join(
                genome_dir, ConfigMetadata.GTDB_SSU_SILVA_OUTPUT_DIR, ConfigMetadata.GTDB_SSU_FNA_FILE)
            metadata_ssu_silva_summary_file = os.path.join(
                genome_dir, ConfigMetadata.GTDB_SSU_SILVA_OUTPUT_DIR, ConfigMetadata.GTDB_SSU_SILVA_SUMMARY_FILE)

            genome_list_taxonomy, ssu_count, ssu_query_id = self._parse_taxonomy_file(
                metadata_ssu_silva_path, ConfigMetadata.GTDB_SSU_SILVA_PREFIX, metadata_ssu_silva_summary_file)
            if genome_list_taxonomy:
                for c, v in genome_list_taxonomy:
                    try:
                        if "blast_subject_id" not in c:
                            v = float(v)
                        self.cur.execute(query_taxonomy, [AsIs(c), v])
                    except:
                        self.cur.execute(query_taxonomy, [AsIs(c), v])
                if ssu_query_id is not None:
                    genome_list_sequence = self._parse_sequence_file(
                        metadata_ssu_fna_silva_path, ConfigMetadata.GTDB_SSU_SILVA_PREFIX, ssu_query_id)
                    for c, v in genome_list_sequence:
                        self.cur.execute(query_sequence, [AsIs(c), v])

            # SILVA LSU metadata saved in metadata_ssu table [HACK: eventually
            # information will only be stored in this table]
            query_taxonomy = "UPDATE metadata_rna SET %s = %s WHERE id = {0}".format(
                db_genome_id)
            metadata_lsu_silva_path = os.path.join(
                genome_dir, ConfigMetadata.GTDB_LSU_SILVA_OUTPUT_DIR, ConfigMetadata.GTDB_LSU_FILE)
            metadata_lsu_fna_silva_path = os.path.join(
                genome_dir, ConfigMetadata.GTDB_LSU_SILVA_OUTPUT_DIR, ConfigMetadata.GTDB_LSU_FNA_FILE)
            metadata_lsu_silva_summary_file = os.path.join(
                genome_dir, ConfigMetadata.GTDB_LSU_SILVA_OUTPUT_DIR, ConfigMetadata.GTDB_LSU_SILVA_SUMMARY_FILE)
            genome_list_taxonomy, lsu_count, lsu_query_id = self._parse_taxonomy_file(
                metadata_lsu_silva_path, ConfigMetadata.GTDB_LSU_SILVA_PREFIX, metadata_lsu_silva_summary_file)
            if genome_list_taxonomy:
                for c, v in genome_list_taxonomy:
                    try:
                        if "blast_subject_id" not in c:
                            v = float(v)
                        self.cur.execute(query_taxonomy, [AsIs(c), v])
                    except:
                        self.cur.execute(query_taxonomy, [AsIs(c), v])
                if lsu_query_id is not None:
                    genome_list_sequence = self._parse_sequence_file(
                        metadata_lsu_fna_silva_path, ConfigMetadata.GTDB_LSU_SILVA_PREFIX, lsu_query_id)
                    for c, v in genome_list_sequence:
                        self.cur.execute(query_sequence, [AsIs(c), v])

            # SILVA 5s LSU metadata saved in metadata_ssu table [HACK:
            # eventually information will only be stored in this table]
            query_taxonomy = "UPDATE metadata_rna SET %s = %s WHERE id = {0}".format(
                db_genome_id)
            metadata_lsu_5s_path = os.path.join(
                genome_dir, ConfigMetadata.GTDB_LSU_5S_OUTPUT_DIR, ConfigMetadata.GTDB_LSU_5S_FILE)
            metadata_lsu_5s_fna_path = os.path.join(
                genome_dir, ConfigMetadata.GTDB_LSU_5S_OUTPUT_DIR, ConfigMetadata.GTDB_LSU_5S_FNA_FILE)
            metadata_lsu_5s_summary_file = os.path.join(
                genome_dir, ConfigMetadata.GTDB_LSU_5S_OUTPUT_DIR, ConfigMetadata.GTDB_LSU_5S_SUMMARY_FILE)
            genome_list_taxonomy, lsu_5s_count, lsu_5s_query_id = self._parse_taxonomy_file(
                metadata_lsu_5s_path, ConfigMetadata.GTDB_LSU_5S_PREFIX, metadata_lsu_5s_summary_file)
            if genome_list_taxonomy:
                for c, v in genome_list_taxonomy:
                    try:
                        if "blast_subject_id" not in c:
                            v = float(v)
                        self.cur.execute(query_taxonomy, [AsIs(c), v])
                    except:
                        self.cur.execute(query_taxonomy, [AsIs(c), v])
                if lsu_5s_query_id is not None:
                    genome_list_sequence = self._parse_sequence_file(
                        metadata_lsu_5s_fna_path, ConfigMetadata.GTDB_LSU_5S_PREFIX, lsu_5s_query_id)
                    for c, v in genome_list_sequence:
                        self.cur.execute(query_sequence, [AsIs(c), v])

            query_gene_ssu = "UPDATE metadata_genes SET ssu_count = %s WHERE id = {0}".format(
                db_genome_id)
            self.cur.execute(query_gene_ssu, (ssu_count,))

            query_gene_lsu = "UPDATE metadata_genes SET lsu_23s_count = %s WHERE id = {0}".format(
                db_genome_id)
            self.cur.execute(query_gene_lsu, (lsu_count,))

            query_gene_lsu = "UPDATE metadata_genes SET lsu_5s_count = %s WHERE id = {0}".format(
                db_genome_id)
            self.cur.execute(query_gene_lsu, (lsu_5s_count,))

            return True
        except psycopg2.Error as e:
            print "error"
            raise GenomeDatabaseError(e.pgerror)
        except:
            print("Unexpected error:", sys.exc_info()[0])
            raise

    def _parse_taxonomy_file(self, metadata_taxonomy_file, prefix, summary_file=None):
        """Parse metadata file with taxonomic information for rRNA genes.

        Parameters
        ----------
        metadata_taxonomy_file : str
          Full path to file containing 16S rRNA metadata.
        Prefix : str
          Prefix to append to metadata fields.
        Fna file : str
          Full path to file containing 16S rRNA sequences.

        Returns
        -------
        list of tuples
          Field, value pairs for all metadata items.
        int
            Number of 16S sequences.
        """

        if not os.path.exists(metadata_taxonomy_file):
            return None, 0, None

        metadata = []
        ssu_query_id = None

        num_lines = sum(1 for line in open(metadata_taxonomy_file))
        if num_lines == 1:  # Empty file
            return metadata, 0, None

        with open(metadata_taxonomy_file) as f:
            header_line = f.readline().rstrip()
            headers = [prefix + '_' + x for x in header_line.split('\t')]

            if prefix == 'lsu_silva_23s':
                for n, i in enumerate(headers):
                    if i == 'lsu_silva_23s_sequence':
                        headers[n] = 'lsu_23s_sequence'
                    elif i == 'lsu_silva_23s_query_id':
                        headers[n] = 'lsu_23s_query_id'
                    elif i == 'lsu_silva_23s_length':
                        headers[n] = 'lsu_23s_length'
            elif prefix == 'ssu_silva':
                for n, i in enumerate(headers):
                    if i == 'ssu_silva_sequence':
                        headers[n] = 'ssu_sequence'
                    elif i == 'ssu_silva_query_id':
                        headers[n] = 'ssu_query_id'
                    elif i == 'ssu_silva_length':
                        headers[n] = 'ssu_length'

            split_headers = header_line.rstrip().split("\t")
            for pos in range(0, len(split_headers)):
                header = split_headers[pos]
                if header == 'query_id':
                    query_id_pos = pos
                    break

            # Report hit to longest 16S rRNA gene. It is possible that
            # the HMMs identified a putative 16S rRNA gene, but that
            # there was no valid BLAST hit.
            longest_query_len = 0
            longest_ssu_hit_info = None
            identified_ssu_genes = 0
            for line in f:
                line_split = line.strip().split('\t')
                query_len = int(line_split[2])
                if query_len > longest_query_len:
                    longest_query_len = query_len
                    longest_ssu_hit_info = line_split
                    ssu_query_id = line_split[query_id_pos]

            if longest_ssu_hit_info:
                metadata = [(headers[i], value)
                            for i, value in enumerate(longest_ssu_hit_info)]

                if summary_file is not None and os.path.exists(summary_file):
                    with open(summary_file) as fsum:
                        header_line = fsum.readline()  # consume header line
                        header_list = [x.strip()
                                       for x in header_line.split('\t')]
                        idx_seq = header_list.index("Sequence length")
                        for line in fsum:
                            identified_ssu_genes += 1
                            sum_list = [x.strip() for x in line.split('\t')]
                            if sum_list[0] == ssu_query_id:
                                if prefix == 'lsu_silva_23s':
                                    metadata.append(
                                        ('lsu_23s_contig_len', sum_list[idx_seq]))
                                elif prefix == 'ssu_silva':
                                    metadata.append(
                                        ('ssu_contig_len', sum_list[idx_seq]))
                                else:
                                    metadata.append(
                                        ('{0}_contig_len'.format(prefix), sum_list[idx_seq]))

            return metadata, identified_ssu_genes, ssu_query_id

    def _parse_sequence_file(self, fna_file, prefix, ssu_query_id):
        metadata = []
        all_genes_dict = read_fasta(fna_file, False)
        sequence = all_genes_dict[ssu_query_id]
        if prefix == 'lsu_silva_23s':
            metadata.append(('lsu_23s_sequence'.format(prefix), sequence))
        elif prefix == 'ssu_silva':
            metadata.append(('ssu_sequence'.format(prefix), sequence))
        else:
            metadata.append(('{0}_sequence'.format(prefix), sequence))
        return metadata

    def _storeCheckM(self, db_genome_id, checkm_results):
        """Store CheckM results for genome.

        Parameters
        ----------
        db_genome_id : str
            Unique database identifer of genome.
        checkm_results : dict
            CheckM metadata.
        """
        try:
            checkm_data = [('checkm_completeness', checkm_results['completeness']),
                           ('checkm_contamination',
                            checkm_results['contamination']),
                           ('checkm_marker_count',
                            checkm_results['marker_count']),
                           ('checkm_marker_lineage',
                            checkm_results['lineage']),
                           ('checkm_genome_count',
                            checkm_results['genome_count']),
                           ('checkm_marker_set_count',
                            checkm_results['set_count']),
                           ('checkm_strain_heterogeneity', checkm_results['heterogeneity'])]

            query = "UPDATE metadata_genes SET %s = %s WHERE id = {0}".format(
                db_genome_id)
            self.cur.executemany(query, [(AsIs(c), v)
                                         for (c, v) in checkm_data])
            return True
        except:
            print("Unexpected error:", sys.exc_info()[0])
            raise

        def run_trna_scan(genome_path, accession, domain):
            user_path = os.path.join(Config.GTDB_GENOME_ROOT, Config.GTDB_GENOME_ROOT, Config.USER_PREFIX, os.path.dirname(
                genome_path), ConfigMetadata.GTDB_TRNA_OUTPUT_DIR)
            genome_file = os.path.join(
                Config.GTDB_GENOME_ROOT, Config.GTDB_GENOME_ROOT, Config.USER_PREFIX, genome_path)
            output_file = os.path.join(user_path, accession + '_trna.tsv')
            log_file = os.path.join(user_path, accession + '_trna.log')
            stats_file = os.path.join(user_path, accession + '_trna_stats.tsv')

            domain_flag = '-B'
            if domain == 'd__Archaea':
                domain_flag = '-A'

            cmd = 'tRNAscan-SE %s -q -Q -o %s -m %s -l %s %s' % (domain_flag,
                                                                 output_file,
                                                                 stats_file,
                                                                 log_file,
                                                                 genome_file)

        def parse_trna_file(genome_id, trna_file):
            """Parse tRNA information."""

            if not os.path.exists(trna_file):
                return
            # parse tRNA summary file
            trna_count = 0
            trna_selenocysteine_count = 0
            trna_aa_count = 0
            read_aa = False
            for line in open(trna_file):
                if line.startswith('tRNAs decoding Standard 20 AA'):
                    trna_count = int(line.split(':')[1])
                elif line.startswith('Selenocysteine tRNAs (TCA)'):
                    trna_selenocysteine_count = int(line.split(':')[1])
                    trna_count += trna_selenocysteine_count
                elif line.startswith('Isotype / Anticodon Counts:'):
                    read_aa = True

                if read_aa:
                    # Parsing lines with the following format:
                    # Ala   : 2      AGC:         GGC:         CGC: 1
                    # TGC: 1
                    if ' : ' in line:
                        line_split = line.split('\t')[0]
                        line_split = map(str.strip, line_split.split(':'))
                        if len(line_split[0]) == 3:  # this is an amino acid
                            if int(line_split[1]) > 0:
                                trna_aa_count += 1

        def store_trna_file(genome_id, trna_count, trna_selenocysteine_count, trna_aa_count):
            query_trna = "UPDATE metadata_nucleotide SET trna_count = {0} , trna_selenocysteine_count ={1} , trna_aa_count = {2} WHERE id = {3}".format(trna_count,
                                                                                                                                                        trna_selenocysteine_count,
                                                                                                                                                        trna_aa_count,
                                                                                                                                                        db_genome_id)
            self.cur.execute(query_trna)

    def _parse_silva_mapping(self, path):
        # type: (str) -> dict
        """ Parses the SILVA taxonomy to rank file and returns it as a dictionary.
        e.g. https://www.arb-silva.de/fileadmin/silva_databases/release_123_1/Exports/taxonomy/tax_slv_ssu_123.1.txt

        Parameters
        ----------
        :param path: str
            The path to the file provided.
        :return: dict
            A dictionary containing the mapping.
        """
        out = dict()
        with open(path, 'r') as f:
            re_hits = re.findall(r'(.+;)\t+[\d]+\t+(\w+)\t.*\n', f.read())
        for hit in re_hits:
            out[hit[0]] = hit[1]  # out[full SILVA taxonomy string] = [rank]
        return out

    def _parse_silva_ssu_aggregate(self, path):
        # type: (str) -> dict
        """ Parses the SSU file which is an aggregate of all BLAST hits against 16s seqs for representatives.

        Parameters
        ----------
        :param path str
            The path to the aggregate file
        :return: dict
            A dictionary containing the BLAST hits for each genome id.
        """
        out = defaultdict(list)
        with open(path, 'r') as f:
            for line in f.readlines():
                line = line.rstrip()
                vals = line.split('\t')
                out[int(vals[0])].append({'query_id': vals[1], 'taxonomy': vals[2], 'length': int(vals[3]),
                                          'blast_subject_id': vals[4], 'blast_evalue': float(vals[5]),
                                          'blast_bitscore': float(vals[6]), 'blast_align_len': int(vals[7]),
                                          'blast_perc_identity': float(vals[8])})
        return out
