diff --git a/scripts/ml/MultModelsCl_dissected.py b/scripts/ml/MultModelsCl_dissected.py new file mode 100644 index 0000000..d8804af --- /dev/null +++ b/scripts/ml/MultModelsCl_dissected.py @@ -0,0 +1,284 @@ +#!/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 + +#%% 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_dissected(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 + 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) + + #---------- + # check 1 + #---------- + foo_df = pd.DataFrame.from_dict(skf_cv_modD, orient ='index') + #foo_df = pd.DataFrame.from_dict(skf_cv_modD) + + #=================== + # Confusion matrix: Not an easy problem to solve! STILL DOING it, 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 + #=================== + y_pred = cross_val_predict(model_pipeline, input_df, target, cv = 10, **njobs) + #_tn, _fp, _fn, _tp = confusion_matrix(y_pred, y).ravel() # internally + tn, fp, fn, tp = confusion_matrix(y_pred, target).ravel() + # create a dict of confusion matrix that can be appended to the one above + # cmD = {'TN' : np.array(tn) + # , 'FP': np.array(fp) + # , 'FN': np.array(fn) + # , 'TP': np.array(tp)} + + cmD = {'TN' : tn + , 'FP': fp + , 'FN': fn + , 'TP': tp} + skf_cv_modD.update(cmD) + + #---------- + # check 2 + #---------- + #foo2_df = pd.DataFrame.from_dict(skf_cv_modD, orient ='index') + #foo_df = pd.DataFrame.from_dict(skf_cv_modD) + + 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 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) + diff --git a/scripts/ml/ml_data_dissected.py b/scripts/ml/ml_data_dissected.py new file mode 100644 index 0000000..4bd588c --- /dev/null +++ b/scripts/ml/ml_data_dissected.py @@ -0,0 +1,699 @@ +#!/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 +#%% 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 +#========================== +my_df_ml = my_df.copy() + +#%% 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 + +#X_aaindexFN = list(aa_df_cols) + +#print('\nTotal no. of features for aaindex:', len(X_aaindexFN)) + +# numerical feature names [NO aa_index] +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't use it to predict [USED it for uq_v1, not v2] + , 'active_site' #[didn't use it for uq_v1] + #, 'gene_name' # will be required for the combined stuff + ] +#---------------------------------------------- +# count numerical and categorical features +#---------------------------------------------- + +print('\nNo. of numerical features:', len(numerical_FN) + , '\nNo. of categorical features:', len(categorical_FN)) + +########################################################################### +#======================= +# Masking columns: +# (mCSM-lig, mCSM-NA, mCSM-ppi2) values for lig_dist >10 +#======================= +# 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 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') + +#=================================================== +# 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[numerical_FN + categorical_FN] +X_bts = blind_test_df[numerical_FN + categorical_FN] + +#------ +# 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] + +print('\n-------------------------------------------------------------' + , '\nSuccessfully split data: UQ [no aa_index but active site included] training' + , '\nactual values: training set' + , '\nimputed values: blind test set' + , '\nTrain data size:', X.shape + , '\nTest data size:', X_bts.shape + , '\ny_train numbers:', yc1 + , '\ny_train ratio:',yc1_ratio + , '\n' + , '\ny_test_numbers:', yc2 + , '\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 diff --git a/scripts/ml/pnca_config_dissected.py b/scripts/ml/pnca_config_dissected.py new file mode 100644 index 0000000..3c3868d --- /dev/null +++ b/scripts/ml/pnca_config_dissected.py @@ -0,0 +1,210 @@ +#!/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/LSHTM_analysis/scripts/ml/') + +#--------------------------- +# Version 1: no AAindex +#from UQ_ML_data import * +#setvars(gene,drug) +#from UQ_ML_data import * +#--------------------------- + +from ml_data_dissected import * +setvars(gene,drug) +from ml_data_dissected import * + +# from YC run_all_ML: run locally +#from UQ_yc_RunAllClfs import run_all_ML + +# TT run all ML clfs: baseline mode +from MultModelsCl_dissected import MultModelsCl_dissected + +############################################################################ +print('\n#####################################################################\n' + , '\nRunning ML analysis: UQ [without AA index but with active site annotations]' + , '\nGene name:', gene + , '\nDrug name:', drug) + +#================== +# Specify outdir +#================== + +outdir_ml = outdir + 'ml/uq_v1/dissected' + +print('\nOutput directory:', outdir_ml) + +#%%########################################################################### +print('\nSanity checks:' + , '\nTotal input features:', len(X.columns) + , '\n' + , '\nTraining data size:', X.shape + , '\nTest data size:', X_bts.shape + , '\n' + , '\nTarget feature numbers (training data):', Counter(y) + , '\nTarget features ratio (training data:', yc1_ratio + , '\n' + , '\nTarget feature numbers (test data):', Counter(y_bts) + , '\nTarget features ratio (test data):', yc2_ratio + + , '\n\n#####################################################################\n') + +print('\n================================================================\n') + +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) ): +if ( len(X.columns) == len(X_ssFN) + 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')