<jupyter_start><jupyter_text>Load Python Modules<jupyter_code>import warnings
warnings.simplefilter(action='ignore')
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals

import os
import numpy as np
import pandas as pd
import shutil
import deepchem as dc
import json

from rdkit import Chem
from sklearn.metrics import roc_auc_score
from  sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import RandomForestClassifier
from tqdm import tnrange, tqdm_notebook
os.chdir('/home/yuke/PythonProject/DrugEmbedding/')
from decode import *
"""
SIDER dataset loader.
"""
def load_sider(featurizer='ECFP', split='index', seed=0):
    #current_dir = os.path.dirname(os.path.realpath(__file__))

    # Load SIDER dataset
    print("About to load SIDER dataset.")
    dataset_file = os.path.join("./data/sider/deepchem/sider.csv.gz")
    dataset = dc.utils.save.load_from_disk(dataset_file)
    print("Columns of dataset: %s" % str(dataset.columns.values))
    print("Number of examples in dataset: %s" % str(dataset.shape[0]))

    # Featurize SIDER dataset
    print("About to featurize SIDER dataset.")
    if featurizer == 'ECFP':
        featurizer = dc.feat.CircularFingerprint(size=1024)
    elif featurizer == 'GraphConv':
        featurizer = dc.feat.ConvMolFeaturizer()

    SIDER_tasks = dataset.columns.values[1:].tolist()
    print("SIDER tasks: %s" % str(SIDER_tasks))
    print("%d tasks in total" % len(SIDER_tasks))

    loader = dc.data.CSVLoader(
      tasks=SIDER_tasks, smiles_field="smiles", featurizer=featurizer)
    dataset = loader.featurize(dataset_file)
    print("%d datapoints in SIDER dataset" % len(dataset))

    # Initialize transformers
    transformers = [
      dc.trans.BalancingTransformer(transform_w=True, dataset=dataset)]
    print("About to transform data")
    for transformer in transformers:
        dataset = transformer.transform(dataset)

    splitters = {'index': dc.splits.IndexSplitter(),
               'random': dc.splits.RandomSplitter(),
               'scaffold': dc.splits.ScaffoldSplitter()}
    splitter = splitters[split]
    train, valid, test = splitter.train_valid_test_split(dataset, seed=seed)

    return SIDER_tasks, (train, valid, test), transformers<jupyter_output><empty_output><jupyter_text>Load Model<jupyter_code>#seed = 100, 99, 98, 97, 95
seed = 101
np.random.seed(seed)
exp_dir = './experiments/KDD/kdd_010'
checkpoint = 'checkpoint_epoch110.model'
#exp_dir = './experiments/EXP_TASK/exp_task_010'
#checkpoint = 'checkpoint_epoch100.model'
config_path = os.path.join(exp_dir, 'configs.json')
checkpoint_path = os.path.join(exp_dir, checkpoint)

with open(config_path, 'r') as fp:
    configs = json.load(fp)
fp.close()

configs['checkpoint'] = checkpoint
configs
model = load_model(configs)
print(model)
exp_dir = './experiments/KDD/kdd_009'
checkpoint = 'checkpoint_epoch110.model'
#exp_dir = './experiments/EXP_TASK/exp_task_009'
#checkpoint = 'checkpoint_epoch100.model'
config_path = os.path.join(exp_dir, 'configs.json')
checkpoint_path = os.path.join(exp_dir, checkpoint)

with open(config_path, 'r') as fp:
    configs_e = json.load(fp)
fp.close()

configs_e['checkpoint'] = checkpoint
configs_e
model_e = load_model(configs_e)
print(model_e)<jupyter_output>EVAE(
  (encoder_rnn): GRU(49, 512, batch_first=True)
  (decoder_rnn): GRU(49, 512, batch_first=True)
  (hidden2mean): Linear(in_features=512, out_features=64, bias=True)
  (hidden2logv): Linear(in_features=512, out_features=64, bias=True)
  (latent2hidden): Linear(in_features=64, out_features=512, bias=True)
  (outputs2vocab): Linear(in_features=512, out_features=49, bias=True)
  (RECON): NLLLoss()
)<jupyter_text>Process Data<jupyter_code>sider_tasks, datasets, transformers = load_sider(split='random', seed=seed)
train_dataset, valid_dataset, test_dataset = datasets
# side effects labels
train_y = train_dataset.y
valid_y = valid_dataset.y
test_y = test_dataset.y
# Finger Print
train_fp = train_dataset.X
valid_fp = valid_dataset.X
test_fp = test_dataset.X
train_smi = train_dataset.ids
valid_smi = valid_dataset.ids
test_smi = test_dataset.ids
def smi2mean(configs, model, smi_lst):
    mean_lst = []
    for i in tnrange(len(smi_lst)):
        # convert to canonical form
        smi = smi_lst[i]
        try:
            smi_can = Chem.MolToSmiles(Chem.MolFromSmiles(smi))
        except:
            smi_can = smi
        mean, _ = smiles2mean(configs, smi_can, model)
        mean_lst.append(mean.squeeze().cpu().detach().numpy())
    return np.array(mean_lst)
