# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: rdkit-env
#     language: python
#     name: rdkit-env
# ---

# +
import wget
import zipfile
import re
import tqdm
import os
import numpy as np

from rdkit import Chem
from rdkit.Chem.Descriptors import MolWt
from rdkit.Chem.MolStandardize import rdMolStandardize
enumerator = rdMolStandardize.TautomerEnumerator()


import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
import matplotlib as mpl
from pylab import cm


# +
########################################################### FUNCTIONS BY ITERATIONS #######################################

def get_stable_tatuomer(Smiles):
    """ Get Sable Tautoer """
    tautomer = ""
    try:
        tautomer = Chem.CanonSmiles(Chem.MolToSmiles(enumerator.Canonicalize(Chem.MolFromSmiles(Smiles))))
    except:
        print(f"Could not find tatuomer for {Smiles}")
        
        
    return tautomer
    
def network_iterations(path,Tautomers=False):
    """ Generates a dictionary that saves molecules by Iteration"""
    Network = dict()    # Define an Empty Dictionary
    Loaded = 0          # Number of molecules loaded
    ### Open file with context handler ###
    with open(path, "r") as handler:
        for line in  handler:                     # Iterate over the lines
            try:
                iteration,smiles = line.split()      # Splite every line: The format is Iteration,moleule  e.g. G1    ccc1c3
                try:
                    if Tautomers:
                        CanonicalSmiles = get_stable_tatuomer(smiles)
                    else:   
                        CanonicalSmiles = Chem.CanonSmiles(smiles)
                except:
                    pass
                
                if CanonicalSmiles != "":
                    try:

                        Network[iteration].append(CanonicalSmiles) # Append molecule to iteration list inside Network dictionary
                        Loaded +=1

                    except:
                        Network[iteration] = list()           # If an exception occurs then iteration list does not exxists, create it
                        Network[iteration].append(CanonicalSmiles) # Append molecule to iteration list inside Network dictionary
                        Loaded +=1
            except:
                pass
    
    print(f"Iterations: {Network.values().__len__()}, Total Moleules: {Loaded}")
        
    return Network


                
def get_max_mass_iterations(Network):
    MaxMass = dict()                                                             # Initiate Dictioanry
    for Iteration,SmilesList in zip(Network.keys(),Network.values()):            # Iterate over Network
        Mass = max([MolWt(Chem.MolFromSmiles(Smiles)) for Smiles in SmilesList]) # Get Maximum mass at each iteration
        MaxMass[Iteration] = Mass                                                # Append it to dictionary
        print(f"Max Molecular Weight for iteration {Iteration}: {Mass}")
        
    
    return MaxMass

###################################################### Functions without database iterations

def database(path,MaxMass,Tautomers=False):
    """ Generates a list that saves smiles strings """
    Ignored = 0
    Network = list()                            # Define an Empty Dictionary
    Limit = max(MaxMass.values())               # Maxmimum molecular weight of network 
    ### Open file with context handler ###
    with open(path, "r") as handler:            # Iterate over handler
        for line in handler:                    # Iterate over the lines
            smiles = line                       # Splite every line
            try:
                if Tautomers:
                    CanonicalSmiles = get_stable_tatuomer(smiles)
                else:   
                    CanonicalSmiles = Chem.CanonSmiles(smiles) # Convert smiles to canonical smiles
                # Check that molecules is not heavier than heavies molecule in network (we will compare apples to apples)
                if MolWt(Chem.MolFromSmiles(CanonicalSmiles)) <= Limit: 
                    if re.match("^[CcOoNn]",CanonicalSmiles): # If Molecules starts with C,N,O
                        Network.append(CanonicalSmiles)       # Append molecule to network if mass is lower than limit
                else:
                    Ignored +=1 # If not add 1 to ignore counter
                    
            except:
                Ignored +=1   # If there is no canonical smiels, than there must be a typo in string, ingore it 
            
    
    print(f"{Ignored} molecules were ignored out of {len(Network)+Ignored}. Leaving a total of  {len(Network)} Molecules.")         
    return Network


