<jupyter_start><jupyter_text>Using ML to Predict Bioactivity of Drugs**Analysis for coronavirus drug data retrieved from [*ChEMBL Database*](https://www.ebi.ac.uk/chembl/)**Submitted by- <NAME> 16235  Installing library which will be used to retreive data<jupyter_code>try:
    import chembl_webresource_client
    from chembl_webresource_client.new_client import new_client
except:
    !pip install chembl_webresource_client
    import chembl_webresource_client
    from chembl_webresource_client.new_client import new_client
import pandas as pd<jupyter_output><empty_output><jupyter_text>retrieving target based data<jupyter_code>target = new_client.target
target_query = target.search('coronavirus')
targets = pd.DataFrame.from_dict(target_query)
targets<jupyter_output><empty_output><jupyter_text>Will be using *SARS coronavirus 3C-like proteinase* Bioactivity Data<jupyter_code>selected_target = targets.target_chembl_id[4]
selected_target
activity = new_client.activity
res = activity.filter(target_chembl_id=selected_target).filter(standard_type="IC50")
data = pd.DataFrame.from_dict(res)
data.columns
data.head()
data_1 = data[data.standard_value.notna()] #removing unavailable data
data.shape<jupyter_output><empty_output><jupyter_text>Preprocesssing  **Labeling compounds as either being active, inactive or intermediate**The bioactivity data is in the IC50 unit. Compounds having values of less than 1000 nM will be considered to be **active** while those greater than 10,000 nM will be considered to be **inactive**. As for those values in between 1,000 and 10,000 nM will be referred to as **intermediate**.<jupyter_code>bioactivity_class = []
for i in data_1.standard_value:
  if float(i) >= 10000:
    bioactivity_class.append("inactive")
  elif float(i) <= 1000:
    bioactivity_class.append("active")
  else:
    bioactivity_class.append("intermediate")
# USING ONLY RELEVANT COLUMNS TO US
selection = ['molecule_chembl_id', 'canonical_smiles', 'standard_value']
data_2 = data_1[selection]
data_3 = pd.concat([data_2,pd.Series(bioactivity_class)], axis=1)
data_3.rename(columns={0:'Bio_class'},inplace=True)
data_3.head()<jupyter_output><empty_output><jupyter_text>Downloading other Liberaries to create descriptors for Data analysis<jupyter_code># uncomment this if You are running in colab. It is for rdkit installation on google colab.
#! wget https://repo.anaconda.com/miniconda/Miniconda3-py37_4.8.2-Linux-x86_64.sh
#! chmod +x Miniconda3-py37_4.8.2-Linux-x86_64.sh
#! bash ./Miniconda3-py37_4.8.2-Linux-x86_64.sh -b -f -p /usr/local
#! conda install -c rdkit rdkit -y
#import sys
#sys.path.append('/usr/local/lib/python3.7/site-packages/')
# This was to check whether we have desired csv present in directory or not.
#import os
#os.getcwd()
#os.chdir('Colab Notebooks')
#os.listdir()
import numpy as np
from rdkit import Chem
from rdkit.Chem import Descriptors, Lipinski
# code used from dataprof 
def descriptors(smiles, verbose=False):

    moldata= []
    for elem in smiles:
        mol=Chem.MolFromSmiles(elem) 
        moldata.append(mol)
       
    baseData= np.arange(1,1)
    i=0  
    for mol in moldata:        
       
        desc_MolWt = Descriptors.MolWt(mol)
        desc_MolLogP = Descriptors.MolLogP(mol)
        desc_NumHDonors = Lipinski.NumHDonors(mol)
        desc_NumHAcceptors = Lipinski.NumHAcceptors(mol)
        # other descriptors are not relatable for Bioactivity predictions, SO 
        # sticking with default 4 of these.

           
        row = np.array([desc_MolWt,
                        desc_MolLogP,
                        desc_NumHDonors,
                        desc_NumHAcceptors])   
    
        if(i==0):
            baseData=row
        else:
            baseData=np.vstack([baseData, row])
        i=i+1      
    
    columnNames=["MW","LogP","NumHDonors","NumHAcceptors"]   
    descriptors = pd.DataFrame(data=baseData,columns=columnNames)
    
    return descriptors
data_desc = descriptors(data_3.canonical_smiles)
data_desc.head()
data_combined = pd.concat([data_3,data_desc], axis = 1)
data_combined.head()<jupyter_output><empty_output><jupyter_text>Converting standard_value IC$_{50}$ (have nM-nanoMolar unit) to Logrithmic pIC$_{50}$.<jupyter_code>data_combined['standard_value'] = data_combined['standard_value'].astype(float)
data_combined.standard_value.describe()
def pIC50(input):
    pIC50 = []

    for i in input['standard_value_norm']:
        molar = i*(10**-9)             # Converts nM to M
        pIC50.append(-np.log10(molar))

    input['pIC50'] = pIC50
    x = input.drop('standard_value_norm', 1)
        
    return x

def norm_value(input):
    norm = []

    for i in input['standard_value']:
        if i > 100000000:
          i = 100000000
        norm.append(i)

    input['standard_value_norm'] = norm
    x = input.drop('standard_value', 1)
        
    return x<jupyter_output><empty_output><jupyter_text>Since moles can't be negative, values more than 100 million are capped to 100 million<jupyter_code>data_normed = norm_value(data_combined)
data_final = pIC50(data_normed)
for i in data_final.pIC50:
    if i<0:
        print(i)
data_final.head()<jupyter_output><empty_output><jupyter_text>We don't want logrithmic negative values therefore capping the max of IC50 data to 100 million<jupyter_code>data_final.to_csv('data_final.csv',index=False)
#data_final = pIC50(df_norm)
#data_final.head()
import seaborn as sns
sns.set(style='ticks')
import matplotlib.pyplot as plt
plt.figure(figsize=(5.5, 5.5))
sns.countplot(x='Bio_class', data=data_final, edgecolor='black')
plt.xlabel('Bioactivity class', fontsize=14, fontweight='bold')
plt.ylabel('Frequency', fontsize=14, fontweight='bold')
#plt.savefig('plot_bioactivity_class.png')<jupyter_output><empty_output><jupyter_text>Scatter plot below<jupyter_code>plt.figure(figsize=(5.5, 5.5))
sns.scatterplot(x='MW', y='LogP', data=data_final, hue='Bio_class', size='pIC50', edgecolor='black', alpha=0.7)
plt.xlabel('MW', fontsize=14, fontweight='bold')
plt.ylabel('LogP', fontsize=14, fontweight='bold')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0)
#plt.savefig('plot_MW_vs_LogP.png')<jupyter_output><empty_output><jupyter_text>pIC50 box Plot<jupyter_code>plt.figure(figsize=(5.5, 5.5))
sns.boxplot(x = 'Bio_class', y = 'pIC50', data = data_final)
plt.xlabel('Bioactivity class', fontsize=14, fontweight='bold')
plt.ylabel('pIC50 value', fontsize=14, fontweight='bold')
#plt.savefig('plot_ic50.png')<jupyter_output><empty_output><jupyter_text>Comparing Active and Inactive Biactivity levels. To see these two groups are different statistically or not using ***Man-whiteny U test***<jupyter_code>data_2class = data_final[data_final.Bio_class!= 'intermediate']
data_2class.shape<jupyter_output><empty_output><jupyter_text>Statistical analysis | Man-whiteny U test~~~This test is used to analyze whether two groups are related or not. Inspired from dataprof.<jupyter_code>def mannwhitney(descriptor, verbose=False):
  # https://machinelearningmastery.com/nonparametric-statistical-significance-tests-in-python/
  from numpy.random import seed
  from numpy.random import randn
  from scipy.stats import mannwhitneyu

