initial adding of ml scripts for baseline models

scripts/ml/MultModelsCl.py

@@ -0,0 +1,272 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Mar  4 15:25:33 2022
+
+@author: tanu
+"""
+#%%
+import os, sys
+import pandas as pd
+import numpy as np
+import pprint as pp
+from copy import deepcopy
+from sklearn import linear_model
+from sklearn import datasets
+from collections import Counter
+
+from sklearn.linear_model import LogisticRegression, LogisticRegressionCV
+from sklearn.linear_model import RidgeClassifier, RidgeClassifierCV, SGDClassifier, PassiveAggressiveClassifier
+
+from sklearn.naive_bayes import BernoulliNB
+from sklearn.neighbors import KNeighborsClassifier
+from sklearn.svm import SVC
+from sklearn.tree import DecisionTreeClassifier, ExtraTreeClassifier
+from sklearn.ensemble import RandomForestClassifier, ExtraTreesClassifier, AdaBoostClassifier, GradientBoostingClassifier, BaggingClassifier
+from sklearn.naive_bayes import GaussianNB
+from sklearn.gaussian_process import GaussianProcessClassifier, kernels
+from sklearn.gaussian_process.kernels import RBF, DotProduct, Matern, RationalQuadratic, WhiteKernel
+
+from sklearn.discriminant_analysis import LinearDiscriminantAnalysis, QuadraticDiscriminantAnalysis
+from sklearn.neural_network import MLPClassifier
+
+from sklearn.svm import SVC
+from xgboost import XGBClassifier
+from sklearn.naive_bayes import MultinomialNB
+from sklearn.preprocessing import StandardScaler, MinMaxScaler, OneHotEncoder
+
+from sklearn.compose import ColumnTransformer
+from sklearn.compose import make_column_transformer
+
+from sklearn.metrics import make_scorer, confusion_matrix, accuracy_score, balanced_accuracy_score, precision_score, average_precision_score, recall_score
+from sklearn.metrics import roc_auc_score, roc_curve, f1_score, matthews_corrcoef, jaccard_score, classification_report
+
+from sklearn.model_selection import train_test_split, cross_validate, cross_val_score
+from sklearn.model_selection import StratifiedKFold,RepeatedStratifiedKFold, RepeatedKFold
+
+from sklearn.pipeline import Pipeline, make_pipeline
+
+from sklearn.feature_selection import RFE, RFECV
+
+import itertools
+import seaborn as sns
+import matplotlib.pyplot as plt
+
+from statistics import mean, stdev, median, mode
+
+from imblearn.over_sampling import RandomOverSampler
+from imblearn.under_sampling import RandomUnderSampler
+from imblearn.over_sampling import SMOTE
+from sklearn.datasets import make_classification
+from imblearn.combine import SMOTEENN
+from imblearn.combine import SMOTETomek
+
+from imblearn.over_sampling import SMOTENC
+from imblearn.under_sampling import EditedNearestNeighbours
+from imblearn.under_sampling import RepeatedEditedNearestNeighbours
+
+from sklearn.model_selection import GridSearchCV
+from sklearn.base import BaseEstimator
+from sklearn.impute import KNNImputer as KNN
+import json
+
+#%% GLOBALS
+rs = {'random_state': 42}
+njobs = {'n_jobs': 10}
+
+scoring_fn =  ({ 'mcc'        : make_scorer(matthews_corrcoef)
+                , 'accuracy'      : make_scorer(accuracy_score)
+                , 'fscore'     : make_scorer(f1_score)
+                , 'precision'  : make_scorer(precision_score)
+                , 'recall'     : make_scorer(recall_score)
+                , 'roc_auc'    : make_scorer(roc_auc_score)
+                , 'jcc'        : make_scorer(jaccard_score)
+            }) 
+  
+skf_cv = StratifiedKFold(n_splits = 10
+                          #, shuffle = False, random_state= None)
+                           , shuffle = True,**rs)
+
+rskf_cv = RepeatedStratifiedKFold(n_splits = 10
+                                  , n_repeats = 3
+                                  , **rs)
+
+mcc_score_fn = {'mcc': make_scorer(matthews_corrcoef)}
+jacc_score_fn = {'jcc': make_scorer(jaccard_score)}
+#%%
+# Multiple Classification - Model Pipeline
+def MultModelsCl(input_df, target, skf_cv
+                       , blind_test_input_df
+                       , blind_test_target
+                       , var_type = ['numerical', 'categorical','mixed']):
+
+    '''
+    @ param input_df: input features 
+    @ type: df with input features WITHOUT the target variable
+    
+    @param target: target (or output) feature
+    @type: df or np.array or Series
+    
+    @param skv_cv: stratifiedK fold int or object to allow shuffle and random state to pass
+    @type: int or StratifiedKfold()
+    
+    @var_type: numerical, categorical and mixed to determine what col_transform to apply (MinMaxScalar and/or one-ho    t encoder)
+    @type: list
+
+    returns
+    Dict containing multiple classification scores for each model and mean of each Stratified Kfold including training
+    '''
+    
+    # Determine categorical and numerical features
+    numerical_ix = input_df.select_dtypes(include=['int64', 'float64']).columns
+    numerical_ix
+    categorical_ix = input_df.select_dtypes(include=['object', 'bool']).columns
+    categorical_ix    
+
+    # Determine preprocessing steps ~ var_type
+    if var_type == 'numerical':
+        t = [('num', MinMaxScaler(), numerical_ix)]
+
+    if var_type == 'categorical':
+        t = [('cat', OneHotEncoder(), categorical_ix)]
+    
+    if var_type == 'mixed':
+        t = [('num', MinMaxScaler(), numerical_ix)
+            , ('cat', OneHotEncoder(), categorical_ix) ]
+        
+    col_transform = ColumnTransformer(transformers = t
+                                       , remainder='passthrough')
+    
+    # Specify multiple Classification models
+    lr      = LogisticRegression(**rs)
+    lrcv    = LogisticRegressionCV(**rs)
+    gnb     = GaussianNB()
+    nb      = BernoulliNB()
+    knn     = KNeighborsClassifier()
+    svc     = SVC(**rs)
+    mlp     = MLPClassifier(max_iter = 500, **rs)
+    dt      = DecisionTreeClassifier(**rs)
+    ets     = ExtraTreesClassifier(**rs)
+    
+    rf      = RandomForestClassifier(**rs, n_estimators = 1000 )
+    rf2     = RandomForestClassifier(
+                          min_samples_leaf = 5
+                          , n_estimators     = 1000
+                          , bootstrap        = True
+                          , oob_score        = True
+                          , **njobs
+                          , **rs
+                          , max_features     = 'auto')
+    xgb = XGBClassifier(**rs, verbosity = 0, use_label_encoder =False)
+                        
+    lda = LinearDiscriminantAnalysis()
+    
+    mnb = MultinomialNB()
+    
+    pa  = PassiveAggressiveClassifier(**rs, **njobs)
+    
+    sgd = SGDClassifier(**rs, **njobs)  
+
+    abc = AdaBoostClassifier(**rs)
+    bc  = BaggingClassifier(**rs, **njobs, bootstrap = True, oob_score = True)
+    et  = ExtraTreeClassifier(**rs)
+    gpc = GaussianProcessClassifier(**rs)
+    gbc = GradientBoostingClassifier(**rs)
+    qda = QuadraticDiscriminantAnalysis()
+    rc  = RidgeClassifier(**rs)
+    rccv  = RidgeClassifierCV(cv = 10)
+    
+    models = [('Logistic Regression'       , lr)
+            , ('Logistic RegressionCV'     , lrcv)
+            , ('Gaussian NB'               , gnb)
+            , ('Naive Bayes'               , nb)
+            , ('K-Nearest Neighbors'       , knn) 
+            , ('SVM'                       , svc) 
+            , ('MLP'                       , mlp) 
+            , ('Decision Tree'             , dt) 
+            , ('Extra Trees'               , ets) 
+            , ('Extra Tree'                , et)
+            , ('Random Forest'             , rf) 
+            , ('Random Forest2'            , rf2) 
+            , ('Naive Bayes'               , nb)
+            , ('XGBoost'                   , xgb)
+            , ('LDA'                       , lda)
+            , ('Multinomial'               , mnb)
+            , ('Passive Aggresive'         , pa)
+            , ('Stochastic GDescent'       , sgd)
+            , ('AdaBoost Classifier'       , abc)
+            , ('Bagging Classifier'        , bc)
+            , ('Gaussian Process'          , gpc)
+            , ('Gradient Boosting'         , gbc)
+            , ('QDA'                       , qda)
+            , ('Ridge Classifier'          , rc)
+            , ('Ridge ClassifierCV'        , rccv)
+            ]
+        
+    mm_skf_scoresD = {}
+     
+    for model_name, model_fn in models:
+        print('\nModel_name:', model_name
+        , '\nModel func:'    , model_fn
+        , '\nList of models:', models)
+    
+        model_pipeline = Pipeline([
+            ('prep'     , col_transform)
+            , ('model'  , model_fn)])
+            
+        print('Running model pipeline:', model_pipeline)
+        skf_cv_mod = cross_validate(model_pipeline
+                              , input_df
+                              , target
+                              , cv = skf_cv
+                              , scoring = scoring_fn
+                              , return_train_score = True)
+        mm_skf_scoresD[model_name] = {}
+        for key, value in skf_cv_mod.items():
+            print('\nkey:', key, '\nvalue:', value)
+            print('\nmean value:', mean(value))
+            mm_skf_scoresD[model_name][key] = round(mean(value),2)
+            #pp.pprint(mm_skf_scoresD)
+            #cvtrain_mcc = mm_skf_scoresD[model_name]['test_mcc']
+            
+    #return(mm_skf_scoresD)
+#%%
+        #=========================
+        # Blind test: BTS results
+        #=========================
+        # Build the final results with all scores for a feature selected model
+        #bts_predict = gscv_fs.predict(blind_test_input_df)
+        model_pipeline.fit(input_df, target)
+        bts_predict = model_pipeline.predict(blind_test_input_df)
+        
+        bts_mcc_score = round(matthews_corrcoef(blind_test_target, bts_predict),2)
+        print('\nMCC on Blind test:'     , bts_mcc_score)
+        print('\nAccuracy on Blind test:', round(accuracy_score(blind_test_target, bts_predict),2))
+        
+        # Diff b/w train and bts test scores
+        #train_test_diff_MCC = cvtrain_mcc - bts_mcc_score
+        # print('\nDiff b/w train and blind test score (MCC):', train_test_diff)
+        
+       
+        # # create a dict with all scores
+        # lr_btsD = { 'model_name': model_name
+        #            , 'bts_mcc':None
+        #                , 'bts_fscore':None
+        #                , 'bts_precision':None
+        #                , 'bts_recall':None
+        #                , 'bts_accuracy':None
+        #                , 'bts_roc_auc':None
+        #                , 'bts_jaccard':None}
+        
+        
+        mm_skf_scoresD[model_name]['bts_mcc']       = bts_mcc_score
+        mm_skf_scoresD[model_name]['bts_fscore']    = round(f1_score(blind_test_target, bts_predict),2)
+        mm_skf_scoresD[model_name]['bts_precision'] = round(precision_score(blind_test_target, bts_predict),2)
+        mm_skf_scoresD[model_name]['bts_recall']    = round(recall_score(blind_test_target, bts_predict),2)
+        mm_skf_scoresD[model_name]['bts_accuracy']  = round(accuracy_score(blind_test_target, bts_predict),2)
+        mm_skf_scoresD[model_name]['bts_roc_auc']   = round(roc_auc_score(blind_test_target, bts_predict),2)
+        mm_skf_scoresD[model_name]['bts_jaccard']   = round(jaccard_score(blind_test_target, bts_predict),2)
+        #mm_skf_scoresD[model_name]['diff_mcc']      = train_test_diff_MCC
+
+    return(mm_skf_scoresD)
+    

scripts/ml/ml_data.py

@@ -0,0 +1,537 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Mar  6 13:41:54 2022
+
+@author: tanu
+"""
+def setvars(gene,drug):
+    #https://stackoverflow.com/questions/51695322/compare-multiple-algorithms-with-sklearn-pipeline
+    import os, sys
+    import pandas as pd
+    import numpy as np
+    print(np.__version__)
+    print(pd.__version__)
+    import pprint as pp
+    from copy import deepcopy
+    from collections import Counter
+    from sklearn.impute import KNNImputer as KNN
+    from imblearn.over_sampling import RandomOverSampler
+    from imblearn.under_sampling import RandomUnderSampler
+    from imblearn.over_sampling import SMOTE
+    from sklearn.datasets import make_classification
+    from imblearn.combine import SMOTEENN
+    from imblearn.combine import SMOTETomek
+    
+    from imblearn.over_sampling import SMOTENC
+    from imblearn.under_sampling import EditedNearestNeighbours
+    from imblearn.under_sampling import RepeatedEditedNearestNeighbours
+    
+    from sklearn.metrics import make_scorer, confusion_matrix, accuracy_score, balanced_accuracy_score, precision_score, average_precision_score, recall_score
+    from sklearn.metrics import roc_auc_score, roc_curve, f1_score, matthews_corrcoef, jaccard_score, classification_report
+
+    from sklearn.model_selection import train_test_split, cross_validate, cross_val_score
+    from sklearn.model_selection import StratifiedKFold,RepeatedStratifiedKFold, RepeatedKFold
+
+    from sklearn.pipeline import Pipeline, make_pipeline
+    #%% GLOBALS
+    rs = {'random_state': 42}
+    njobs = {'n_jobs': 10}
+
+    scoring_fn =  ({ 'mcc'         : make_scorer(matthews_corrcoef)
+                    , 'accuracy'   : make_scorer(accuracy_score)
+                    , 'fscore'     : make_scorer(f1_score)
+                    , 'precision'  : make_scorer(precision_score)
+                    , 'recall'     : make_scorer(recall_score)
+                    , 'roc_auc'    : make_scorer(roc_auc_score)
+                    , 'jcc'        : make_scorer(jaccard_score)
+                }) 
+      
+    skf_cv = StratifiedKFold(n_splits = 10
+                              #, shuffle = False, random_state= None)
+                               , shuffle = True,**rs)
+
+    rskf_cv = RepeatedStratifiedKFold(n_splits = 10
+                                      , n_repeats = 3
+                                      , **rs)
+
+    mcc_score_fn  = {'mcc': make_scorer(matthews_corrcoef)}
+    jacc_score_fn = {'jcc': make_scorer(jaccard_score)}
+
+    #%% FOR LATER: Combine ED logo data
+    #%% FOR LATER: active aa site annotations **DONE on 15/05/2022 as part of generating merged_dfs
+    ###########################################################################
+    rs = {'random_state': 42}
+    njobs = {'n_jobs': 10}
+    homedir = os.path.expanduser("~")
+    
+    #==============
+    # directories
+    #==============
+    datadir = homedir + '/git/Data/'
+    indir   = datadir + drug + '/input/'
+    outdir  = datadir + drug + '/output/'
+    
+    #=======
+    # input
+    #=======
+    infile_ml1 = outdir + gene.lower() + '_merged_df3.csv' 
+    #infile_ml2 = outdir + gene.lower() + '_merged_df2.csv'
+    
+    my_df = pd.read_csv(infile_ml1, index_col = 0)
+    my_df.dtypes
+    my_df_cols = my_df.columns
+    
+    geneL_basic     = ['pnca']
+    geneL_na        = ['gid']
+    geneL_na_ppi2   = ['rpob']
+    geneL_ppi2      = ['alr', 'embb', 'katg']
+    #%% get cols
+    mycols = my_df.columns
+    
+    # # change from numberic to 
+    # num_type = ['int64', 'float64']
+    # cat_type = ['object', 'bool']
+    
+    # if my_df['active_aa_pos'].dtype in num_type:
+    #     my_df['active_aa_pos'] = my_df['active_aa_pos'].astype(object)
+    #     my_df['active_aa_pos'].dtype
+    
+    # FIXME: if this is not structural, remove from source..
+    # Drop NA where numerical cols have them
+    if gene.lower() in geneL_na_ppi2:
+        #D1148 get rid of
+        na_index = my_df['mutationinformation'].index[my_df['mcsm_na_affinity'].apply(np.isnan)]
+        my_df = my_df.drop(index=na_index)
+    
+    # FIXME: either impute or remove!
+    # for embb (L114M, F115L, V123L, V125I, V131M) delete for now
+    if gene.lower() in ['embb']:
+        na_index = my_df['mutationinformation'].index[my_df['ligand_distance'].apply(np.isnan)]
+        my_df = my_df.drop(index=na_index)# RERUN embb with the 5 values now present
+       
+    ###########################################################################
+    #%% Add lineage calculation columns
+    #FIXME: Check if this can be imported from config?
+    total_mtblineage_uc = 8
+    lineage_colnames = ['lineage_list_all', 'lineage_count_all', 'lineage_count_unique', 'lineage_list_unique', 'lineage_multimode']
+    #bar = my_df[lineage_colnames]
+    my_df['lineage_proportion']      = my_df['lineage_count_unique']/my_df['lineage_count_all']
+    my_df['dist_lineage_proportion'] = my_df['lineage_count_unique']/total_mtblineage_uc
+    ###########################################################################
+    #%% AA property change
+    #--------------------
+    # Water prop change
+    #--------------------
+    my_df['water_change'] = my_df['wt_prop_water'] + str('_to_') + my_df['mut_prop_water']
+    my_df['water_change'].value_counts()
+    
+    water_prop_changeD = {
+        'hydrophobic_to_neutral'          : 'change'
+        , 'hydrophobic_to_hydrophobic'    : 'no_change'
+        , 'neutral_to_neutral'            : 'no_change'
+        , 'neutral_to_hydrophobic'        : 'change'
+        , 'hydrophobic_to_hydrophilic'    : 'change'
+        , 'neutral_to_hydrophilic'        : 'change'
+        , 'hydrophilic_to_neutral'        : 'change'
+        , 'hydrophilic_to_hydrophobic'    : 'change'
+        , 'hydrophilic_to_hydrophilic'    : 'no_change'
+    }
+
+    my_df['water_change'] = my_df['water_change'].map(water_prop_changeD)
+    my_df['water_change'].value_counts()
+    
+    #--------------------
+    # Polarity change
+    #--------------------
+    my_df['polarity_change'] = my_df['wt_prop_polarity'] + str('_to_') + my_df['mut_prop_polarity']
+    my_df['polarity_change'].value_counts()
+
+    polarity_prop_changeD = {
+        'non-polar_to_non-polar'     : 'no_change'
+        , 'non-polar_to_neutral'     : 'change'  
+        , 'neutral_to_non-polar'     : 'change'  
+        , 'neutral_to_neutral'       : 'no_change'  
+        , 'non-polar_to_basic'       : 'change'  
+        , 'acidic_to_neutral'        : 'change'  
+        , 'basic_to_neutral'         : 'change'  
+        , 'non-polar_to_acidic'      : 'change'  
+        , 'neutral_to_basic'         : 'change'  
+        , 'acidic_to_non-polar'      : 'change'  
+        , 'basic_to_non-polar'       : 'change'
+        , 'neutral_to_acidic'        : 'change'
+        , 'acidic_to_acidic'         : 'no_change'
+        , 'basic_to_acidic'          : 'change'
+        , 'basic_to_basic'           : 'no_change'
+        , 'acidic_to_basic'          : 'change'}
+
+    my_df['polarity_change'] = my_df['polarity_change'].map(polarity_prop_changeD)
+    my_df['polarity_change'].value_counts()
+    
+    #--------------------
+    # Electrostatics change
+    #--------------------
+    my_df['electrostatics_change'] = my_df['wt_calcprop'] + str('_to_') + my_df['mut_calcprop']
+    my_df['electrostatics_change'].value_counts()
+
+    calc_prop_changeD = {
+            'non-polar_to_non-polar'     : 'no_change'
+            , 'non-polar_to_polar'       : 'change'
+            , 'polar_to_non-polar'       : 'change'
+            , 'non-polar_to_pos'         : 'change'
+            , 'neg_to_non-polar'         : 'change'
+            , 'non-polar_to_neg'         : 'change'
+            , 'pos_to_polar'             : 'change'
+            , 'pos_to_non-polar'         : 'change'
+            , 'polar_to_polar'           : 'no_change'
+            , 'neg_to_neg'               : 'no_change'
+            , 'polar_to_neg'             : 'change'
+            , 'pos_to_neg'               : 'change'
+            , 'pos_to_pos'               : 'no_change'
+            , 'polar_to_pos'             : 'change'
+            , 'neg_to_polar'             : 'change'
+            , 'neg_to_pos'               : 'change'
+    }
+
+    my_df['electrostatics_change'] = my_df['electrostatics_change'].map(calc_prop_changeD)
+    my_df['electrostatics_change'].value_counts()
+
+    #--------------------    
+    # Summary change: Create a combined column summarising these three cols
+    #--------------------
+    detect_change = 'change'
+    check_prop_cols = ['water_change', 'polarity_change', 'electrostatics_change']
+    #my_df['aa_prop_change'] = (my_df.values == detect_change).any(1).astype(int)
+    my_df['aa_prop_change'] = (my_df[check_prop_cols].values == detect_change).any(1).astype(int)
+    my_df['aa_prop_change'].value_counts()
+    my_df['aa_prop_change'].dtype
+    
+    my_df['aa_prop_change'] = my_df['aa_prop_change'].map({1:'change'
+                                                           , 0: 'no_change'})
+
+    my_df['aa_prop_change'].value_counts()
+    my_df['aa_prop_change'].dtype
+    
+    #%% IMPUTE values for OR [check script for exploration: UQ_or_imputer]
+    #or_cols = ['or_mychisq', 'log10_or_mychisq', 'or_fisher']
+    sel_cols = ['mutationinformation', 'or_mychisq', 'log10_or_mychisq']
+    or_cols = ['or_mychisq', 'log10_or_mychisq']
+    
+    print("count of NULL values before imputation\n")
+    print(my_df[or_cols].isnull().sum())
+    
+    my_dfI = pd.DataFrame(index = my_df['mutationinformation'] )
+    
+        
+    my_dfI = pd.DataFrame(KNN(n_neighbors= 5, weights="uniform").fit_transform(my_df[or_cols])
+                          , index =  my_df['mutationinformation']
+                          , columns = or_cols )
+    my_dfI.columns = ['or_rawI', 'logorI']
+    my_dfI.columns
+    my_dfI = my_dfI.reset_index(drop = False) # prevents old index from being added as a column
+    my_dfI.head()
+    print("count of NULL values AFTER imputation\n")
+    print(my_dfI.isnull().sum())
+    
+    #-------------------------------------------
+    # OR df Merge: with original based on index
+    #-------------------------------------------
+    my_df['index_bm'] = my_df.index
+    mydf_imputed = pd.merge(my_df
+                        , my_dfI
+                        , on = 'mutationinformation')
+    mydf_imputed = mydf_imputed.set_index(['index_bm'])
+    
+    my_df['log10_or_mychisq'].isna().sum()
+    mydf_imputed['log10_or_mychisq'].isna().sum()
+    mydf_imputed['logorI'].isna().sum()
+    
+    len(my_df.columns)
+    len(mydf_imputed.columns)  
+    
+    #-----------------------------------------
+    # REASSIGN my_df after imputing OR values
+    #-----------------------------------------
+    my_df = mydf_imputed.copy()
+
+    #%%########################################################################
+    #==========================
+    #     Data for ML
+    #==========================
+    my_df_ml = my_df.copy()
+    
+    #==========================
+    #     BLIND test set
+    #==========================
+    # Separate blind test set
+    my_df_ml[drug].isna().sum()
+    
+    blind_test_df = my_df_ml[my_df_ml[drug].isna()]
+    blind_test_df.shape
+    
+    training_df =  my_df_ml[my_df_ml[drug].notna()]
+    training_df.shape
+    
+    # Target1: dst
+    training_df[drug].value_counts()
+    training_df['dst_mode'].value_counts()
+    
+    #%% Build X: input for ML
+    common_cols_stabiltyN = ['ligand_distance'
+               , 'ligand_affinity_change'
+               , 'duet_stability_change'
+               , 'ddg_foldx'
+               , 'deepddg'
+               , 'ddg_dynamut2'
+               , 'mmcsm_lig'
+               , 'contacts']
+    
+    # Build stability columns ~ gene
+    if gene.lower() in geneL_basic:
+        X_stabilityN = common_cols_stabiltyN
+        cols_to_mask = ['ligand_affinity_change']
+        
+    if gene.lower() in geneL_ppi2:
+    #    X_stabilityN = common_cols_stabiltyN + ['mcsm_ppi2_affinity' , 'interface_dist'] 
+        geneL_ppi2_st_cols = ['mcsm_ppi2_affinity', 'interface_dist'] 
+        X_stabilityN = common_cols_stabiltyN + geneL_ppi2_st_cols
+        cols_to_mask = ['ligand_affinity_change', 'mcsm_ppi2_affinity']
+    
+    if gene.lower() in geneL_na:
+    #    X_stabilityN = common_cols_stabiltyN + ['mcsm_na_affinity'] 
+        geneL_na_st_cols =  ['mcsm_na_affinity'] 
+        X_stabilityN = common_cols_stabiltyN + geneL_na_st_cols
+        cols_to_mask = ['ligand_affinity_change', 'mcsm_na_affinity']
+    
+    if gene.lower() in geneL_na_ppi2:
+    #    X_stabilityN = common_cols_stabiltyN + ['mcsm_na_affinity'] + ['mcsm_ppi2_affinity', 'interface_dist'] 
+        geneL_na_ppi2_st_cols = ['mcsm_na_affinity'] + ['mcsm_ppi2_affinity', 'interface_dist'] 
+        X_stabilityN = common_cols_stabiltyN + geneL_na_ppi2_st_cols
+        cols_to_mask = ['ligand_affinity_change', 'mcsm_na_affinity', 'mcsm_ppi2_affinity']
+    
+
+    X_foldX_cols = [ 'electro_rr', 'electro_mm', 'electro_sm', 'electro_ss'
+    , 'disulfide_rr', 'disulfide_mm', 'disulfide_sm', 'disulfide_ss'
+    , 'hbonds_rr', 'hbonds_mm', 'hbonds_sm', 'hbonds_ss'
+    , 'partcov_rr', 'partcov_mm', 'partcov_sm', 'partcov_ss'
+    , 'vdwclashes_rr', 'vdwclashes_mm', 'vdwclashes_sm', 'vdwclashes_ss'
+    , 'volumetric_rr', 'volumetric_mm', 'volumetric_ss'
+    ]
+    
+    X_str =  ['rsa'
+               #, 'asa'
+               , 'kd_values'
+               , 'rd_values']    
+    
+    X_ssFN = X_stabilityN + X_str + X_foldX_cols
+    
+    X_evolFN =  ['consurf_score'
+               , 'snap2_score'
+               , 'provean_score']
+        
+    X_genomic_mafor =  ['maf'
+                    , 'logorI'
+                    # , 'or_rawI'
+                    # , 'or_mychisq'
+                    # , 'or_logistic'
+                    # , 'or_fisher'
+                    # , 'pval_fisher'
+                    ]
+    
+    X_genomic_linegae  = ['lineage_proportion'
+                          , 'dist_lineage_proportion'
+                          #, 'lineage' # could be included as a category but it has L2;L4  formatting
+                          , 'lineage_count_all'
+                          , 'lineage_count_unique'
+                          ]
+    
+    X_genomicFN = X_genomic_mafor + X_genomic_linegae
+    
+    #%% Construct numerical and categorical column names
+    # numerical feature names
+#    numerical_FN = common_cols_stabiltyN + foldX_cols + X_strFN + X_evolFN + X_genomicFN 
+    
+    numerical_FN = X_ssFN  + X_evolFN + X_genomicFN 
+
+    #categorical feature names
+    categorical_FN = ['ss_class'
+                # , 'wt_prop_water'
+                # , 'mut_prop_water'
+                # , 'wt_prop_polarity'
+                # , 'mut_prop_polarity'
+                # , 'wt_calcprop'
+                # , 'mut_calcprop'
+                , 'aa_prop_change'
+                , 'electrostatics_change'
+                , 'polarity_change'
+                , 'water_change'
+                , 'drtype_mode_labels' # beware then you can use it to predict 
+                #, 'active_aa_pos' # TODO?
+                 ]
+    ###########################################################################
+    #=======================
+    # Masking columns:
+    # (mCSM-lig, mCSM-NA, mCSM-ppi2) values for lig_dist >10
+    #=======================
+    #%% Masking columns 
+    # my_df_ml['mutationinformation'][my_df['ligand_distance']>10].value_counts()
+    # my_df_ml.groupby('mutationinformation')['ligand_distance'].apply(lambda x: (x>10)).value_counts()
+    
+    # my_df_ml.loc[(my_df_ml['ligand_distance'] > 10), 'ligand_affinity_change'] = 0
+    # (my_df_ml['ligand_affinity_change'] == 0).sum()
+    
+    my_df_ml['mutationinformation'][my_df_ml['ligand_distance']>10].value_counts()
+    my_df_ml.groupby('mutationinformation')['ligand_distance'].apply(lambda x: (x>10)).value_counts()
+    my_df_ml.loc[(my_df_ml['ligand_distance'] > 10), cols_to_mask].value_counts()
+    
+    # mask the column ligand distance > 10
+    my_df_ml.loc[(my_df_ml['ligand_distance'] > 10), cols_to_mask] = 0
+    (my_df_ml['ligand_affinity_change'] == 0).sum()
+    
+    mask_check = my_df_ml[['mutationinformation', 'ligand_distance'] + cols_to_mask]  
+    
+    # write file for check
+    mask_check.sort_values(by = ['ligand_distance'], ascending = True, inplace = True)
+    mask_check.to_csv(outdir + 'ml/' + gene.lower() + '_mask_check.csv')
+
+    #%% extracting dfs based on numerical, categorical column names
+    #----------------------------------
+    # WITHOUT the target var included
+    #----------------------------------
+    num_df = training_df[numerical_FN]
+    num_df.shape
+    
+    cat_df = training_df[categorical_FN]
+    cat_df.shape
+    
+    all_df = training_df[numerical_FN + categorical_FN]
+    all_df.shape
+    
+    #------------------------------
+    # WITH the target var included:
+        #'wtgt': with target
+    #------------------------------
+    # drug and dst_mode should be the same thing
+    num_df_wtgt = training_df[numerical_FN + ['dst_mode']]
+    num_df_wtgt.shape
+    
+    cat_df_wtgt = training_df[categorical_FN + ['dst_mode']]
+    cat_df_wtgt.shape
+    
+    all_df_wtgt = training_df[numerical_FN + categorical_FN + ['dst_mode']]
+    all_df_wtgt.shape
+    #%%########################################################################
+    #============
+    # ML data
+    #============
+    #------
+    # X: Training and Blind test (BTS)
+    #------
+    X     = all_df_wtgt[numerical_FN + categorical_FN] # training data ALL
+    X_bts = blind_test_df[numerical_FN + categorical_FN] # blind test data ALL
+    #X = all_df_wtgt[numerical_FN] # training numerical only
+    #X_bts = blind_test_df[numerical_FN] # blind test data numerical
+    
+    #------
+    # y
+    #------
+    y = all_df_wtgt['dst_mode'] # training data y
+    y_bts = blind_test_df['dst_mode'] # blind data test y
+
+    #X_bts_wt = blind_test_df[numerical_FN + ['dst_mode']] 
+    
+    # Quick check
+    #(X['ligand_affinity_change']==0).sum() == (X['ligand_distance']>10).sum()
+    for i in range(len(cols_to_mask)):
+        ind = i+1
+        print('\nindex:', i, '\nind:', ind)
+        print('\nMask count check:'
+              , (my_df_ml[cols_to_mask[i]]==0).sum() == (my_df_ml['ligand_distance']>10).sum()
+              )
+    
+    print('Original Data\n', Counter(y)
+          , 'Data dim:', X.shape)
+    
+    ###############################################################################
+    #%% 
+    ############################################################################
+    #                               RESAMPLING
+    ###############################################################################
+    #------------------------------
+    # Simple Random oversampling
+    # [Numerical + catgeorical]
+    #------------------------------
+    oversample = RandomOverSampler(sampling_strategy='minority')
+    X_ros, y_ros = oversample.fit_resample(X, y)
+    print('Simple Random OverSampling\n', Counter(y_ros))
+    print(X_ros.shape)
+    
+    #------------------------------
+    # Simple Random Undersampling
+    # [Numerical + catgeorical]
+    #------------------------------
+    undersample = RandomUnderSampler(sampling_strategy='majority')
+    X_rus, y_rus = undersample.fit_resample(X, y)
+    print('Simple Random UnderSampling\n', Counter(y_rus))
+    print(X_rus.shape)
+    
+    #------------------------------
+    # Simple combine ROS and RUS
+    # [Numerical + catgeorical]
+    #------------------------------
+    oversample = RandomOverSampler(sampling_strategy='minority')
+    X_ros, y_ros = oversample.fit_resample(X, y)
+    undersample = RandomUnderSampler(sampling_strategy='majority')
+    X_rouC, y_rouC = undersample.fit_resample(X_ros, y_ros)
+    print('Simple Combined Over and UnderSampling\n',  Counter(y_rouC))
+    print(X_rouC.shape)
+    
+    #------------------------------
+    # SMOTE_NC: oversampling 
+    # [numerical + categorical]
+    #https://stackoverflow.com/questions/47655813/oversampling-smote-for-binary-and-categorical-data-in-python
+    #------------------------------
+    # Determine categorical and numerical features
+    numerical_ix = X.select_dtypes(include=['int64', 'float64']).columns
+    numerical_ix
+    num_featuresL = list(numerical_ix)
+    numerical_colind = X.columns.get_indexer(list(numerical_ix) )
+    numerical_colind
+    
+    categorical_ix = X.select_dtypes(include=['object', 'bool']).columns
+    categorical_ix    
+    categorical_colind = X.columns.get_indexer(list(categorical_ix))
+    categorical_colind
+    
+    k_sm = 5 # 5 is deafult
+    sm_nc = SMOTENC(categorical_features=categorical_colind, k_neighbors = k_sm, **rs, **njobs)
+    X_smnc, y_smnc = sm_nc.fit_resample(X, y)
+    print('SMOTE_NC OverSampling\n', Counter(y_smnc))
+    print(X_smnc.shape)
+    globals().update(locals()) # TROLOLOLOLOLOLS
+    #print("i did a horrible hack :-)")
+    ###############################################################################
+    #%% SMOTE RESAMPLING for NUMERICAL ONLY*
+    # #------------------------------
+    # # SMOTE: Oversampling
+    # # [Numerical ONLY]
+    # #------------------------------
+    # k_sm = 1
+    # sm = SMOTE(sampling_strategy = 'auto', k_neighbors = k_sm, **rs)
+    # X_sm, y_sm = sm.fit_resample(X, y)
+    # print(X_sm.shape)
+    # print('SMOTE OverSampling\n', Counter(y_sm))
+    # y_sm_df = y_sm.to_frame()
+    # y_sm_df.value_counts().plot(kind = 'bar')
+    
+    # #------------------------------
+    # # SMOTE: Over + Undersampling COMBINED
+    # # [Numerical ONLY]
+    # #-----------------------------
+    # sm_enn = SMOTEENN(enn=EditedNearestNeighbours(sampling_strategy='all', **rs, **njobs ))
+    # X_enn, y_enn = sm_enn.fit_resample(X, y)
+    # print(X_enn.shape)
+    # print('SMOTE Over+Under Sampling combined\n', Counter(y_enn))
+    
+    ###############################################################################
+    # TODO: Find over and undersampling JUST for categorical data

