"""Provide a (semi) automated mapping between GenBank and ontologies.

The goal is to provide an ontology namespace and term for each item
of a standard GenBank file.
"""
from __future__ import with_statement
import os
import collections
import csv

from rdflib.Graph import Graph
#from Mgh.SemanticLiterature import sparta
import sparta

def main(ft_file, so_ft_map_file, so_file, dc_file, rdfs_file):
    # -- get all of the keys we need to map from GenBank
    # pulled by hand from the Biopython parser
    header_keys = ['ACCESSION', 'AUTHORS', 'COMMENT', 'CONSRTM', 'DBLINK',
            'DBSOURCE', 'DEFINITION', 'JOURNAL', 'KEYWORDS', 'MEDLINE', 'NID',
            'ORGANISM', 'PID', 'PROJECT', 'PUBMED',
            'REMARK', 'SEGMENT', 'SOURCE', 'TITLE', 'VERSION',
            'VERSION;gi', 'LOCUS;date']
    feature_keys, qual_keys = parse_feature_table(ft_file)

    # -- terms to map them two, along with dictionaries for the existing or
    # hand defined mappings
    dc_terms = parse_dc_terms(dc_file)
    dc_ontology = OntologyGroup("http://purl.org/dc/terms/", dc_terms)
    so_terms = parse_so_terms(so_file)
    so_ontology = OntologyGroup("http://purl.org/obo/owl/SO", so_terms)
    so_ft_map = parse_so_ft_map(so_ft_map_file)
    so_ontology.add_map('feature', 
            'http://www.sequenceontology.org/mappings/FT_SO_map.txt', so_ft_map)
    rdfs_terms = parse_rdfs_terms(rdfs_file)
    rdfs_ontology = OntologyGroup('http://www.w3.org/2000/01/rdf-schema#',
            rdfs_terms)
    me_ft_map = dict(
            rep_origin = 'origin_of_replication',
            unsure = 'sequence_uncertainty',
            conflict = 'sequence_conflict',
            GC_signal = 'GC_rich_promoter_region',
            mat_peptide = 'mature_protein_region',
            C_region = 'C_cluster',
            J_segment = 'J_gene',
            N_region = '',
            S_region = '',
            V_region = 'V_cluster',
            V_segment = 'V_gene'
            )
    me_ft_dc_map = dict(
            old_sequence = 'replaces')
    me_hd_map = {'ACCESSION': 'databank_entry'}
    me_hd_dc_map = {'DEFINITION' : 'description',
                    'VERSION' : 'hasVersion',
                    'VERSION;gi' : 'identifier',
                    'KEYWORDS' : 'subject',
                    'LOCUS;date' : 'created',
                    'PUBMED' : 'relation',
                    'JOURNAL' : 'source',
                    'AUTHORS' : 'contributor',
                    'CONSRTM' : 'creator',
                    'ORGANISM' : '',
                    'SEGMENT' : 'isPartOf',
                    'DBLINK' : 'relation',
                    'DBSOURCE' : 'relation',
                    'MEDLINE' : 'relation',
                    'NID' : 'relation',
                    'PID' : 'relation',
                    'PROJECT' : 'relation'}
    me_ql_map = dict(
            bio_material = 'biomaterial_region',
            bound_moiety = 'bound_by_factor',
            codon_start = 'coding_start',
            direction = 'direction_attribute',
            experiment = 'experimental_result_region',
            macronuclear = 'macronuclear_sequence',
            map = 'fragment_assembly',
            mobile_element = 'integrated_mobile_genetic_element',
            mod_base = 'modified_base_site',
            mol_type = 'sequence_attribute',
            ncRNA_class = 'ncRNA',
            organelle = 'organelle_sequence',
            PCR_primers = 'primer',
            proviral = 'proviral_region',
            pseudo = 'pseudogene',
            rearranged = 'rearranged_at_DNA_level',
            satellite = 'satellite_DNA',
            segment = 'gene_segment',
            rpt_family = 'repeat_family',
            rpt_type = 'repeat_unit',
            rpt_unit_range = 'repeat_region',
            rpt_unit_seq = 'repeat_component',
            tag_peptide = 'cleaved_peptide_region',
            trans_splicing = 'trans_spliced',
            translation = 'polypeptide',
            )
    me_ql_dc_map = dict(
            citation = 'bibliographicCitation',
            gene_synonym = 'alternative',
            identified_by = 'creator',
            label = 'alternative',
            locus_tag = 'alternative',
            old_locus_tag = 'replaces',
            replace = 'isReplacedBy',
            db_xref = 'relation',
            compare = 'relation',
            EC_number = 'relation',
            protein_id = 'relation',
            product = 'alternative',
            standard_name = 'alternative',
            number = 'coverage',
            function = 'description',
            )
    me_ql_rdfs_map = dict(
            note = 'comment',
            )
    ql_make_no_sense_list = ['number']
    so_ontology.add_map('header', 'Brad', me_hd_map)
    so_ontology.add_map('feature', 'Brad', me_ft_map)
    so_ontology.add_map('qualifier', 'Brad', me_ql_map)
    dc_ontology.add_map('header', 'Brad', me_hd_dc_map)
    dc_ontology.add_map('feature', 'Brad', me_ft_dc_map)
    dc_ontology.add_map('qualifier', 'Brad', me_ql_dc_map)
    rdfs_ontology.add_map('qualifier', 'Brad', me_ql_rdfs_map)

    # -- write out the mappings in each of the categories
    with open("genbank_ontology_map.txt", "w") as out_handle:
        out_writer = csv.writer(out_handle, delimiter="\t")
        out_writer.writerow(['gb section', 'identifier', 'ontology',
            'namespace', 'evidence'])
        match_keys_to_ontology('header', header_keys, [so_ontology, dc_ontology,
            rdfs_ontology], out_writer)
        match_keys_to_ontology('feature', feature_keys, [so_ontology,
            dc_ontology, rdfs_ontology], out_writer)
        match_keys_to_ontology('qualifier', qual_keys, [so_ontology,
            dc_ontology, rdfs_ontology], out_writer)