# seed the random number generator
  seed(1)

# actives and inactives
  selection = [descriptor, 'Bio_class']
  df = data_2class[selection]
  active = df[df.Bio_class == 'active']
  active = active[descriptor]

  selection = [descriptor, 'Bio_class']
  df = data_2class[selection]
  inactive = df[df.Bio_class == 'inactive']
  inactive = inactive[descriptor]

  selection = [descriptor, 'Bio_class']
  df = data_final[selection]
  intermediate = df[df.Bio_class == 'intermediate']
  intermediate = intermediate[descriptor]

# compare samples
  stat, p = mannwhitneyu(active, inactive)
  #print('Statistics=%.3f, p=%.3f' % (stat, p))

# interpret
  alpha = 0.05
  if p > alpha:
    interpretation = 'Same distribution (fail to reject H0)'
  else:
    interpretation = 'Different distribution (reject H0)'
  
  results = pd.DataFrame({'Descriptor':descriptor,
                          'Statistics':stat,
                          'p':p,
                          'alpha':alpha,
                          'Interpretation':interpretation}, index=[0])
  
  return results
mannwhitney('pIC50')<jupyter_output><empty_output><jupyter_text>MW<jupyter_code>plt.figure(figsize=(5.5, 5.5))
sns.boxplot(x = 'Bio_class', y = 'MW', data = data_2class)
plt.xlabel('Bioactivity class', fontsize=14, fontweight='bold')
plt.ylabel('MW', fontsize=14, fontweight='bold')
#plt.savefig('plot_MW_ic50.png')
mannwhitney('MW')
plt.figure(figsize=(5.5, 5.5))
sns.boxplot(x = 'Bio_class', y = 'LogP', data = data_2class)
plt.xlabel('Bioactivity class', fontsize=14, fontweight='bold')
plt.ylabel('LogP', fontsize=14, fontweight='bold')
#plt.savefig('plot_LogP_IC50.png')
mannwhitney('LogP')
plt.figure(figsize=(5.5, 5.5))
sns.boxplot(x = 'Bio_class', y = 'NumHDonors', data = data_2class)
plt.xlabel('Bioactivity class', fontsize=14, fontweight='bold')
plt.ylabel('NumHDonors', fontsize=14, fontweight='bold')
#plt.savefig('plot_NumHDonors_IC50.png')
mannwhitney('NumHDonors')
plt.figure(figsize=(5.5, 5.5))
sns.boxplot(x = 'Bio_class', y = 'NumHAcceptors', data = data_2class)
plt.xlabel('Bioactivity class', fontsize=14, fontweight='bold')
plt.ylabel('NumHAcceptors', fontsize=14, fontweight='bold')
#plt.savefig('plot_NumHAcceptors_IC50.png')
mannwhitney('NumHAcceptors')<jupyter_output><empty_output><jupyter_text>Interpretation of Results From Box Plots**pIC50 values**Taking a look at pIC50 values, the actives and inactives displayed statistically significant difference, which is to be expected since threshold values   (IC50  10,000 nM = Inactives, corresponding to pIC50 > 6 = Actives and pIC50 < 5 = Inactives) were used to define actives and inactives.**Lipinski Descriptors**Of the 4 Lipinski's descriptors (MW, LogP, NumHDonors and NumHAcceptors), only **LogP** exhibited no difference between the actives and inactives while the other 3 descriptors **(MW, NumHDonors and NumHAcceptors)** shows statistically significant difference between actives and inactives.  Using Descriptor generatorThis piece of code is taken from [dataprof](https://www.youtube.com/watch?v=zD2focOkQ48).<jupyter_code>if 'padel.sh' in os.listdir():
    pass
