makeFigures.py

# -*- coding: utf-8 -*-
"""
Created on Thu Feb  6 09:46:15 2020

@author: Tiago
"""
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
import csv
from rdkit import Chem
from rdkit.Chem import Crippen
from rdkit.Chem import Descriptors as desc
from rdkit.Chem import QED
from utils import load_config,generate2file,smiles2mol
from sascorer_calculator import SAscore
from prediction import Predictor
from keras.models import Sequential
from model import Model 

colors = ['#ff7f0e', '#1f77b4', '#d62728', '#2ca02c', '#9467bd']  # orange, blue, green, red, purple
           
def properties_violin(filepaths,labels,pred_type):
    """
    Function that extracts a pd dataframe with QED, pIC50, and SAscore from the
    specified filepaths.
    """
    properties = []
      
    for i,fname in enumerate(filepaths):
        with open(filepaths[i], 'r') as f:
            reader = csv.reader(f)
        
            it = iter(reader)
#            next(it, None)  # skip first item.    
            for row in it:
                    if pred_type == 'pIC50':
                        properties.append([labels[i],'pIC50 for KOR',float(row[1])]) 
                    if i != 0:
                        properties.append([labels[i],'SA score',float(row[2])])
                        try:
                            mol = Chem.MolFromSmiles(row[0])
                            q = QED.qed(mol)
    #                        x, y = desc.MolWt(mol), Crippen.MolLogP(mol)
    #                        properties.append([labels[i],'Molecular weight',x])
    #                        properties.append([labels[i],'logP',y])
                            properties.append([labels[i],'QED',q])
                            
                        except: 
                            print("Non-Canonical SMILES: " + row[0])
                    else:
                        
                        try:
                            mole = smiles2mol(row[0])
                            prediction_sas = SAscore(mole)
                            properties.append([labels[i],'SA score',float(prediction_sas[0])])
                            mol = Chem.MolFromSmiles(row[0])
                            q = QED.qed(mol)
    #                        x, y = desc.MolWt(mol), Crippen.MolLogP(mol)
    #                        properties.append([labels[i],'Molecular weight',x])
    #                        properties.append([labels[i],'logP',y])
                            properties.append([labels[i],'QED',q])
                        except: 
                            print("Non-Canonical SMILES: " + row[0])
                    
    df = pd.DataFrame(properties, columns=['Sets', 'Property', 'Value'])
    return df

def properties_mw_logp(filepaths): 
    """
    Function used to extract a pd dataframe with MW and logP from the specified 
    filepaths.
    """   
    properties = []
      
    for i,fname in enumerate(filepaths):
        with open(filepaths[i], 'r') as f:
            reader = csv.reader(f)
        
            it = iter(reader)
            if not ("generated" in fname):
                for row in it:
                    try:
                        properties.append([float(row[2]),float(row[3]),i])
                    except:
                        print("")
            else:
                for row in it:
                    try:
                        mol = Chem.MolFromSmiles(row[0])
                        x, y = desc.MolWt(mol), Crippen.MolLogP(mol)
                        properties.append([x,y,i])
                    except: 
                        print("Non-Canonical SMILES: " + row[0])

    df = pd.DataFrame(properties[2000:2355], columns=['MW', 'logP', 'Label'])
    return df


def violin_plot(pred_identifier):
    """ 
    Function that graphs a violin plot for the physicochemical properties comparison.
    It compares molecules generated by unbiased and optimized models
    """
    config_file = 'configReinforce.json' 
    configReinforce,exp_time=load_config(config_file)
    
    if pred_identifier == 'pIC50':
         # Load the predictor model
        predictor = Predictor(configReinforce,'kor')
    else:
        predictor = None
    
    biased_generator = Sequential()
    biased_generator = Model(configReinforce)
    unbiased_generator = Sequential()
    unbiased_generator = Model(configReinforce)
    biased_generator.model.load_weights(configReinforce.model_name_biased+".h5")
    unbiased_generator.model.load_weights(configReinforce.model_name_unbiased)

    # generate with unbiased
    generate2file(predictor,unbiased_generator,configReinforce,50,True)
    generate2file(predictor,biased_generator,configReinforce,50,False)
#     
    plt.figure(figsize=(15, 5))

    sns.set(rc={"axes.facecolor":"white",
            "axes.grid":False,
            'axes.labelsize':20,
            'figure.figsize':(20.0, 10.0),
            'xtick.labelsize':15,
            'ytick.labelsize':15})
    