class OntologyGroup:
    def __init__(self, namespace, terms):
        self.ns = namespace
        self.terms = terms
        self._maps = collections.defaultdict(lambda: [])

    def add_map(self, key_type, origin, key_map):
        """Add a mapping of keys to terms within this ontology.
        """
        self._maps[key_type].append((origin, key_map))

    def normalized_terms(self):
        """Retrieve the terms all lower cased and with extra items removed.
        """
        lower_so_terms = {}
        for term in self.terms:
            lower_so_terms[self._normal_term(term)] = term
        return lower_so_terms

    def _normal_term(self, term):
        return term.replace("_", "").lower() 

    def match_key_to_ontology(self, key_type, cur_key):
        normal_terms = self.normalized_terms()
        # try to get it from a dictionary
        for map_origin, check_map in self._maps[key_type]:
            try:
                match_key = check_map[cur_key]
            except KeyError:
                match_key = None
            if (match_key and
                    self._normal_term(match_key) in normal_terms.keys()):
                return match_key, map_origin
        # try to match it by name
        if self._normal_term(cur_key) in normal_terms.keys():
            match_key = normal_terms[self._normal_term(cur_key)]
            return match_key, 'name_match'
        #in_terms = [x for x in normal_terms.keys() if
        #        x.find(self._normal_term(cur_key)) != -1]
        #if len(in_terms) > 0:
        #    print '***', cur_key, in_terms, self.ns
        # could not find a match
        return None, ''

def match_keys_to_ontology(key_type, keys, ontologies, out_writer):
    no_matches = []
    for cur_key in keys:
        found_match = False
        for ontology in ontologies:
            match_key, origin = ontology.match_key_to_ontology(key_type, cur_key)
            if match_key is not None:
                out_writer.writerow([key_type, cur_key, match_key, ontology.ns,
                    origin])
                found_match = True
                break
        if not found_match:
            no_matches.append(cur_key)
    for no_key in no_matches:
        out_writer.writerow([key_type, no_key])

def _parse_terms_from_rdf(in_file, prefix):
    graph = Graph()
    graph.parse(in_file)
    sparta_store = sparta.ThingFactory(graph)
    subjects = []
    for subj in graph.subjects():
        if subj not in subjects:
            subjects.append(subj)
    terms = []
    for subj in subjects:
        rdf_item = sparta_store(subj)
        str_item = str(rdf_item.get_id().replace(prefix + "_", ""))
        if str_item:
            terms.append(str_item)
    terms.sort()
    #print terms
    return terms

def parse_rdfs_terms(rdfs_file):
    return _parse_terms_from_rdf(rdfs_file, "rdfs")

def parse_dc_terms(dc_file):
    """Retrieve a list of Dublin core terms from the RDF file.
    """
    return _parse_terms_from_rdf(dc_file, "dcterms")

def parse_so_terms(so_file):
    """Retrieve all available Sequence Ontology terms from the file.
    """
    so_terms = []
    with open(so_file) as in_handle:
        for line in in_handle:
            if line.find('name:') == 0:
                name = line[5:].strip()
                so_terms.append(name)
    return so_terms

