<jupyter_start><jupyter_text>MAHOMES holdout T-metal-site evaluationThis notebook trains an ML model using the metal ion site data-set and evaluates the model's predictions on the T-metal-site set. The models and scalers are saved during this.<jupyter_code># libraries
import numpy as np
import pandas as pd
import joblib

# scale features
from sklearn import preprocessing
from sklearn import impute
# classifier
from sklearn.ensemble import ExtraTreesClassifier
# scoring metrics
from sklearn.metrics import confusion_matrix, matthews_corrcoef

# custom scripts
import sys
sys.path.insert(0, "%s" % "CV/")
import GetFeatureSet as GetFeatureSet

# allow fancey printed strings 
# from https://stackoverflow.com/questions/8924173/how-do-i-print-bold-text-in-python/8930747
class color:
    PURPLE = '\033[95m'
    CYAN = '\033[96m'
    DARKCYAN = '\033[36m'
    BLUE = '\033[94m'
    GREEN = '\033[92m'
    YELLOW = '\033[93m'
    RED = '\033[91m'
    BOLD = '\033[1m'
    UNDERLINE = '\033[4m'
    END = '\033[0m'

print(color.BOLD + 'Hello World !' + color.END)
# Can be used to produce a new saved version of MAHOMES (ModelsForMAHOMES/) incase 
# the gitHub pickles don't work (generated using python 3.6.10 and MacOS 10.15.7)
save_models = False<jupyter_output><empty_output><jupyter_text>Read in data-set and T-metal-sites<jupyter_code>sites = pd.read_csv("publication_sites/sites_calculated_features.txt")
sites = sites.set_index('SITE_ID',drop=True)

# The following labels need to be changed after looking over literature (see Feehan, Franklin, Slusky 2021)
change_site_labels = ["5zb8_0", "6aci_0", "6oq7_0", "6pjv_1", "6q55_0",
                      "6q55_2", "6rmg_0", "6rtg_0", "6rw0_0", "6v77_0"]
# The following sites are removed due to unkopwn correct labels (see Feehan, Franklin, Slusky 2021)
sites.loc[sites.index.isin(change_site_labels), 'Catalytic']=True
remove_sites = ["6mf0_1", "6okh_0", "6qwo_0", "6r9n_0"]
sites=sites.loc[~sites.index.isin(remove_sites)]