#    sns.set(style="white", palette="colorblind", color_codes=True)
    df = properties_violin(['Generated/smiles_prop_original.smi','Generated/smiles_prop_biased.smi'], ['Original generator', 'Fine-tuned generator'],pred_identifier)
    sns.violinplot(x='Property', y='Value', hue='Sets', data=df, linewidth=1, split=True, bw=1, legend = False)
    sns.despine(left=True)
    plt.ylim([-3, 15])
#    viol_plot.ax.legend(loc=2)

def mw_logp_plot():  
    """
    Function that graphs a chemical space comparison based on logP ~ MW
    """
    fig = plt.figure(figsize=(12, 12))
    lab = ['Chembl Dataset', 'Fine-tuned dataset']
    ax1 = fig.add_subplot(221)
    df = properties_mw_logp(['data/tyk_clean.csv', 'data/generated.smi'])

    group0, group1 = df[df.Label == 0], df[df.Label== 1]
    ax1.scatter(group0.MW, group0.logP, s=10, marker='o', label=lab[0], c='', edgecolor=colors[1])
    ax1.scatter(group1.MW, group1.logP, s=10, marker='o', label=lab[1], c='', edgecolor=colors[3])
    ax1.set(ylabel='LogP', xlabel='Molecular Weight')
    ax1.legend(loc='lower right')

    
def performance_barplot():
    """
    Function that graphs two bar plots with the QSAR model metrics (MSE and Q2)
    """
    objects = ('RNN', 'SVM', 'FCNN', 'RF', 'KNN')
    y_pos = np.arange(len(objects))
    
    performance_mse = [0.0224, 0.0306, 0.0397, 0.0255, 0.037826]
    plt.bar(y_pos, performance_mse, align='center', alpha=0.5, color=['blue' ,'yellow', 'red', 'green','orange'],  edgecolor='black')
    plt.xticks(y_pos, objects)
    plt.ylabel('MSE')
    plt.title('QSAR models evaluation: Mean Squared Error')
    plt.show()
    

    performance_q2 = [80.12, 74.857, 66.004, 78.96, 69.144]   
    plt.bar(y_pos, performance_q2, align='center', alpha=0.5, color=['blue' ,'yellow', 'red', 'green','orange'],  edgecolor='black')
    plt.xticks(y_pos, objects)
    plt.ylabel('Q2 Score')
    plt.title('QSAR models evaluation: Q-squared')
    plt.show()
    
    performance_rmse = [0.1496, 0.1749, 0.1992, 0.1597, 0.1945]   
    plt.bar(y_pos, performance_rmse, align='center', alpha=0.5, color=['blue' ,'yellow', 'red', 'green','orange'],  edgecolor='black')
    plt.xticks(y_pos, objects)
    plt.ylabel('RMSE')
    plt.title('QSAR models evaluation: Root Mean Squared Error')
    plt.show()
    
    performance_ccc = [89.862, 84.43, 80.1, 88.2, 82.34]   
    plt.bar(y_pos, performance_ccc, align='center', alpha=0.5, color=['blue' ,'yellow', 'red', 'green','orange'],  edgecolor='black')
    plt.xticks(y_pos, objects)
    plt.ylabel('CCC')
    plt.title('QSAR models evaluation: Concordance Correlation Coefficient')
    plt.show()

def regression_plot(y_true,y_pred):
    """
    Function that graphs a scatter plot and the respective regression line to 
    evaluate the QSAR models.
    Parameters
    ----------
    y_true: True values from the label
    y_pred: Predictions obtained from the model
    Returns
    -------
    This function returns a scatter plot.
    """
    fig, ax = plt.subplots()
    ax.scatter(y_true, y_pred)
    ax.plot([np.min(y_true), np.max(y_true)], [np.min(y_true), np.max(y_true)], 'k--', lw=4)
    ax.set_xlabel('True')
    ax.set_ylabel('Predicted')
    plt.show()
    

def plot_MO():
    """
    This function plots a scatter diagram with multi-objective results 
    Parameters
    ----------
    cumulative_rewards_qed: List with the previous averaged rewards for QED property
    cumulative_rewards_kor: List with the previous averaged rewards for KOR property
    cumulative_rewards: List with the previous scalarized combined rewards
    previous_weights: List with the all the previous choosen weights
    Returns
    -------
    This function plots the scatter diagram with all results for each weights assignment
    """

    # Linear scalarization