scripts/ml/pnca_config.py

@@ -0,0 +1,188 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Sat May 28 05:25:30 2022
+
+@author: tanu
+"""
+
+import os
+
+gene  = 'pncA'
+drug  = 'pyrazinamide'
+#total_mtblineage_uc = 8
+
+homedir = os.path.expanduser("~")
+os.chdir( homedir + '/git/ML_AI_training/')
+
+#---------------------------
+# Version 1: no AAindex
+#from UQ_ML_data import *
+#setvars(gene,drug)
+#from UQ_ML_data import *
+#---------------------------
+
+from UQ_ML_data2 import *
+setvars(gene,drug)
+from UQ_ML_data2 import *
+
+# from YC run_all_ML: run locally
+#from UQ_yc_RunAllClfs import run_all_ML
+
+# TT run all ML clfs: baseline mode
+from UQ_MultModelsCl import MultModelsCl
+
+#%%###########################################################################
+
+print('\n#####################################################################\n')
+
+print('TESTING cmd:'
+      , '\nGene name:', gene
+      , '\nDrug name:', drug
+      , '\nTotal input features:', X.shape
+      , '\n', Counter(y))
+
+print('Strucutral features (n):'
+      , len(X_ssFN)
+      , '\nThese are:'
+      , '\nCommon stablity features:', X_stabilityN
+      , '\nFoldX columns:', X_foldX_cols
+      , '\nOther struc columns:', X_str
+      , '\n================================================================\n')
+
+print('AAindex features (n):'
+      , len(X_aaindexFN)
+      , '\nThese are:\n'
+      , X_aaindexFN
+      , '\n================================================================\n')      
+
+print('Evolutionary features (n):'
+      , len(X_evolFN)
+      , '\nThese are:\n'
+      , X_evolFN
+      , '\n================================================================\n')
+
+print('Genomic features (n):'
+      , len(X_genomicFN)
+      , '\nThese are:\n'
+      , X_genomic_mafor, '\n'
+      , X_genomic_linegae
+      , '\n================================================================\n')
+
+print('Categorical features (n):'
+      , len(categorical_FN)
+      , '\nThese are:\n'
+      , categorical_FN
+      , '\n================================================================\n')
+
+if ( len(X.columns) ==  len(X_ssFN) + len(X_aaindexFN) + len(X_evolFN) + len(X_genomicFN) + len(categorical_FN) ):
+    print('\nPass: No. of features match')
+else:
+    sys.exit('\nFail: Count of feature mismatch')
+
+print('\n#####################################################################\n')
+
+################################################################################
+#==================
+# Baseline models 
+#==================
+mm_skf_scoresD = MultModelsCl(input_df = X
+                                        , target = y
+                                        , var_type = 'mixed'
+                                        , skf_cv = skf_cv
+                                        , blind_test_input_df = X_bts
+                                        , blind_test_target = y_bts)
+
+baseline_all = pd.DataFrame(mm_skf_scoresD)
+baseline_all = baseline_all.T
+#baseline_train = baseline_all.filter(like='train_', axis=1)
+baseline_CT = baseline_all.filter(like='test_', axis=1)
+baseline_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True)
+
+baseline_BT = baseline_all.filter(like='bts_', axis=1)
+baseline_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True)
+
+# Write csv
+baseline_CT.to_csv(outdir + 'ml/' + gene.lower() + '_baseline_CT_allF.csv')
+baseline_BT.to_csv(outdir + 'ml/' + gene.lower() + '_baseline_BT_allF.csv')
+
+
+#%% SMOTE NC: Oversampling [Numerical + categorical]
+mm_skf_scoresD7 = MultModelsCl(input_df = X_smnc
+                                        , target = y_smnc
+                                        , var_type = 'mixed'
+                                        , skf_cv = skf_cv
+                                        , blind_test_input_df = X_bts
+                                        , blind_test_target = y_bts)
+smnc_all = pd.DataFrame(mm_skf_scoresD7)
+smnc_all = smnc_all.T
+
+smnc_CT = smnc_all.filter(like='test_', axis=1)
+smnc_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True)
+
+smnc_BT = smnc_all.filter(like='bts_', axis=1)
+smnc_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True)
+
+# Write csv
+smnc_CT.to_csv(outdir + 'ml/' + gene.lower() + '_smnc_CT_allF.csv')
+smnc_BT.to_csv(outdir + 'ml/' + gene.lower() + '_smnc_BT_allF.csv')
+
+#%% ROS: Numerical + categorical
+mm_skf_scoresD3 = MultModelsCl(input_df = X_ros
+                                        , target = y_ros
+                                        , var_type = 'mixed'
+                                        , skf_cv = skf_cv
+                                        , blind_test_input_df = X_bts
+                                        , blind_test_target = y_bts)
+ros_all = pd.DataFrame(mm_skf_scoresD3)
+ros_all = ros_all.T
+
+ros_CT  = ros_all.filter(like='test_', axis=1)
+ros_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True)
+
+ros_BT  = ros_all.filter(like='bts_', axis=1)
+ros_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True)
+
+# Write csv
+ros_CT.to_csv(outdir + 'ml/' + gene.lower() + '_ros_CT_allF.csv')
+ros_BT.to_csv(outdir + 'ml/' + gene.lower() + '_ros_BT_allF.csv')
+
+#%% RUS: Numerical + categorical
+mm_skf_scoresD4 = MultModelsCl(input_df = X_rus
+                                        , target = y_rus
+                                        , var_type = 'mixed'
+                                        , skf_cv = skf_cv
+                                        , blind_test_input_df = X_bts
+                                        , blind_test_target = y_bts)
+rus_all = pd.DataFrame(mm_skf_scoresD4)
+rus_all = rus_all.T
+
+rus_CT = rus_all.filter(like='test_', axis=1)
+rus_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True)
+
+rus_BT = rus_all.filter(like='bts_' , axis=1)
+rus_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True)
+
+# Write csv
+rus_CT.to_csv(outdir + 'ml/' + gene.lower() + '_rus_CT_allF.csv')
+rus_BT.to_csv(outdir + 'ml/' + gene.lower() + '_rus_BT_allF.csv')
+
+#%% ROS + RUS Combined: Numerical + categorical
+mm_skf_scoresD8 = MultModelsCl(input_df = X_rouC
+                                        , target = y_rouC
+                                        , var_type = 'mixed'
+                                        , skf_cv = skf_cv
+                                        , blind_test_input_df = X_bts
+                                        , blind_test_target = y_bts)
+rouC_all = pd.DataFrame(mm_skf_scoresD8)
+rouC_all = rouC_all.T
+
+rouC_CT = rouC_all.filter(like='test_', axis=1)
+rouC_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True)
+
+rouC_BT = rouC_all.filter(like='bts_', axis=1)
+rouC_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True)
+
+# Write csv
+rouC_CT.to_csv(outdir + 'ml/' + gene.lower() + '_rouC_CT_allF.csv')
+rouC_BT.to_csv(outdir + 'ml/' + gene.lower() + '_rouC_BT_allF.csv')