print(color.BOLD + "All features:" + color.END)
print("sites: %s \tcolumns: %s"%(sites.shape[0], sites.shape[1]))
sizes = sites.groupby(["Set", "Catalytic"]).size()
print(sizes)<jupyter_output>[1mAll features:[0m
sites: 3981 	columns: 485
Set   Catalytic
data  False        2636
      True          829
test  False         345
      True          171
dtype: int64<jupyter_text>Normalize data<jupyter_code>def get_scaled_features():
    # seperate the sets (only dataset will be used to set scaling)
    data = sites.loc[sites.Set == "data"].copy()
    Tsites = sites.loc[sites.Set == "test"].copy()

    #split for scaling into categorical and not categorical
    not_ctg_geom = ("geom_gRMSD", "geom_MaxgRMSDDev","geom_val", "geom_nVESCUM","geom_AtomRMSD", "geom_AvgO", "geom_AvgN", "geom_AvgS", "geom_AvgOther", "geom_Charge")
    geom = [name for name in data if name.startswith("geom")]
    
    ctg_data = [x for x in geom if not x in not_ctg_geom]
    ctg_data.extend(["Set", 'Catalytic'])

    ## scale cont. features
    cont_scaler = preprocessing.RobustScaler(quantile_range=(20,80))
    #Fit scaler to X, then transform it
    data_nonctg = data[data.columns.difference(ctg_data)]#so that I can have columns
    data_scaled = pd.DataFrame(cont_scaler.fit_transform(data_nonctg), columns=data_nonctg.columns, index=data_nonctg.index)
    
    #scale the test set based on the scale of the training set
    Tsites_nonctg  = Tsites[Tsites.columns.difference(ctg_data)]
    Tsites_scaled = pd.DataFrame(cont_scaler.transform(Tsites_nonctg), columns=Tsites_nonctg.columns, index=Tsites_nonctg.index)
    
    #replace continuous feature null values with mean
    cont_imputer = impute.SimpleImputer(strategy="mean")
    data_scaled = pd.DataFrame(cont_imputer.fit_transform(data_scaled), columns=data_scaled.columns, index=data_scaled.index)
    Tsites_scaled = pd.DataFrame(cont_imputer.transform(Tsites_scaled), columns=Tsites_scaled.columns, index=Tsites_scaled.index)
    
    if save_models==True:
        joblib.dump(cont_scaler, "ModelsForMAHOMES/ContVarScaler.pkl")
        joblib.dump(cont_imputer, "ModelsForMAHOMES/ContVarImpute.pkl")
    
    #remove groupID so that it also isn't MinMax scaled either
    ctg_data.remove("Set");ctg_data.remove("Catalytic");
    #transform categorical data to [0,1] interval
    if len(data.columns.intersection(ctg_data)) > 0:
        ctg_scaler = preprocessing.MinMaxScaler()
        
        # fit the scaler to the data-set (training) and scale
        data_ctg = data[data.columns.intersection(ctg_data)]
        data_ctg_scaled = pd.DataFrame(ctg_scaler.fit_transform(data_ctg), columns=data_ctg.columns, index=data_ctg.index)
        
        #scale the test set based on the scale of the training set
        Tsites_ctg = Tsites[Tsites.columns.intersection(ctg_data)]
        Tsites_ctg_scaled = pd.DataFrame(ctg_scaler.transform(Tsites_ctg), columns=Tsites_ctg.columns, index=Tsites_ctg.index) 
        
        #replace categoric features null values with median value
        ctg_imputer = impute.SimpleImputer(strategy="median")
        data_ctg_scaled = pd.DataFrame(ctg_imputer.fit_transform(data_ctg_scaled), columns=data_ctg_scaled.columns, index=data_ctg_scaled.index)
        Tsites_ctg_scaled = pd.DataFrame(ctg_imputer.transform(Tsites_ctg_scaled), columns=Tsites_ctg_scaled.columns, index=Tsites_ctg_scaled.index)
        
        #concatenate the scaled categorical data to the robustly scaled data
        data_scaled = pd.merge(data_scaled, data_ctg_scaled, left_index=True, right_index=True)
        Tsites_scaled = pd.merge(Tsites_scaled, Tsites_ctg_scaled, left_index=True, right_index=True)
        
        if save_models==True:
            joblib.dump(ctg_scaler, "ModelsForMAHOMES/CtgVarScaler.pkl")
            joblib.dump(ctg_imputer, "ModelsForMAHOMES/CtgVarImpute.pkl")
            
    data_scaled = pd.merge(data_scaled, data['Catalytic'], left_index=True, right_index=True)
    Tsites_scaled = pd.merge(Tsites_scaled, Tsites['Catalytic'], left_index=True, right_index=True)
    
    return(data_scaled, Tsites_scaled)

data_scaled, Tsites_scaled = get_scaled_features()

print(color.BOLD + "All scaled data-set features:" + color.END)
print("sites: %s \tcolumns: %s"%(data_scaled.shape[0], data_scaled.shape[1]))
print(data_scaled.groupby(["Catalytic"]).size())

print(color.BOLD + "\nAll scaled T-metal-site features:" + color.END)
print("sites: %s \tcolumns: %s"%(Tsites_scaled.shape[0], Tsites_scaled.shape[1]))
print(Tsites_scaled.groupby(["Catalytic"]).size())<jupyter_output>[1mAll scaled data-set features:[0m
sites: 3465 	columns: 484
Catalytic
False    2636
True      829
dtype: int64
[1m
All scaled T-metal-site features:[0m
sites: 516 	columns: 484
Catalytic
False    345
True     171
dtype: int64<jupyter_text>Evaluate performace using T-metal-sites<jupyter_code>## returns relevent data-set data for training ML model
def get_training_data(feature_set, random_seed):
    ## random under sample data-set (1+:3-)
    X_Cat = data_scaled[data_scaled['Catalytic']==True].copy()
    X_nonCat = data_scaled[data_scaled['Catalytic']==False].copy()
    X_nonCat = X_nonCat.sample(n=len(X_Cat)*3, axis=0, random_state=random_seed)
    X_prep = X_Cat.append(X_nonCat)
    
    ## seperate target value
    y = X_prep['Catalytic']; del X_prep['Catalytic']
    
    ## only return features in specific feature set
    X = GetFeatureSet.feature_subset(X_prep, feature_set, noBSA=True)
    
    return(X, y)

## number of iterations to improve reproducability
num_rand_seeds = 10 # 10 provides 3 decimal level reproducability across my machines)
def evaluate_model_with_Tsite(clf, feature_set):
    ## prepare test-set
    testX = Tsites_scaled.copy(); 
    testY = testX['Catalytic']; del testX['Catalytic']
    testX = GetFeatureSet.feature_subset(testX, feature_set, noBSA=True)
    
    ## get multiple predictions for test-set w/ diff random seeds
    test_site_preds = {'actual': pd.Series(testY, index=testX.index)}
    for rand_seed in range(0,num_rand_seeds):
        # get undersampled training data for feature set 
        X, y = get_training_data(feature_set, rand_seed)
        print("random_seed = %s"%(rand_seed), end="\t")
        print("(num. training sites= %s (%s+ : %s-) \tnum. features: %s)"%(X.shape[0], len(y[y==True]),len(y[y==False]), X.shape[1]))
        
        ## train classifier and make test-set predictions
        clf.set_params(random_state=rand_seed)
        clf.fit(X, y)
        test_preds = clf.predict(testX)
        test_site_preds['prediction_%s'%(rand_seed)]= pd.Series(test_preds, index=testX.index)
        if save_models==True:
            joblib.dump(clf, "ModelsForMAHOMES/MAHOMES%s.pkl"%(rand_seed))
        
        ## output results for this random seed to get an idea of prediction variation levels
        TN, FP, FN, TP = confusion_matrix(testY, test_preds).ravel()
        mcc = matthews_corrcoef(testY, test_preds)
        print("\tTP=%s \tTN=%s \tFP=%s \tFN=%s"%(TP, TN, FP, FN))

    ## calcualte the average of all random seed predictions
    test_predictions = pd.DataFrame(test_site_preds)
    test_predictions['prediction']=0
    for rand_seed in range(0,num_rand_seeds):
        test_predictions['prediction']+=test_predictions['prediction_%s'%(rand_seed)] 
    test_predictions['prediction']=test_predictions['prediction']/num_rand_seeds
    
    ## make final prediction
    test_predictions['bool_pred']=False
    test_predictions.loc[test_predictions['prediction']>=0.5, 'bool_pred']=True
    
    return(test_predictions)

## return result metrics for final predictions
def check_result_metrics(alg, feat_set, prediction_df):
    mcc = matthews_corrcoef(prediction_df['actual'], prediction_df['bool_pred'])
    TN, FP, FN, TP = confusion_matrix(prediction_df['actual'], prediction_df['bool_pred']).ravel()
    
    TPR=(TP/(TP+FN))*100
    TNR=(TN/(TN+FP))*100
    acc=((TP+TN)/(TP+TN+FP+FN))*100
    Prec=(TP/(TP+FP))*100
    return(pd.DataFrame([[alg, feat_set, acc, mcc, TPR, TNR, Prec]],
        columns=['Algorithm', 'Feature Set', 'Accuracy', 'MCC', 'Recall', 'TrueNegRate', 'Precision']))
# specify algorithm and feature set for MAHOMES (top model from outer CV results)
MAHOMES_alg = "ExtraTrees"
MAHOMES_feature_set = "AllMeanSph"

## set extra trees classifier with optimal parameters found during inner CV
MAHOMES_clf = ExtraTreesClassifier(n_estimators=500, min_samples_split=3, max_depth=None,
                           criterion="gini",bootstrap=False, class_weight=None,
                           max_features=None)


MAHOMES_Tsite_predictions = evaluate_model_with_Tsite(MAHOMES_clf, MAHOMES_feature_set)

## display final results
scores = check_result_metrics(MAHOMES_alg, MAHOMES_feature_set,  MAHOMES_Tsite_predictions)
display(scores.round(2))<jupyter_output>random_seed = 0	(num. training sites= 3316 (829+ : 2487-) 	num. features: 181)
	TP=153 	TN=333 	FP=12 	FN=18
random_seed = 1	(num. training sites= 3316 (829+ : 2487-) 	num. features: 181)
	TP=150 	TN=333 	FP=12 	FN=21
random_seed = 2	(num. training sites= 3316 (829+ : 2487-) 	num. features: 181)
	TP=154 	TN=329 	FP=16 	FN=17
random_seed = 3	(num. training sites= 3316 (829+ : 2487-) 	num. features: 181)
	TP=153 	TN=332 	FP=13 	FN=18
random_seed = 4	(num. training sites= 3316 (829+ : 2487-) 	num. features: 181)
	TP=154 	TN=331 	FP=14 	FN=17
random_seed = 5	(num. training sites= 3316 (829+ : 2487-) 	num. features: 181)
	TP=153 	TN=330 	FP=15 	FN=18
random_seed = 6	(num. training sites= 3316 (829+ : 2487-) 	num. features: 181)
	TP=153 	TN=331 	FP=14 	FN=18
random_seed = 7	(num. training sites= 3316 (829+ : 2487-) 	num. features: 181)
	TP=154 	TN=331 	FP=14 	FN=17
random_seed = 8	(num. training sites= 3316 (829+ : 2487-) 	num. features: 181)
	TP=151 	TN=333 	FP=12 	FN=20
random_seed = 9	(nu[...]
{'fraction_non_alphanumeric': 0.12010819034739244, 'fraction_numerical': 0.03786662158735525, 'mean_word_length': 4.298701298701299, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 16, 'https://': 1, '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': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '229068', 'license': 'MIT', 'n_tokens_mistral': 4687, 'n_tokens_neox': 4266, 'n_words': 928, 'path': 'MAHOMES_eval_T-metal-site.ipynb', 'repo_name': 'SluskyLab/MAHOMES', 'size': '11831'}