def parse_so_ft_map(so_ft_map_file):
    """Parse out mappings between feature keys and SO.
    """
    so_ft_map = {}
    with open(so_ft_map_file) as in_handle:
        in_handle.readline()
        for line in in_handle:
            parts = line.split()
            if parts[1] not in ['undefined']:
                so_ft_map[parts[0]] = parts[1]
    return so_ft_map

def parse_feature_table(ft_file):
    """Parse all available features and qualifiers from the FT definition.

    This is ugly and parses it straight out of the HTML but this is much easier
    than trying to get it from the specs.
    """
    feature_keys = []
    qual_keys = []
    with open(ft_file) as ft_handle:
        in_feature_region = False
        for line in ft_handle:
            if in_feature_region:
                if line.strip() == "":
                    in_feature_region = False
                else:
                    qual_key, feature_key = line.strip().split()
                    qual_keys.append(qual_key)
                    feature_keys.append(feature_key)
            elif line.find('QUALIFIER FEATURE KEY') == 0:
                in_feature_region = True
    qual_keys = list(set(qual_keys))
    qual_keys = [k.replace('/', '') for k in qual_keys]
    feature_keys = list(set(feature_keys))
    qual_keys.sort()
    feature_keys.sort()
    return feature_keys, qual_keys

if __name__ == "__main__":
    ft_file = os.path.join("feature_table", "FT_index.html")
    so_ft_map_file = os.path.join("feature_table", "FT_SO_map.txt")
    so_file = os.path.join("ontologies", "so.obo")
    dc_file = os.path.join("ontologies", "dcterms.rdf")
    rdfs_file = os.path.join("ontologies", "rdf-schema.rdf")
    main(ft_file, so_ft_map_file, so_file, dc_file, rdfs_file)

"""Provide a (semi) automated mapping between GenBank and ontologies.

The goal is to provide an ontology namespace and term for each item
of a standard GenBank file.
"""
from __future__ import with_statement
import os
import collections
import csv

from rdflib.Graph import Graph
#from Mgh.SemanticLiterature import sparta
import sparta