#    cumulative_rewards_qed = [1.2076, 1.203, 1.2106, 1.2071, 1.1976, 1.196, 1.1918, 1.1957, 1.1913, 1.1589, 1.199]
#    cumulative_rewards_kor = [2.328, 2.467, 2.498, 2.693, 2.696, 2.638, 2.597, 2.7578, 2.6952,  2.874, 2.676]	
#    previous_weights = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]

	# Weight search - linear scalarization
#    cumulative_rewards_qed = [1.2076, 1.199, 1.196, 1.1914, 1.2178, 1.199, 1.197, 1.209, 1.1989, 1.2073,1.2056]
#    cumulative_rewards_kor = [2.328, 2.676, 2.638, 2.7067, 2.468, 2.453, 2.574, 2.54, 2.5288, 2.5651,2.597]	
#    previous_weights = [0, 1, 0.5, 0.75, 0.25, 0.375, 0.438, 0.406, 0.391, 0.383, 0.379]
    
    # Chebyshev scalarization
#    cumulative_rewards_qed = [1.2076, 1.2129, 1.212, 1.2111, 1.19456, 1.2184, 1.1984, 1.1985, 1.20449, 1.187,1.1796 ]
#    cumulative_rewards_kor = [2.328, 2.383, 2.453, 2.4498, 2.6202, 2.5195, 2.6927, 2.654, 2.63835, 2.778,2.783]	
#    previous_weights = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]    
    
    plt.clf()

    # plot the chart
    plt.scatter(cumulative_rewards_kor,cumulative_rewards_qed)

    # zip joins x and y coordinates in pairs
    i = 0
    for x,y in zip(cumulative_rewards_kor,cumulative_rewards_qed):
        
        label = '(' + "{:.2f}".format(previous_weights[i]) +", "+ "{:.2f}".format(1-previous_weights[i]) + ')' 

        # this method is called for each point
        plt.annotate(label, # this is the text
                     (x,y), # this is the point to label
                    textcoords="offset points", # how to position the text
                     xytext=(0,10), # distance from text to points (x,y)
                     ha='center') # horizontal alignment can be left, right or center
        i+=1
    plt.xlabel("Mean KOR reward")
    plt.ylabel("Mean QED reward")
    #plt.xticks(np.arange(0,10,1))
    #plt.yticks(np.arange(0,5,0.5))

    plt.show()
    
    
def compute_hypervolume():
    """
    This function computes the hypervolume indicator based on sampled solutions
    obtained by three scalarizations techniques
    """
    refPoint = [2.4,1.15]
    
 #   # Linear
#    obj1 = [2.498,2.874,2.7578,2.693]
#    obj2 = [1.2106,1.1589,1.1957,1.2071]
    
  #   Linear ws
    obj1 = [2.676,2.7067,2.468,2.5651,2.597,2.54]
    obj2 = [1.199,1.1914,1.2178,1.2073,1.2056,1.209]
 
# # Chebyshev
#    obj1 = [2.5195,2.783,2.778,2.63835,2.6927]
#    obj2 = [1.2184,1.1796,1.187,1.20449,1.1984]
    
    
    # code
    obj1.sort()
    obj2.sort(reverse = True)
    hv = 0
    
    for i,solution1 in enumerate(obj1):
        if i == 0:
            hv = (solution1-refPoint[0])*(obj2[i]-refPoint[1])
        else:
            hv+= (solution1-obj1[i-1])*(obj2[i]-refPoint[1])
    print("Hypervolume: " + str(hv))
    
    
def main():
    """
    Main of function that plots several figures to help the interpretation of the 
    results obtained from both the Predictor and Generator.
    """
    predictor = 'pIC50' # pIC50 or sas          
    y_true = [1, 2, 3,4,5,2,2,1,3]
    y_pred = [1.1, 1.6, 3, 4.1, 4.5, 1.8, 2,1.6,3.1]
#    regression_plot(y_true,y_pred)
    performance_barplot()
#    violin_plot(predictor)
#    mw_logp_plot()
#    plot_MO()
#    plot_rapido()
#    compute_hypervolume()
if __name__ == '__main__':
    main()