#### Compare program ###
def compare_networks_database(mod,database,verbose=True,plot=True):
    matched = dict()
    missed = dict()

    for mol2 in database:
        for iterations in  mod.keys():
            #print(f" Currently working on Iteration: {iterations}")
            total_set = set(mod[iterations])
            for mol in total_set:
                #### Capture Matches ####
                if mol == mol2:
                    try:
                        matched[iterations].add(mol) # Append molecule to iteration list inside Network dictionary

                    except:
                        matched[iterations] = set()    # If an exception occurs then iteration list does not exxists, create it
                        matched[iterations].add(mol)   # Append molecule to iteration list inside Network dictionary
    
    
    for iterations in  mod.keys():
        total_set    = set(mod[iterations])         # Define Universe for iteration
        matched_set  = matched[iterations]          # Get matched set 
        missed[iterations] = total_set -matched_set # Get Intersection
        
    
    if verbose:
        matched_list = list()
        matched_database =list()
        unmatched_list = list()
        print()
        print("Summary (Unique Molecules)")
        print("-".center(75,"-"))
        for iterations in  mod.keys():
            print(f"Iteration {iterations}")  
            matched_set = matched[iterations]
            total_set   = set(mod[iterations])
            matched_percent = len(matched_set)/len(total_set)*100
            match_d=  len(matched_set)/len(database)*100
            matched_database.append(match_d)
            print(f"Molecules Matches: {len(matched_set)} of total {len(total_set)}, percentage: {matched_percent} %")
            matched_list.append(matched_percent)
            unmatched = missed[iterations]
            missed_percent = len(unmatched)/len(total_set)*100
            #{len(missed[iterations])}
            print(f"Molecules Missed: {len(unmatched)} of total {len(total_set)} , percent: {missed_percent} %")
            unmatched_list.append(missed_percent)
            print("-".center(75,"-"))
            print()
            
    if plot:
        
      
        Matches_sorted = np.log10(np.array([ len(a[1]) for a in  sorted(matched.items())]))
        Iterations = range(len(Matches_sorted))
        labels =  [ "G" + str(i) for i in Iterations]
        Missed_sorted = np.log10(np.array([ len(a[1]) for a in  sorted(missed.items())]))
        Total_sorted  = Matches_sorted + Missed_sorted
        Network_sorted = np.log10(np.array([ len(a[1]) for a in  sorted(Network.items())]))

        plt.figure(figsize=(10,7))
        plt.title("Matched/Unmatched molecules through iterations")
        plt.xlabel("Iterations",fontsize=12)
        plt.ylabel("Number of Molecules (log10)",fontsize=12)
        plt.plot(Matches_sorted,"--g",label="Matched molecules")
        plt.plot(Network_sorted,"--y",label="Network molecules")
        i=0
        for match,percent in zip(Matches_sorted,  matched_database):
            plt.text(i-0.1,match+0.2,str(round(percent,2)) + "%",fontsize=12)
            i+=1
        plt.xticks(Iterations,labels)
        plt.legend(loc="best")
        plt.grid()
        plt.show()


    return matched,missed,matched_list,unmatched_list
# -

path = "Datasets/glucose_degradation_output2.txt"
Network =  network_iterations(path,Tautomers=False)

path = "Datasets/glucose_degradation_output2.txt"
Network =  network_iterations(path,Tautomers=True)

MaxMass =  get_max_mass_iterations(Network)

path =  "Datasets/HarrisonDataAugust2019_molescules.csv"
Database = database(path,MaxMass,Tautomers=False)

# Font
mpl.rcParams['font.family'] ='Times New Roman'
plt.rcParams['font.size'] = 16
plt.rcParams['axes.linewidth'] = 1


matched,missed,matched_list,unmatched_list = compare_networks_database(Network ,Database,verbose=True)

len(matched_l)

plt.bar(list(MaxMass.keys()),list(MaxMass.values()), width=0.8, bottom=None, align='center')
plt.show()

MaxMass.values()


def DownloadHMDB(Verbose):
    """Downloads an instance of the HMDB.sdf Database if a local copy is not avialble"""
    if Verbose:
        print("Oops, you Don't have a local coppy of HMDB.sdf")
        print("A new instance of the Database will be dowloaded for you")
    url = "https://hmdb.ca/system/downloads/current/structures.zip"
    root = os.getcwd()
    filename = wget.download(url)
    with zipfile.ZipFile(filename, 'r') as zip_ref:
        zip_ref.extractall(os.getcwd())
    if Verbose:
        print("Almost Done...")
    os.rename(r'structures.sdf',r'HMDB.sdf')
    os.chdir(root)
    if Verbose:
         print("Done")


DownloadHMDB(Verbose=True)     # Download sdf
  

path = "HMDB.sdf"
hmdb=Chem.SDMolSupplier(path)  # Load sdf

hmdbsmiles = []
for mol in hmdb:
    try:
        hmdbsmiles.append(Chem.MolToSmiles(mol)) 
    except:
        pass
    

len(hmdbsmiles)


def database_smiles(SMILES,MaxMass,Tautomers=False):
    """ Generates a list that saves smiles strings """
    Ignored = 0
    Network = list()                            # Define an Empty Dictionary
    Limit = max(MaxMass.values())               # Maxmimum molecular weight of network 
    ### Open file with context handler ###
    for smiles in SMILES:                    # Iterate over the lines
        try:
            if Tautomers:
                CanonicalSmiles = get_stable_tatuomer(smiles)
            else:   
                CanonicalSmiles = Chem.CanonSmiles(smiles) # Convert smiles to canonical smiles
            # Check that molecules is not heavier than heavies molecule in network (we will compare apples to apples)
            if MolWt(Chem.MolFromSmiles(CanonicalSmiles)) <= Limit: 
                if re.match("^[CcOoNn]",CanonicalSmiles): # If Molecules starts with C,N,O
                    Network.append(CanonicalSmiles)       # Append molecule to network if mass is lower than limit
            else:
                Ignored +=1 # If not add 1 to ignore counter

        except:
            Ignored +=1   # If there is no canonical smiels, than there must be a typo in string, ingore it 


    print(f"{Ignored} molecules were ignored out of {len(Network)+Ignored}. Leaving a total of  {len(Network)} Molecules.")         
    return Network


Database2 = database_smiles(hmdbsmiles,MaxMass,Tautomers=False)

matched,missed,matched_list,unmatched_list = compare_networks_database(Network ,Database2,verbose=True)





StarcoderdataJupyterScriptsDedupFiltered