def main(ft_file, so_ft_map_file, so_file, dc_file, rdfs_file):
    # -- get all of the keys we need to map from GenBank
    # pulled by hand from the Biopython parser
    header_keys = ['ACCESSION', 'AUTHORS', 'COMMENT', 'CONSRTM', 'DBLINK',
            'DBSOURCE', 'DEFINITION', 'JOURNAL', 'KEYWORDS', 'MEDLINE', 'NID',
            'ORGANISM', 'PID', 'PROJECT', 'PUBMED',
            'REMARK', 'SEGMENT', 'SOURCE', 'TITLE', 'VERSION',
            'VERSION;gi', 'LOCUS;date']
    feature_keys, qual_keys = parse_feature_table(ft_file)

    # -- terms to map them two, along with dictionaries for the existing or
    # hand defined mappings
    dc_terms = parse_dc_terms(dc_file)
    dc_ontology = OntologyGroup("http://purl.org/dc/terms/", dc_terms)
    so_terms = parse_so_terms(so_file)
    so_ontology = OntologyGroup("http://purl.org/obo/owl/SO", so_terms)
    so_ft_map = parse_so_ft_map(so_ft_map_file)
    so_ontology.add_map('feature', 
            'http://www.sequenceontology.org/mappings/FT_SO_map.txt', so_ft_map)
    rdfs_terms = parse_rdfs_terms(rdfs_file)
    rdfs_ontology = OntologyGroup('http://www.w3.org/2000/01/rdf-schema#',
            rdfs_terms)
    me_ft_map = dict(
            rep_origin = 'origin_of_replication',
            unsure = 'sequence_uncertainty',
            conflict = 'sequence_conflict',
            GC_signal = 'GC_rich_promoter_region',
            mat_peptide = 'mature_protein_region',
            C_region = 'C_cluster',
            J_segment = 'J_gene',
            N_region = '',
            S_region = '',
            V_region = 'V_cluster',
            V_segment = 'V_gene'
            )
    me_ft_dc_map = dict(
            old_sequence = 'replaces')
    me_hd_map = {'ACCESSION': 'databank_entry'}
    me_hd_dc_map = {'DEFINITION' : 'description',
                    'VERSION' : 'hasVersion',
                    'VERSION;gi' : 'identifier',
                    'KEYWORDS' : 'subject',
                    'LOCUS;date' : 'created',
                    'PUBMED' : 'relation',
                    'JOURNAL' : 'source',
                    'AUTHORS' : 'contributor',
                    'CONSRTM' : 'creator',
                    'ORGANISM' : '',
                    'SEGMENT' : 'isPartOf',
                    'DBLINK' : 'relation',
                    'DBSOURCE' : 'relation',
                    'MEDLINE' : 'relation',
                    'NID' : 'relation',
                    'PID' : 'relation',
                    'PROJECT' : 'relation'}
    me_ql_map = dict(
            bio_material = 'biomaterial_region',
            bound_moiety = 'bound_by_factor',
            codon_start = 'coding_start',
            direction = 'direction_attribute',
            experiment = 'experimental_result_region',
            macronuclear = 'macronuclear_sequence',
            map = 'fragment_assembly',
            mobile_element = 'integrated_mobile_genetic_element',
            mod_base = 'modified_base_site',
            mol_type = 'sequence_attribute',
            ncRNA_class = 'ncRNA',
            organelle = 'organelle_sequence',
            PCR_primers = 'primer',
            proviral = 'proviral_region',
            pseudo = 'pseudogene',
            rearranged = 'rearranged_at_DNA_level',
            satellite = 'satellite_DNA',
            segment = 'gene_segment',
            rpt_family = 'repeat_family',
            rpt_type = 'repeat_unit',
            rpt_unit_range = 'repeat_region',
            rpt_unit_seq = 'repeat_component',
            tag_peptide = 'cleaved_peptide_region',
            trans_splicing = 'trans_spliced',
            translation = 'polypeptide',
            )
    me_ql_dc_map = dict(
            citation = 'bibliographicCitation',
            gene_synonym = 'alternative',
            identified_by = 'creator',
            label = 'alternative',
            locus_tag = 'alternative',
            old_locus_tag = 'replaces',
            replace = 'isReplacedBy',
            db_xref = 'relation',
            compare = 'relation',
            EC_number = 'relation',
            protein_id = 'relation',
            product = 'alternative',
            standard_name = 'alternative',
            number = 'coverage',
            function = 'description',
            )
    me_ql_rdfs_map = dict(
            note = 'comment',
            )
    ql_make_no_sense_list = ['number']
    so_ontology.add_map('header', 'Brad', me_hd_map)
    so_ontology.add_map('feature', 'Brad', me_ft_map)
    so_ontology.add_map('qualifier', 'Brad', me_ql_map)
    dc_ontology.add_map('header', 'Brad', me_hd_dc_map)
    dc_ontology.add_map('feature', 'Brad', me_ft_dc_map)
    dc_ontology.add_map('qualifier', 'Brad', me_ql_dc_map)
    rdfs_ontology.add_map('qualifier', 'Brad', me_ql_rdfs_map)

    # -- write out the mappings in each of the categories
    with open("genbank_ontology_map.txt", "w") as out_handle:
        out_writer = csv.writer(out_handle, delimiter="\t")
        out_writer.writerow(['gb section', 'identifier', 'ontology',
            'namespace', 'evidence'])
        match_keys_to_ontology('header', header_keys, [so_ontology, dc_ontology,
            rdfs_ontology], out_writer)
        match_keys_to_ontology('feature', feature_keys, [so_ontology,
            dc_ontology, rdfs_ontology], out_writer)
        match_keys_to_ontology('qualifier', qual_keys, [so_ontology,
            dc_ontology, rdfs_ontology], out_writer)

class OntologyGroup:
    def __init__(self, namespace, terms):
        self.ns = namespace
        self.terms = terms
        self._maps = collections.defaultdict(lambda: [])

    def add_map(self, key_type, origin, key_map):
        """Add a mapping of keys to terms within this ontology.
        """
        self._maps[key_type].append((origin, key_map))

    def normalized_terms(self):
        """Retrieve the terms all lower cased and with extra items removed.
        """
        lower_so_terms = {}
        for term in self.terms:
            lower_so_terms[self._normal_term(term)] = term
        return lower_so_terms

    def _normal_term(self, term):
        return term.replace("_", "").lower() 

    def match_key_to_ontology(self, key_type, cur_key):
        normal_terms = self.normalized_terms()
        # try to get it from a dictionary
        for map_origin, check_map in self._maps[key_type]:
            try:
                match_key = check_map[cur_key]
            except KeyError:
                match_key = None
            if (match_key and
                    self._normal_term(match_key) in normal_terms.keys()):
                return match_key, map_origin
        # try to match it by name
        if self._normal_term(cur_key) in normal_terms.keys():
            match_key = normal_terms[self._normal_term(cur_key)]
            return match_key, 'name_match'
        #in_terms = [x for x in normal_terms.keys() if
        #        x.find(self._normal_term(cur_key)) != -1]
        #if len(in_terms) > 0:
        #    print '***', cur_key, in_terms, self.ns
        # could not find a match
        return None, ''