else:
    ! wget https://github.com/dataprofessor/bioinformatics/raw/master/padel.zip
    ! wget https://github.com/dataprofessor/bioinformatics/raw/master/padel.sh
if 'PaDEL-Descriptor' in os.listdir():
    pass
else:
    ! unzip padel.zip
selection = ['canonical_smiles','molecule_chembl_id']
data_selection = data_final[selection]
data_selection.to_csv('molecule.smi', sep='\t', index=False, header=False)
! cat molecule.smi | head -5  #this code won't word on jupyter offline
! cat molecule.smi | wc -l<jupyter_output>133<jupyter_text>Calculate PaDEL fingerprints<jupyter_code>! bash padel.sh<jupyter_output>Processing CHEMBL187579 in molecule.smi (1/133). 
Processing CHEMBL188487 in molecule.smi (2/133). 
Processing CHEMBL185698 in molecule.smi (3/133). Average speed: 2.27 s/mol.
Processing CHEMBL426082 in molecule.smi (4/133). Average speed: 1.15 s/mol.
Processing CHEMBL365134 in molecule.smi (6/133). Average speed: 0.82 s/mol.
Processing CHEMBL187717 in molecule.smi (5/133). Average speed: 0.99 s/mol.
Processing CHEMBL187598 in molecule.smi (7/133). Average speed: 0.75 s/mol.
Processing CHEMBL190743 in molecule.smi (8/133). Average speed: 0.71 s/mol.
Processing CHEMBL365469 in molecule.smi (9/133). Average speed: 0.63 s/mol.
Processing CHEMBL191575 in molecule.smi (11/133). Average speed: 0.56 s/mol.
Processing CHEMBL188983 in molecule.smi (10/133). Average speed: 0.62 s/mol.
Processing CHEMBL194398 in molecule.smi (13/133). Average speed: 0.57 s/mol.
Processing CHEMBL370923 in molecule.smi (12/133). Average speed: 0.57 s/mol.
Processing CHEMBL196635 in molecule.smi (14/133). Average sp[...]<jupyter_text>Preparing X and Y Matrices<jupyter_code>data_X = pd.read_csv('descriptors_output.csv')
data_X = data_X.drop(columns=['Name'])
data_X.head()
data_Y = data_final['pIC50']
data_Y
final_dataset = pd.concat([data_X,data_Y], axis=1)
final_dataset.head()
final_dataset.to_csv('Final_dataset.csv', index=False)<jupyter_output><empty_output><jupyter_text>Training ModelNow we have our cleaned and preprocessed data. So we will train few models and use the best performing one. Input and output features<jupyter_code>X,Y = data_X, data_Y
X.shape, Y.shape
from sklearn.model_selection import train_test_split
from sklearn import linear_model, svm
from sklearn.ensemble import RandomForestRegressor 
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2)
X_train.shape, Y_train.shape, X_test.shape, Y_test.shape
classifiers = [
    svm.SVR(),
    linear_model.SGDRegressor(),
    linear_model.BayesianRidge(),
    linear_model.LassoLars(),
    linear_model.ARDRegression(),
    linear_model.PassiveAggressiveRegressor(),
    linear_model.TheilSenRegressor(),
    linear_model.LinearRegression(),
    RandomForestRegressor()]

