scripts/UniKP_Km.py

import _bootstrap
import torch
from model.build_vocab import WordVocab
from model.pretrain_trfm import TrfmSeq2seq
from model.utils import split
import json
from transformers import T5EncoderModel, T5Tokenizer
import re
import gc
from sklearn import metrics
from sklearn.ensemble import ExtraTreesRegressor
import numpy as np
import pandas as pd
from sklearn.model_selection import KFold
from sklearn.metrics import mean_squared_error, mean_absolute_error
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split
import random
import pickle
import math
from sklearn.model_selection import train_test_split
import xgboost
import lightgbm
from sklearn.model_selection import GridSearchCV
from project_paths import (

    DEGREE_SMILES_PATH, KCAT_DATASET_PATH, KCAT_KM_SAMPLES_PATH, KM_TEST_PATH,

    PH_SMILES_PATH, PROT_T5_MODEL, TRFM_PATH, UNIKP_MODEL_DIR, VOCAB_PATH,

)



def smiles_to_vec(Smiles):
    pad_index = 0
    unk_index = 1
    eos_index = 2
    sos_index = 3
    mask_index = 4
    vocab = WordVocab.load_vocab(VOCAB_PATH)
    def get_inputs(sm):
        seq_len = 220
        sm = sm.split()
        if len(sm)>218:
            print('SMILES is too long ({:d})'.format(len(sm)))
            sm = sm[:109]+sm[-109:]
        ids = [vocab.stoi.get(token, unk_index) for token in sm]
        ids = [sos_index] + ids + [eos_index]
        seg = [1]*len(ids)
        padding = [pad_index]*(seq_len - len(ids))
        ids.extend(padding), seg.extend(padding)
        return ids, seg
    def get_array(smiles):
        x_id, x_seg = [], []
        for sm in smiles:
            a,b = get_inputs(sm)
            x_id.append(a)
            x_seg.append(b)
        return torch.tensor(x_id), torch.tensor(x_seg)
    trfm = TrfmSeq2seq(len(vocab), 256, len(vocab), 4)
    trfm.load_state_dict(torch.load(TRFM_PATH))
    trfm.eval()
    x_split = [split(sm) for sm in Smiles]
    xid, xseg = get_array(x_split)
    X = trfm.encode(torch.t(xid))
    return X


def Seq_to_vec(Sequence):
    for i in range(len(Sequence)):
        if len(Sequence[i]) > 1000:
            Sequence[i] = Sequence[i][:500] + Sequence[i][-500:]
    sequences_Example = []
    for i in range(len(Sequence)):
        zj = ''
        for j in range(len(Sequence[i]) - 1):
            zj += Sequence[i][j] + ' '
        zj += Sequence[i][-1]
        sequences_Example.append(zj)
    tokenizer = T5Tokenizer.from_pretrained(PROT_T5_MODEL, do_lower_case=False)
    model = T5EncoderModel.from_pretrained(PROT_T5_MODEL)
    gc.collect()
    print(torch.cuda.is_available())
    # 'cuda:0' if torch.cuda.is_available() else
    device = torch.device('cuda:1' if torch.cuda.is_available() else 'cpu')
    model = model.to(device)
    model = model.eval()
    features = []
    for i in range(len(sequences_Example)):
        print('For sequence ', str(i+1))
        sequences_Example_i = sequences_Example[i]
        sequences_Example_i = [re.sub(r"[UZOB]", "X", sequences_Example_i)]
        ids = tokenizer.batch_encode_plus(sequences_Example_i, add_special_tokens=True, padding=True)
        input_ids = torch.tensor(ids['input_ids']).to(device)
        attention_mask = torch.tensor(ids['attention_mask']).to(device)
        with torch.no_grad():
            embedding = model(input_ids=input_ids, attention_mask=attention_mask)
        embedding = embedding.last_hidden_state.cpu().numpy()
        for seq_num in range(len(embedding)):
            seq_len = (attention_mask[seq_num] == 1).sum()
            seq_emd = embedding[seq_num][:seq_len - 1]
            features.append(seq_emd)
    features_normalize = np.zeros([len(features), len(features[0][0])], dtype=float)
    for i in range(len(features)):
        for k in range(len(features[0][0])):
            for j in range(len(features[i])):
                features_normalize[i][k] += features[i][j][k]
            features_normalize[i][k] /= len(features[i])
    return features_normalize


def Kcat_predict(Ifeature_ini, Label_ini):
    Myseed = random.randint(0, 1000)
    print(Myseed)
    Ifeature, Ifeature_test, Label, Label_test = train_test_split(Ifeature_ini, Label_ini, test_size=0.2,
                                                                  random_state=Myseed)
    model = ExtraTreesRegressor()
    model.fit(Ifeature, Label)
    Pre_label = model.predict(Ifeature_test)
    Pcc = np.corrcoef(Label_test, Pre_label)[1][0]
    RMSE = np.sqrt(mean_squared_error(Label_test, Pre_label))
    MAE = mean_absolute_error(Label_test, Pre_label)
    r2 = r2_score(Label_test, Pre_label)
    print(r2, Pcc, RMSE, MAE)


if __name__ == '__main__':
    # Dataset Load
    with open(KM_TEST_PATH, 'rb') as file:
        datasets = pickle.load(file)
    # datasets = datasets
    # print(datasets)
    sequence = datasets['Sequence']
    smiles = datasets['smiles']
    Label = datasets['log10_KM']
    print(len(smiles), len(Label))
    # smiles_input = smiles_to_vec(smiles)
    # sequence_input = Seq_to_vec(sequence)
    # feature = np.concatenate((smiles_input, sequence_input), axis=1)
    # with open("Km/Km_features_11722_PreKcat.pkl", "wb") as f:
    #     pickle.dump(feature, f)
    with open("Km/Km_features_11722_PreKcat.pkl", "rb") as f:
        feature = pickle.load(f)
    feature = np.array(feature)
    Label = np.array(Label)
    Kcat_predict(feature, Label)