def match_keys_to_ontology(key_type, keys, ontologies, out_writer):
    no_matches = []
    for cur_key in keys:
        found_match = False
        for ontology in ontologies:
            match_key, origin = ontology.match_key_to_ontology(key_type, cur_key)
            if match_key is not None:
                out_writer.writerow([key_type, cur_key, match_key, ontology.ns,
                    origin])
                found_match = True
                break
        if not found_match:
            no_matches.append(cur_key)
    for no_key in no_matches:
        out_writer.writerow([key_type, no_key])

def _parse_terms_from_rdf(in_file, prefix):
    graph = Graph()
    graph.parse(in_file)
    sparta_store = sparta.ThingFactory(graph)
    subjects = []
    for subj in graph.subjects():
        if subj not in subjects:
            subjects.append(subj)
    terms = []
    for subj in subjects:
        rdf_item = sparta_store(subj)
        str_item = str(rdf_item.get_id().replace(prefix + "_", ""))
        if str_item:
            terms.append(str_item)
    terms.sort()
    #print terms
    return terms

def parse_rdfs_terms(rdfs_file):
    return _parse_terms_from_rdf(rdfs_file, "rdfs")

def parse_dc_terms(dc_file):
    """Retrieve a list of Dublin core terms from the RDF file.
    """
    return _parse_terms_from_rdf(dc_file, "dcterms")

def parse_so_terms(so_file):
    """Retrieve all available Sequence Ontology terms from the file.
    """
    so_terms = []
    with open(so_file) as in_handle:
        for line in in_handle:
            if line.find('name:') == 0:
                name = line[5:].strip()
                so_terms.append(name)
    return so_terms

def parse_so_ft_map(so_ft_map_file):
    """Parse out mappings between feature keys and SO.
    """
    so_ft_map = {}
    with open(so_ft_map_file) as in_handle:
        in_handle.readline()
        for line in in_handle:
            parts = line.split()
            if parts[1] not in ['undefined']:
                so_ft_map[parts[0]] = parts[1]
    return so_ft_map

def parse_feature_table(ft_file):
    """Parse all available features and qualifiers from the FT definition.

    This is ugly and parses it straight out of the HTML but this is much easier
    than trying to get it from the specs.
    """
    feature_keys = []
    qual_keys = []
    with open(ft_file) as ft_handle:
        in_feature_region = False
        for line in ft_handle:
            if in_feature_region:
                if line.strip() == "":
                    in_feature_region = False
                else:
                    qual_key, feature_key = line.strip().split()
                    qual_keys.append(qual_key)
                    feature_keys.append(feature_key)
            elif line.find('QUALIFIER FEATURE KEY') == 0:
                in_feature_region = True
    qual_keys = list(set(qual_keys))
    qual_keys = [k.replace('/', '') for k in qual_keys]
    feature_keys = list(set(feature_keys))
    qual_keys.sort()
    feature_keys.sort()
    return feature_keys, qual_keys

if __name__ == "__main__":
    ft_file = os.path.join("feature_table", "FT_index.html")
    so_ft_map_file = os.path.join("feature_table", "FT_SO_map.txt")
    so_file = os.path.join("ontologies", "so.obo")
    dc_file = os.path.join("ontologies", "dcterms.rdf")
    rdfs_file = os.path.join("ontologies", "rdf-schema.rdf")
    main(ft_file, so_ft_map_file, so_file, dc_file, rdfs_file)

Provide a (semi) automated mapping between GenBank and ontologies.  The goal is to provide an ontology namespace and term for each item of a standard GenBank file. #from Mgh.SemanticLiterature import sparta # -- get all of the keys we need to map from GenBank # pulled by hand from the Biopython parser # -- terms to map them two, along with dictionaries for the existing or # hand defined mappings #', # -- write out the mappings in each of the categories Add a mapping of keys to terms within this ontology. Retrieve the terms all lower cased and with extra items removed. # try to get it from a dictionary # try to match it by name #in_terms = [x for x in normal_terms.keys() if #        x.find(self._normal_term(cur_key)) != -1] #if len(in_terms) > 0: #    print '***', cur_key, in_terms, self.ns # could not find a match #print terms Retrieve a list of Dublin core terms from the RDF file. Retrieve all available Sequence Ontology terms from the file. Parse out mappings between feature keys and SO. Parse all available features and qualifiers from the FT definition.      This is ugly and parses it straight out of the HTML but this is much easier     than trying to get it from the specs.