for item in classifiers:
    #print(item,'\n')
    clf = item
    clf.fit(X_train, Y_train)
    #print(clf.predict(X_test),'\n')
    print(f" r2 score of {str(item).split('(')[0]} is {clf.score(X_test, Y_test):.4} \n")<jupyter_output>r2 score of SVR is 0.2464 

 r2 score of SGDRegressor is 0.1369 

 r2 score of BayesianRidge is 0.3175 

 r2 score of LassoLars is -0.06408 

 r2 score of ARDRegression is 0.0175 

 r2 score of PassiveAggressiveRegressor is -0.1793 

 r2 score of TheilSenRegressor is -0.6508 

 r2 score of LinearRegression is -2.684e+17 

 r2 score of RandomForestRegressor is 0.005313<jupyter_text>From above testing it is clear that [**BayesianRidge**](https://scikit-learn.org/stable/auto_examples/linear_model/plot_bayesian_ridge.htmlbayesian-ridge-regression) method is having highest score, although it is also not good.<jupyter_code>clf = linear_model.BayesianRidge()
clf.fit(X_train,Y_train)
Y_pred = clf.predict(X_test)
sns.set(color_codes=True)
sns.set_style("white")

ax = sns.regplot(Y_test, Y_pred, scatter_kws={'alpha':0.4})
ax.set_xlabel('Experimental pIC50', fontsize='large', fontweight='bold')
ax.set_ylabel('Predicted pIC50', fontsize='large', fontweight='bold')
ax.set_xlim(0, 12)
ax.set_ylim(0, 12)
ax.figure.set_size_inches(5, 5)
plt.show<jupyter_output>/usr/local/lib/python3.6/dist-packages/seaborn/_decorators.py:43: FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.
  FutureWarning
{'fraction_non_alphanumeric': 0.11882873680626489, 'fraction_numerical': 0.03663602315287708, 'mean_word_length': 5.0386513157894735, 'pattern_counts': {'":': 0, '<': 73, '<?xml version=': 0, '>': 74, 'https://': 7, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '589307', 'license': 'Apache-2.0', 'n_tokens_mistral': 5580, 'n_tokens_neox': 5063, 'n_words': 1154, 'path': 'Bioactivity/Using_ML_to_predict_bioActivity_of_covid_drugs.ipynb', 'repo_name': 'axxe007/ML', 'size': '14692'}