diff --git a/scripts/ml/MultModelsCl_dissected.py b/scripts/ml/MultModelsCl_dissected.py deleted file mode 100644 index 5a88c05..0000000 --- a/scripts/ml/MultModelsCl_dissected.py +++ /dev/null @@ -1,318 +0,0 @@ -#!/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 - -# added -from sklearn.model_selection import train_test_split, cross_validate, cross_val_score, LeaveOneOut, KFold, RepeatedKFold, cross_val_predict - -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 -import argparse -import re - -#%% GLOBALS -rs = {'random_state': 42} -njobs = {'n_jobs': 10} - -scoring_fn = ({ 'mcc' : make_scorer(matthews_corrcoef) - , 'fscore' : make_scorer(f1_score) - , 'precision' : make_scorer(precision_score) - , 'recall' : make_scorer(recall_score) - , 'accuracy' : make_scorer(accuracy_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 - , add_cm = True # adds confusion matrix based on cross_val_predict - , add_yn = True # adds target var class numbers - , 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 - #====================================================== - models = [('Logistic Regression' , LogisticRegression(**rs) ) - , ('Logistic RegressionCV' , LogisticRegressionCV(**rs) ) - , ('Gaussian NB' , GaussianNB() ) - , ('Naive Bayes' , BernoulliNB() ) - # , ('K-Nearest Neighbors' , KNeighborsClassifier() ) - # , ('SVC' , SVC(**rs) ) - # , ('MLP' , MLPClassifier(max_iter = 500, **rs) ) - # , ('Decision Tree' , DecisionTreeClassifier(**rs) ) - # , ('Extra Trees' , ExtraTreesClassifier(**rs) ) - # , ('Extra Tree' , ExtraTreeClassifier(**rs) ) - # , ('Random Forest' , RandomForestClassifier(**rs, n_estimators = 1000 ) ) - # , ('Random Forest2' , RandomForestClassifier(min_samples_leaf = 5 - # , n_estimators = 1000 - # , bootstrap = True - # , oob_score = True - # , **njobs - # , **rs - # , max_features = 'auto') ) - # , ('XGBoost' , XGBClassifier(**rs, verbosity = 0, use_label_encoder =False) ) - # , ('LDA' , LinearDiscriminantAnalysis() ) - # , ('Multinomial' , MultinomialNB() ) - # , ('Passive Aggresive' , PassiveAggressiveClassifier(**rs, **njobs) ) - # , ('Stochastic GDescent' , SGDClassifier(**rs, **njobs) ) - # , ('AdaBoost Classifier' , AdaBoostClassifier(**rs) ) - # , ('Bagging Classifier' , BaggingClassifier(**rs, **njobs, bootstrap = True, oob_score = True) ) - # , ('Gaussian Process' , GaussianProcessClassifier(**rs) ) - # , ('Gradient Boosting' , GradientBoostingClassifier(**rs) ) - # , ('QDA' , QuadraticDiscriminantAnalysis() ) - # , ('Ridge Classifier' , RidgeClassifier(**rs) ) - # , ('Ridge ClassifierCV' , RidgeClassifierCV(cv = 10) ) - ] - - mm_skf_scoresD = {} - - print('\n==============================================================\n' - , '\nRunning several classification models (n):', len(models) - ,'\nList of models:') - for m in models: - print(m) - print('\n================================================================\n') - - index = 1 - for model_name, model_fn in models: - print('\nRunning classifier:', index - , '\nModel_name:' , model_name - , '\nModel func:' , model_fn) - index = index+1 - - model_pipeline = Pipeline([ - ('prep' , col_transform) - , ('model' , model_fn)]) - - print('\nRunning model pipeline:', model_pipeline) - skf_cv_modD = cross_validate(model_pipeline - , input_df - , target - , cv = skf_cv - , scoring = scoring_fn - , return_train_score = True) - - ####################################################################### - #====================================================== - # Option: Add confusion matrix from cross_val_predict - # Understand and USE with caution - # cross_val_score, cross_val_predict, "Passing these predictions into an evaluation metric may not be a valid way to measure generalization performance. Results can differ from cross_validate and cross_val_score unless all tests sets have equal size and the metric decomposes over samples." - # https://stackoverflow.com/questions/65645125/producing-a-confusion-matrix-with-cross-validate - #====================================================== - if add_cm: - - #----------------------------------------------------------- - # Initialise dict of Confusion Matrix (cm) - #----------------------------------------------------------- - cmD = {} - - # Calculate cm - y_pred = cross_val_predict(model_pipeline, input_df, target, cv = skf_cv, **njobs) - #_tn, _fp, _fn, _tp = confusion_matrix(y_pred, y).ravel() # internally - tn, fp, fn, tp = confusion_matrix(y_pred, target).ravel() - - # Build dict - - cmD = {'TN' : tn - , 'FP': fp - , 'FN': fn - , 'TP': tp} - #--------------------------------- - # Update cv dict with cmD and tbtD - #---------------------------------- - skf_cv_modD.update(cmD) - else: - skf_cv_modD = skf_cv_modD - ####################################################################### - #============================================= - # Option: Add targety numbers for data - #============================================= - if add_yn: - - #----------------------------------------------------------- - # Initialise dict of target numbers: training and blind (tbt) - #----------------------------------------------------------- - tbtD = {} - - # training y - tyn = Counter(target) - tyn_neg = tyn[0] - tyn_pos = tyn[1] - - # blind test y - btyn = Counter(blind_test_target) - btyn_neg = btyn[0] - btyn_pos = btyn[1] - - # Build dict - tbtD = {'trainingY_neg' : tyn_neg - , 'trainingY_pos' : tyn_pos - , 'blindY_neg' : btyn_neg - , 'blindY_pos' : btyn_pos} - - #--------------------------------- - # Update cv dict with cmD and tbtD - #---------------------------------- - skf_cv_modD.update(tbtD) - else: - skf_cv_modD = skf_cv_modD - - ####################################################################### - #============================== - # Extract mean values for CV - #============================== - mm_skf_scoresD[model_name] = {} - - for key, value in skf_cv_modD.items(): - print('\nkey:', key, '\nvalue:', value) - print('\nmean value:', np.mean(value)) - mm_skf_scoresD[model_name][key] = round(np.mean(value),2) - - #return(mm_skf_scoresD) -#%% - #========================= - # Blind test: BTS results - #========================= - # Build the final results with all scores for the 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) - - 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_jcc'] = round(jaccard_score(blind_test_target, bts_predict),2) - #mm_skf_scoresD[model_name]['diff_mcc'] = train_test_diff_MCC - - return(mm_skf_scoresD) - diff --git a/scripts/ml/ml_data_dissected.py b/scripts/ml/ml_data_dissected.py deleted file mode 100644 index d1daa2c..0000000 --- a/scripts/ml/ml_data_dissected.py +++ /dev/null @@ -1,791 +0,0 @@ -#!/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 - import argparse - import re - #%% 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 - #%% DONE: 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("~") - - geneL_basic = ['pnca'] - geneL_na = ['gid'] - geneL_na_ppi2 = ['rpob'] - geneL_ppi2 = ['alr', 'embb', 'katg'] - - #num_type = ['int64', 'float64'] - num_type = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64'] - cat_type = ['object', 'bool'] - - #============== - # directories - #============== - datadir = homedir + '/git/Data/' - indir = datadir + drug + '/input/' - outdir = datadir + drug + '/output/' - - #======= - # input - #======= - - #--------- - # File 1 - #--------- - infile_ml1 = outdir + gene.lower() + '_merged_df3.csv' - #infile_ml2 = outdir + gene.lower() + '_merged_df2.csv' - - my_features_df = pd.read_csv(infile_ml1, index_col = 0) - my_features_df = my_features_df .reset_index(drop = True) - my_features_df.index - - my_features_df.dtypes - mycols = my_features_df.columns - - #--------- - # File 2 - #--------- - infile_aaindex = outdir + 'aa_index/' + gene.lower() + '_aa.csv' - aaindex_df = pd.read_csv(infile_aaindex, index_col = 0) - aaindex_df.dtypes - - #----------- - # check for non-numerical columns - #----------- - if any(aaindex_df.dtypes==object): - print('\naaindex_df contains non-numerical data') - - aaindex_df_object = aaindex_df.select_dtypes(include = cat_type) - print('\nTotal no. of non-numerial columns:', len(aaindex_df_object.columns)) - - expected_aa_ncols = len(aaindex_df.columns) - len(aaindex_df_object.columns) - - #----------- - # Extract numerical data only - #----------- - print('\nSelecting numerical data only') - aaindex_df = aaindex_df.select_dtypes(include = num_type) - - #--------------------------- - # aaindex: sanity check 1 - #--------------------------- - if len(aaindex_df.columns) == expected_aa_ncols: - print('\nPASS: successfully selected numerical columns only for aaindex_df') - else: - print('\nFAIL: Numbers mismatch' - , '\nExpected ncols:', expected_aa_ncols - , '\nGot:', len(aaindex_df.columns)) - - #--------------- - # check for NA - #--------------- - print('\nNow checking for NA in the remaining aaindex_cols') - c1 = aaindex_df.isna().sum() - c2 = c1.sort_values(ascending=False) - print('\nCounting aaindex_df cols with NA' - , '\nncols with NA:', sum(c2>0), 'columns' - , '\nDropping these...' - , '\nOriginal ncols:', len(aaindex_df.columns) - ) - aa_df = aaindex_df.dropna(axis=1) - - print('\nRevised df ncols:', len(aa_df.columns)) - - c3 = aa_df.isna().sum() - c4 = c3.sort_values(ascending=False) - - print('\nChecking NA in revised df...') - - if sum(c4>0): - sys.exit('\nFAIL: aaindex_df still contains cols with NA, please check and drop these before proceeding...') - else: - print('\nPASS: cols with NA successfully dropped from aaindex_df' - , '\nProceeding with combining aa_df with other features_df') - - #--------------------------- - # aaindex: sanity check 2 - #--------------------------- - expected_aa_ncols2 = len(aaindex_df.columns) - sum(c2>0) - if len(aa_df.columns) == expected_aa_ncols2: - print('\nPASS: ncols match' - , '\nExpected ncols:', expected_aa_ncols2 - , '\nGot:', len(aa_df.columns)) - else: - print('\nFAIL: Numbers mismatch' - , '\nExpected ncols:', expected_aa_ncols2 - , '\nGot:', len(aa_df.columns)) - - # Important: need this to identify aaindex cols - aa_df_cols = aa_df.columns - print('\nTotal no. of columns in clean aa_df:', len(aa_df_cols)) - - ############################################################################### - #%% Combining my_features_df and aaindex_df - #=========================== - # Merge my_df + aaindex_df - #=========================== - - if aa_df.columns[aa_df.columns.isin(my_features_df.columns)] == my_features_df.columns[my_features_df.columns.isin(aa_df.columns)]: - print('\nMerging on column: mutationinformation') - - if len(my_features_df) == len(aa_df): - expected_nrows = len(my_features_df) - print('\nProceeding to merge, expected nrows in merged_df:', expected_nrows) - else: - sys.exit('\nNrows mismatch, cannot merge. Please check' - , '\nnrows my_df:', len(my_features_df) - , '\nnrows aa_df:', len(aa_df)) - - #----------------- - # Reset index: mutationinformation - # Very important for merging - #----------------- - aa_df = aa_df.reset_index() - - expected_ncols = len(my_features_df.columns) + len(aa_df.columns) - 1 # for the no. of merging col - - #----------------- - # Merge: my_features_df + aa_df - #----------------- - merged_df = pd.merge(my_features_df - , aa_df - , on = 'mutationinformation') - - #--------------------------- - # aaindex: sanity check 3 - #--------------------------- - if len(merged_df.columns) == expected_ncols: - print('\nPASS: my_features_df and aa_df successfully combined' - , '\nnrows:', len(merged_df) - , '\nncols:', len(merged_df.columns)) - else: - sys.exit('\nFAIL: could not combine my_features_df and aa_df' - , '\nCheck dims and merging cols!') - - #-------- - # Reassign so downstream code doesn't need to change - #-------- - my_df = merged_df.copy() - - #%% Data: my_df - # Check if non structural pos have crept in - # IDEALLY remove from source! But for rpoB do it here - # 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) - - # FIXED: complete data for all muts inc L114M, F115L, V123L, V125I, V131M - # 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) - - # # Sanity check for non-structural positions - # print('\nChecking for non-structural postions') - # na_index = my_df['mutationinformation'].index[my_df['ligand_distance'].apply(np.isnan)] - # if len(na_index) > 0: - # print('\nNon-structural positions detected for gene:', gene.lower() - # , '\nTotal number of these detected:', len(na_index) - # , '\These are at index:', na_index - # , '\nOriginal nrows:', len(my_df) - # , '\nDropping these...') - # my_df = my_df.drop(index=na_index) - # print('\nRevised nrows:', len(my_df)) - # else: - # print('\nNo non-structural positions detected for gene:', gene.lower() - # , '\nnrows:', len(my_df)) - - - ########################################################################### - #%% 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 - ########################################################################### - #%% Active site annotation column - # change from numberic to categorical - - if my_df['active_site'].dtype in num_type: - my_df['active_site'] = my_df['active_site'].astype(object) - my_df['active_site'].dtype - #%% 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] - #-------------------- - # Impute OR values - #-------------------- - #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=3, 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() # should be 0 - - len(my_df.columns) - len(mydf_imputed.columns) - - #----------------------------------------- - # REASSIGN my_df after imputing OR values - #----------------------------------------- - my_df = mydf_imputed.copy() - - if my_df['logorI'].isna().sum() == 0: - print('\nPASS: OR values imputed, data ready for ML') - else: - sys.exit('\nFAIL: something went wrong, Data not ready for ML. Please check upstream!') - - #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - #--------------------------------------- - # TODO: try other imputation like MICE - #--------------------------------------- - #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - #%% Data for ML - #========================== - # Data for ML - #========================== - my_df_ml = my_df.copy() - - # Build column names to mask for affinity chanhes - if gene.lower() in geneL_basic: - #X_stabilityN = common_cols_stabiltyN - gene_affinity_colnames = []# not needed as its the common ones - cols_to_mask = ['ligand_affinity_change'] - - if gene.lower() in geneL_ppi2: - gene_affinity_colnames = ['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: - gene_affinity_colnames = ['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: - gene_affinity_colnames = ['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'] - - #======================= - # Masking columns: - # (mCSM-lig, mCSM-NA, mCSM-ppi2) values for lig_dist >10 - #======================= - 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 mcsm affinity related columns where 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') - #=================================================== - ############################################################################### - #%% Feature groups (FG): Build X for Input ML - ############################################################################ - #=========================== - # FG1: Evolutionary features - #=========================== - X_evolFN = ['consurf_score' - , 'snap2_score' - , 'provean_score'] - - ############################################################################### - #======================== - # FG2: Stability features - #======================== - #-------- - # common - #-------- - X_common_stability_Fnum = [ - 'duet_stability_change' - , 'ddg_foldx' - , 'deepddg' - , 'ddg_dynamut2' - , 'contacts'] - #-------- - # FoldX - #-------- - X_foldX_Fnum = [ '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_stability_FN = X_common_stability_Fnum + X_foldX_Fnum - - ############################################################################### - #=================== - # FG3: Affinity features - #=================== - common_affinity_Fnum = ['ligand_distance' - , 'ligand_affinity_change' - , 'mmcsm_lig'] - - # if gene.lower() in geneL_basic: - # X_affinityFN = common_affinity_Fnum - # else: - # X_affinityFN = common_affinity_Fnum + gene_affinity_colnames - - X_affinityFN = common_affinity_Fnum + gene_affinity_colnames - - ############################################################################### - #============================ - # FG4: Residue level features - #============================ - #----------- - # AA index - #----------- - X_aaindex_Fnum = list(aa_df_cols) - print('\nTotal no. of features for aaindex:', len(X_aaindex_Fnum)) - - #----------------- - # surface area - # depth - # hydrophobicity - #----------------- - X_str_Fnum = ['rsa' - #, 'asa' - , 'kd_values' - , 'rd_values'] - - #--------------------------- - # Other aa properties - # active site indication - #--------------------------- - X_aap_Fcat = ['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' - , 'active_site'] - - - X_resprop_FN = X_aaindex_Fnum + X_str_Fnum + X_aap_Fcat - ############################################################################### - #======================== - # FG5: Genomic features - #======================== - X_gn_mafor_Fnum = ['maf' - , 'logorI' - # , 'or_rawI' - # , 'or_mychisq' - # , 'or_logistic' - # , 'or_fisher' - # , 'pval_fisher' - ] - - X_gn_linegae_Fnum = ['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_gn_Fcat = ['drtype_mode_labels' # beware then you can't use it to predict [USED it for uq_v1, not v2] - #, 'gene_name' # will be required for the combined stuff - ] - - X_genomicFN = X_gn_mafor_Fnum + X_gn_linegae_Fnum + X_gn_Fcat - ############################################################################### - #======================== - # FG6 collapsed: Structural : Atability + Affinity + ResidueProp - #======================== - X_structural_FN = X_stability_FN + X_affinityFN + X_resprop_FN - - ############################################################################### - #======================== - # BUILDING all features - #======================== - all_featuresN = X_evolFN + X_structural_FN + X_genomicFN - - ############################################################################### - #%% Define training and test data - #====================================================== - # Training and BLIND test set [UQ]: actual vs imputed - # No aa index but active_site included - # dst with actual values : training set - # dst with imputed values : blind test - #====================================================== - my_df_ml[drug].isna().sum() #'na' ones are the blind_test set - - 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 - - # Target 1: dst_mode - training_df[drug].value_counts() - training_df['dst_mode'].value_counts() - - #################################################################### - #============ - # ML data - #============ - #------ - # X: Training and Blind test (BTS) - #------ - X = training_df[all_featuresN] - X_bts = blind_test_df[all_featuresN] - - #------ - # y - #------ - y = training_df['dst_mode'] - y_bts = blind_test_df['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) - - yc1 = Counter(y) - yc1_ratio = yc1[0]/yc1[1] - - yc2 = Counter(y_bts) - yc2_ratio = yc2[0]/yc2[1] - - ############################################################################### - #====================================================== - # Determine categorical and numerical features - #====================================================== - numerical_cols = X.select_dtypes(include=['int64', 'float64']).columns - numerical_cols - categorical_cols = X.select_dtypes(include=['object', 'bool']).columns - categorical_cols - - ################################################################################ - # IMPORTANT sanity checks - if len(X.columns) == len(X_evolFN) + len(X_stability_FN) + len(X_affinityFN) + len(X_resprop_FN) + len(X_genomicFN): - print('\nPASS: ML data with input features, training and test generated...' - , '\n\nTotal no. of input features:' , len(X.columns) - , '\n--------No. of numerical features:' , len(numerical_cols) - , '\n--------No. of categorical features:' , len(categorical_cols) - - , '\n\nTotal no. of evolutionary features:' , len(X_evolFN) - - , '\n\nTotal no. of stability features:' , len(X_stability_FN) - , '\n--------Common stabilty cols:' , len(X_common_stability_Fnum) - , '\n--------Foldx cols:' , len(X_foldX_Fnum) - - , '\n\nTotal no. of affinity features:' , len(X_affinityFN) - , '\n--------Common affinity cols:' , len(common_affinity_Fnum) - , '\n--------Gene specific affinity cols:' , len(gene_affinity_colnames) - - , '\n\nTotal no. of residue level features:', len(X_resprop_FN) - , '\n--------AA index cols:' , len(X_aaindex_Fnum) - , '\n--------Residue Prop cols:' , len(X_str_Fnum) - , '\n--------AA change Prop cols:' , len(X_aap_Fcat) - - , '\n\nTotal no. of genomic features:' , len(X_genomicFN) - , '\n--------MAF+OR cols:' , len(X_gn_mafor_Fnum) - , '\n--------Lineage cols:' , len(X_gn_linegae_Fnum) - , '\n--------Other cols:' , len(X_gn_Fcat) - ) - else: - print('\nFAIL: numbers mismatch' - , '\nExpected:',len(X_evolFN) + len(X_stability_FN) + len(X_affinityFN) + len(X_resprop_FN) + len(X_genomicFN) - , '\nGot:', len(X.columns)) - sys.exit() - ############################################################################### - print('\n-------------------------------------------------------------' - , '\nSuccessfully split data: ALL features' - , '\nactual values: training set' - , '\nimputed values: blind test set' - - , '\n\nTotal data size:', len(X) + len(X_bts) - - , '\n\nTrain data size:', X.shape - , '\ny_train numbers:', yc1 - - , '\n\nTest data size:', X_bts.shape - , '\ny_test_numbers:', yc2 - - , '\n\ny_train ratio:',yc1_ratio - , '\ny_test ratio:', yc2_ratio - , '\n-------------------------------------------------------------' - ) - - ########################################################################### - #%% - ########################################################################### - # 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