train_mu = smi2mean(configs, model, train_smi)
valid_mu = smi2mean(configs, model, valid_smi)
test_mu = smi2mean(configs, model, test_smi)
train_mu_e = smi2mean(configs_e, model_e, train_smi)
valid_mu_e = smi2mean(configs_e, model_e, valid_smi)
test_mu_e = smi2mean(configs_e, model_e, test_smi)<jupyter_output><empty_output><jupyter_text>Custom Distance Functions<jupyter_code>def lor_dist(z1, z2):
    m = z1*z2
    lor_prod = m[1:].sum() - m[0]
    x = - lor_prod
    x = np.where(x<1.0, 1.0+1e-6, x)
    return np.log(x + np.sqrt(x**2 - 1))<jupyter_output><empty_output><jupyter_text>KNN Classifier<jupyter_code>def get_scores(y_true, y_pred):
    roc_aucs = []
    n_data, n_dim = y_true.shape
    for j in range(n_dim):
        true_label = y_true[:, j]
        pred_label = y_pred[:, j]
        roc_aucs.append(roc_auc_score(true_label, pred_label))
    return np.array(roc_aucs)
N_val = 11
W_val = 'distance'<jupyter_output><empty_output><jupyter_text>Fit and Predict  Finger Prints<jupyter_code># Finger Print Representations
neigh_fp = KNeighborsClassifier(n_neighbors=N_val, algorithm='brute', metric='rogerstanimoto', weights=W_val)
neigh_fp.fit(train_fp, train_y)
pred_fp = neigh_fp.predict_proba(test_fp) # 27 * size(x) * 2

# covert predicted probability to an array
n_data, n_dim = test_y.shape
test_y_pred_fp = np.zeros((n_data, n_dim))
for j in range(n_dim):
    test_y_pred_fp[:, j] = pred_fp[j][:, 1]
    
#test_y_pred_fp = neigh_fp.predict(test_fp)<jupyter_output><empty_output><jupyter_text>Lorentz Embeddings<jupyter_code>neigh_lor = KNeighborsClassifier(n_neighbors=N_val, algorithm='brute', metric=lor_dist, weights=W_val)
neigh_lor.fit(train_mu, train_y)
pred_lor = neigh_lor.predict_proba(test_mu) # 27 * size(x) * 2

# covert predicted probability to an array
n_data, n_dim = test_y.shape
test_y_pred_lor = np.zeros((n_data, n_dim))
for j in range(n_dim):
    test_y_pred_lor[:, j] = pred_lor[j][:, 1]    
#test_y_pred_lor = neigh_lor.predict(test_mu)<jupyter_output><empty_output><jupyter_text>Euclidean Embeddings<jupyter_code>neigh_euc = KNeighborsClassifier(n_neighbors=N_val, algorithm='brute', weights=W_val)
neigh_euc.fit(train_mu_e, train_y)
pred_euc = neigh_euc.predict_proba(test_mu_e) # 27 * size(x) * 2

# covert predicted probability to an array
n_data, n_dim = test_y.shape
test_y_pred_euc = np.zeros((n_data, n_dim))
for j in range(n_dim):
    test_y_pred_euc[:, j] = pred_euc[j][:, 1]    
#test_y_pred_lor = neigh_lor.predict(test_mu)<jupyter_output><empty_output><jupyter_text>Evaluation<jupyter_code>get_scores(test_y, test_y_pred_fp).mean()
get_scores(test_y, test_y_pred_lor).mean()
get_scores(test_y, test_y_pred_euc).mean()<jupyter_output><empty_output><jupyter_text>RF  Finger Print<jupyter_code>rf_fp = RandomForestClassifier(class_weight="balanced", n_estimators=500)
rf_fp.fit(train_fp, train_y)
rf_pred_fp = rf_fp.predict_proba(test_fp) # 27 * size(x) * 2

# covert predicted probability to an array
n_data, n_dim = test_y.shape
test_y_rf_pred_fp = np.zeros((n_data, n_dim))
for j in range(n_dim):
    test_y_rf_pred_fp[:, j] = rf_pred_fp[j][:, 1]    
#test_y_pred_lor = neigh_lor.predict(test_mu)
get_scores(test_y, test_y_rf_pred_fp).mean()<jupyter_output><empty_output><jupyter_text>Lorentz Embeddings<jupyter_code>rf_lor = RandomForestClassifier(class_weight="balanced", n_estimators=500)
rf_lor.fit(train_mu, train_y)
rf_pred_lor = rf_lor.predict_proba(test_mu) # 27 * size(x) * 2

# covert predicted probability to an array
n_data, n_dim = test_y.shape
test_y_rf_pred_lor = np.zeros((n_data, n_dim))
for j in range(n_dim):
    test_y_rf_pred_lor[:, j] = rf_pred_lor[j][:, 1]    
#test_y_pred_lor = neigh_lor.predict(test_mu)
get_scores(test_y, test_y_rf_pred_lor).mean()<jupyter_output><empty_output><jupyter_text>Euclidean Embeddings<jupyter_code>rf_euc = RandomForestClassifier(class_weight="balanced", n_estimators=500)
rf_euc.fit(train_mu_e, train_y)
rf_pred_euc = rf_euc.predict_proba(test_mu_e) # 27 * size(x) * 2

# covert predicted probability to an array
n_data, n_dim = test_y.shape
test_y_rf_pred_euc = np.zeros((n_data, n_dim))
for j in range(n_dim):
    test_y_rf_pred_euc[:, j] = rf_pred_euc[j][:, 1]    
#test_y_pred_lor = neigh_lor.predict(test_mu)
get_scores(test_y, test_y_rf_pred_euc).mean()<jupyter_output><empty_output><jupyter_text>Deepchem Models<jupyter_code>"""
Script that trains Sklearn multitask models on the sider dataset
@Author <NAME>, <NAME>
"""
from sklearn.ensemble import RandomForestClassifier

metric = dc.metrics.Metric(dc.metrics.roc_auc_score, np.mean,
                           mode="classification")

def model_builder(model_dir):
    sklearn_model = RandomForestClassifier(class_weight="balanced", n_estimators=100)
    return dc.models.SklearnModel(sklearn_model, model_dir)

model = dc.models.SingletaskToMultitask(sider_tasks, model_builder)

# Fit trained model
model.fit(train_dataset)
model.save()

print("About to evaluate model")
train_scores = model.evaluate(train_dataset, [metric], transformers)
valid_scores = model.evaluate(valid_dataset, [metric], transformers)
test_scores = model.evaluate(test_dataset, [metric], transformers)

print("Train scores")
print(train_scores)

print("Validation scores")
print(valid_scores)

print("Test scores")
print(test_scores)<jupyter_output>About to initialize singletask to multitask model
Initializing directory for task Hepatobiliary disorders
Initializing directory for task Metabolism and nutrition disorders
Initializing directory for task Product issues
Initializing directory for task Eye disorders
Initializing directory for task Investigations
Initializing directory for task Musculoskeletal and connective tissue disorders
Initializing directory for task Gastrointestinal disorders
Initializing directory for task Social circumstances
Initializing directory for task Immune system disorders
Initializing directory for task Reproductive system and breast disorders
Initializing directory for task Neoplasms benign, malignant and unspecified (incl cysts and polyps)
Initializing directory for task General disorders and administration site conditions
Initializing directory for task Endocrine disorders
Initializing directory for task Surgical and medical procedures
Initializing directory for task Vascular disorders
Initializing d[...]<jupyter_text>Results<jupyter_code>#seed = 100, 99, 98, 97, 95

knn_fp_auc = [ 0.6489, 0.6077, 0.6403, 0.5954, 0.6458]
knn_lor_auc = [0.6752, 0.6724, 0.6710, 0.6691, 0.6919]
knn_euc_auc = [0.6815, 0.6887, 0.6857, 0.6804, 0.7213]
rf_fp_auc = [0.6826, 0.6684, 0.6665, 0.6783, 0.6497]
rf_lor_auc = [0.7123, 0.7227, 0.6968, 0.6928, 0.6983]
rf_euc_auc = [0.7099, 0.7018, 0.6825, 0.6906, 0.6905]

# no ATC information, only chemical structures
rf_lor_chem_auc = [0.6431, 0.6387, 0.6093, 0.5854, 0.6375]
rf_euc_chem_auc = [0.6359, 0.6140, 0.6292, 0.5955, 0.6215]
rf_fp_arr = np.array(rf_fp_auc)
rf_lor_arr = np.array(rf_lor_auc)
rf_euc_arr = np.array(rf_euc_auc)
print('FP mean:' + str(rf_fp_arr.mean()) + ', std:' + str(rf_fp_arr.std()))
print('Lorentz mean:' + str(rf_lor_arr.mean()) + ', std:' + str(rf_lor_arr.std()))
print('Euclidean mean:' + str(rf_euc_arr.mean()) + ', std:' + str(rf_euc_arr.std()))<jupyter_output>Euclidean mean:0.69506, std:0.009634230638717335<jupyter_text>Stats<jupyter_code>from scipy import stats<jupyter_output><empty_output><jupyter_text>Paired T-test between RF+FP and RF+LDE<jupyter_code>stats.ttest_rel(rf_fp_arr,rf_lor_arr)<jupyter_output><empty_output><jupyter_text>Paired T-test between RF+LDE and RF+EUC<jupyter_code>stats.ttest_rel(rf_euc_arr,rf_lor_arr)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11664890564800394, 'fraction_numerical': 0.030740224608574472, 'mean_word_length': 6.00746697300402, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 68, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '389362', 'license': 'Apache-2.0', 'n_tokens_mistral': 4938, 'n_tokens_neox': 4598, 'n_words': 883, 'path': 'DrugEmbedding/notebooks/SIDER.ipynb', 'repo_name': 'batmanlab/drugEmbedding', 'size': '12216'}