diff --git a/scripts/ml/MultModelsCl.py b/scripts/ml/MultModelsCl.py index 74e2482..5a88c05 100755 --- a/scripts/ml/MultModelsCl.py +++ b/scripts/ml/MultModelsCl.py @@ -41,6 +41,9 @@ 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 @@ -69,18 +72,20 @@ 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) +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 @@ -98,6 +103,8 @@ jacc_score_fn = {'jcc': make_scorer(jaccard_score)} 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']): ''' @@ -116,14 +123,18 @@ def MultModelsCl(input_df, target, skf_cv 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)] @@ -138,42 +149,42 @@ def MultModelsCl(input_df, target, skf_cv , remainder='passthrough') #====================================================== - # Specify multiple Classification models + # 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) ) + # , ('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:') @@ -198,8 +209,74 @@ def MultModelsCl(input_df, target, skf_cv , target , cv = skf_cv , scoring = scoring_fn - , return_train_score = True) + , 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 #============================== @@ -207,15 +284,15 @@ def MultModelsCl(input_df, target, skf_cv for key, value in skf_cv_modD.items(): print('\nkey:', key, '\nvalue:', value) - print('\nmean value:', mean(value)) - mm_skf_scoresD[model_name][key] = round(mean(value),2) - + 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 + # 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) @@ -225,28 +302,16 @@ def MultModelsCl(input_df, target, skf_cv 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 + # 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]['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/MultModelsCl_dissected.py b/scripts/ml/MultModelsCl_dissected.py index b93ada6..5a88c05 100644 --- a/scripts/ml/MultModelsCl_dissected.py +++ b/scripts/ml/MultModelsCl_dissected.py @@ -72,6 +72,8 @@ 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} @@ -98,7 +100,7 @@ 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 +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 @@ -299,6 +301,10 @@ def MultModelsCl_dissected(input_df, target, skf_cv 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) diff --git a/scripts/ml/Mult_dissected_CALL.py b/scripts/ml/Mult_dissected_CALL.py index 000e302..bc4bb07 100644 --- a/scripts/ml/Mult_dissected_CALL.py +++ b/scripts/ml/Mult_dissected_CALL.py @@ -29,46 +29,74 @@ score_type_ordermapD = { 'mcc' : 1 , 'fit_time' : 16 , 'score_time' : 17 } - - +############################################################################### #================== -# Baseline models +# Specify outdir #================== -# cm_di2 = MultModelsCl_dissected(input_df = X -# , target = y -# , var_type = 'mixed' -# , skf_cv = skf_cv -# , blind_test_input_df = X_bts -# , blind_test_target = y_bts -# , add_cm = True -# , add_yn = True) -# baseline_all2 = pd.DataFrame(cm_di2) -# baseline_all2T = baseline_all2.T -# baseline_CTBT2 = baseline_all2T.filter(regex = 'test_.*|bts_.*|TN|FP|FN|TP|.*_neg|.*_pos' , axis = 1) +outdir_ml = outdir + 'ml/uq_v1/fgs/' +print('\nOutput directory:', outdir_ml) +outFile = outdir_ml + gene.lower() + '_baseline_FG.csv' +#================== +# other vars +#================== +tts_split_name = 'original' +sampling_type_name = 'none' + +############################################################################### #================ -# Stability cols -#================ +# Evolutionary +# X_evolFN +#================ +feature_gp_nameEV = 'evolutionary' +n_featuresEV = len(X_evolFN) +scores_mmEV = MultModelsCl_dissected(input_df = X[X_evolFN] + , target = y + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts[X_evolFN] + , blind_test_target = y_bts + , add_cm = True + , add_yn = True) -#================ -# Affinity cols -#================ +baseline_allEV = pd.DataFrame(scores_mmEV) +baseline_EV = baseline_allEV.filter(regex = 'bts_.*|test_.*|.*_time|TN|FP|FN|TP|.*_neg|.*_pos', axis = 0) +baseline_EV = baseline_EV.reset_index() +baseline_EV.rename(columns = {'index': 'original_names'}, inplace = True) -#================ -# Residue level -#================ +# Indicate whether BT or CT +bt_pattern = re.compile(r'bts_.*') +baseline_EV['data_source'] = baseline_EV.apply(lambda row: 'BT' if bt_pattern.search(row.original_names) else 'CV' , axis = 1) +baseline_EV['score_type'] = baseline_EV['original_names'].str.replace('bts_|test_', '', regex = True) +score_type_uniqueN = set(baseline_EV['score_type']) +cL1 = list(score_type_ordermapD.keys()) +cL2 = list(score_type_uniqueN) + +if set(cL1).issubset(cL2): + print('\nPASS: sorting df by score that is mapped onto the order I want') + baseline_EV['score_order'] = baseline_EV['score_type'].map(score_type_ordermapD) + baseline_EV.sort_values(by = ['data_source', 'score_order'], ascending = [True, True], inplace = True) +else: + sys.exit('\nFAIL: could not sort df as score mapping for ordering failed') + +baseline_EV['feature_group'] = feature_gp_nameEV +baseline_EV['sampling_type'] = sampling_type_name +baseline_EV['tts_split'] = tts_split_name +baseline_EV['n_features'] = n_featuresEV +############################################################################### #================ # Genomics # X_genomicFN #================ -feature_gp_name = 'genomics' +feature_gp_nameGN = 'genomics' +n_featuresGN = len(X_genomicFN) -scores_mm_gn = MultModelsCl_dissected(input_df = X[X_genomicFN] +scores_mmGN = MultModelsCl_dissected(input_df = X[X_genomicFN] , target = y , var_type = 'mixed' , skf_cv = skf_cv @@ -77,9 +105,9 @@ scores_mm_gn = MultModelsCl_dissected(input_df = X[X_genomicFN] , add_cm = True , add_yn = True) -baseline_all_gn = pd.DataFrame(scores_mm_gn) +baseline_allGN = pd.DataFrame(scores_mmGN) -baseline_GN = baseline_all_gn.filter(regex = 'bts_.*|test_.*|.*_time|TN|FP|FN|TP|.*_neg|.*_pos', axis = 0) +baseline_GN = baseline_allGN.filter(regex = 'bts_.*|test_.*|.*_time|TN|FP|FN|TP|.*_neg|.*_pos', axis = 0) baseline_GN = baseline_GN.reset_index() baseline_GN.rename(columns = {'index': 'original_names'}, inplace = True) @@ -100,47 +128,340 @@ if set(cL1).issubset(cL2): else: sys.exit('\nFAIL: could not sort df as score mapping for ordering failed') -baseline_GN['feature_group'] = feature_gp_name - -#------------- -# Blind test -#------------- -baseline_BT = baseline_all_gn.filter(regex = 'bts_', axis = 0) -baseline_BT = baseline_BT.reset_index() -baseline_BT.rename(columns = {'index': 'original_names'}, inplace = True) -baseline_BT['score_type'] = baseline_BT['original_names'] -baseline_BT['score_type'] = baseline_BT['score_type'].str.replace('bts_*', '', regex = True) -baseline_BT['data_source'] = 'BT_score' +baseline_GN['feature_group'] = feature_gp_nameGN +baseline_GN['sampling_type'] = sampling_type_name +baseline_GN['tts_split'] = tts_split_name +baseline_GN['n_features'] = n_featuresGN +############################################################################### +#all_featuresN = X_evolFN + X_structural_FN + X_genomicFN +# X_structural_FN = X_stability_FN + X_affinityFN + X_resprop_FN +# X_resprop_FN = X_aaindex_Fnum + X_str_Fnum + X_aap_Fcat +#================ +# Structural cols +# X_structural_FN +#================ +feature_gp_nameSTR = 'structural' +n_featuresSTR = len(X_structural_FN) -#-------- -# CV -#-------- -baseline_CT = baseline_all_gn.filter(regex = '.*_time|test_.*|TN|FP|FN|TP|.*_neg|.*_pos', axis = 0) -baseline_CT = baseline_CT.reset_index() -baseline_CT.rename(columns = {'index': 'original_names'}, inplace = True) -baseline_CT['score_type'] = baseline_CT['original_names'] -baseline_CT['score_type'] = baseline_CT['score_type'].str.replace('test_*', '', regex = True) -baseline_CT['data_source'] = 'CT_score' +scores_mmSTR = MultModelsCl_dissected(input_df = X[X_structural_FN] + , target = y + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts[X_structural_FN] + , blind_test_target = y_bts + , add_cm = True + , add_yn = True) -#---------------------- -# rpow bind: CT and BT -#---------------------- -if all(baseline_BT.columns == baseline_CT.columns): - print('\nPASS: Colnames match, proceeding to row bind for data:', feature_gp_name - , '\nDim of df1 (BT):', baseline_BT.shape - , '\nDim of df2 (CT):', baseline_CT.shape) - comb_df_gn = pd.concat([baseline_BT, baseline_CT], axis = 0, ignore_index = True) - comb_df_gn['feature_group'] = feature_gp_name - print('\nDim of combined df:', comb_df_gn.shape) +baseline_allSTR = pd.DataFrame(scores_mmSTR) + +baseline_STR = baseline_allSTR.filter(regex = 'bts_.*|test_.*|.*_time|TN|FP|FN|TP|.*_neg|.*_pos', axis = 0) +baseline_STR = baseline_STR.reset_index() +baseline_STR.rename(columns = {'index': 'original_names'}, inplace = True) + +# Indicate whether BT or CT +bt_pattern = re.compile(r'bts_.*') +baseline_STR['data_source'] = baseline_STR.apply(lambda row: 'BT' if bt_pattern.search(row.original_names) else 'CV' , axis = 1) + +baseline_STR['score_type'] = baseline_STR['original_names'].str.replace('bts_|test_', '', regex = True) + +score_type_uniqueN = set(baseline_STR['score_type']) +cL1 = list(score_type_ordermapD.keys()) +cL2 = list(score_type_uniqueN) + +if set(cL1).issubset(cL2): + print('\nPASS: sorting df by score that is mapped onto the order I want') + baseline_STR['score_order'] = baseline_STR['score_type'].map(score_type_ordermapD) + baseline_STR.sort_values(by = ['data_source', 'score_order'], ascending = [True, True], inplace = True) else: - print('\nFAIL: colnames mismatch, cannot combine') + sys.exit('\nFAIL: could not sort df as score mapping for ordering failed') -# good way but I don't like to have to rearrange the columns later -#frames_tocombine = [baseline_BT, baseline_CT] -#common_cols = list(set.intersection(*(set(df.columns) for df in frames_tocombine))) -#a = pd.concat([df[common_cols] for df in frames_tocombine], ignore_index=True) +baseline_STR['feature_group'] = feature_gp_nameSTR +baseline_STR['sampling_type'] = sampling_type_name +baseline_STR['tts_split'] = tts_split_name +baseline_STR['n_features'] = n_featuresSTR +############################################################################## +#================ +# Stability cols +# X_stability_FN +#================ +feature_gp_nameSTB = 'stability' +n_featuresSTB = len(X_stability_FN) + +scores_mmSTB = MultModelsCl_dissected(input_df = X[X_stability_FN] + , target = y + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts[X_stability_FN] + , blind_test_target = y_bts + , add_cm = True + , add_yn = True) + +baseline_allSTB = pd.DataFrame(scores_mmSTB) + +baseline_STB = baseline_allSTB.filter(regex = 'bts_.*|test_.*|.*_time|TN|FP|FN|TP|.*_neg|.*_pos', axis = 0) +baseline_STB = baseline_STB.reset_index() +baseline_STB.rename(columns = {'index': 'original_names'}, inplace = True) + +# Indicate whether BT or CT +bt_pattern = re.compile(r'bts_.*') +baseline_STB['data_source'] = baseline_STB.apply(lambda row: 'BT' if bt_pattern.search(row.original_names) else 'CV' , axis = 1) + +baseline_STB['score_type'] = baseline_STB['original_names'].str.replace('bts_|test_', '', regex = True) + +score_type_uniqueN = set(baseline_STB['score_type']) +cL1 = list(score_type_ordermapD.keys()) +cL2 = list(score_type_uniqueN) + +if set(cL1).issubset(cL2): + print('\nPASS: sorting df by score that is mapped onto the order I want') + baseline_STB['score_order'] = baseline_STB['score_type'].map(score_type_ordermapD) + baseline_STB.sort_values(by = ['data_source', 'score_order'], ascending = [True, True], inplace = True) +else: + sys.exit('\nFAIL: could not sort df as score mapping for ordering failed') + +baseline_STB['feature_group'] = feature_gp_nameSTB +baseline_STB['sampling_type'] = sampling_type_name +baseline_STB['tts_split'] = tts_split_name +baseline_STB['n_features'] = n_featuresSTB ############################################################################### #================ -# Evolution +# Affinity cols +# X_affinityFN #================ +feature_gp_nameAFF = 'affinity' +n_featuresAFF = len(X_affinityFN) +scores_mmAFF = MultModelsCl_dissected(input_df = X[X_affinityFN] + , target = y + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts[X_affinityFN] + , blind_test_target = y_bts + , add_cm = True + , add_yn = True) + +baseline_allAFF = pd.DataFrame(scores_mmAFF) + +baseline_AFF = baseline_allAFF.filter(regex = 'bts_.*|test_.*|.*_time|TN|FP|FN|TP|.*_neg|.*_pos', axis = 0) +baseline_AFF = baseline_AFF.reset_index() +baseline_AFF.rename(columns = {'index': 'original_names'}, inplace = True) + +# Indicate whether BT or CT +bt_pattern = re.compile(r'bts_.*') +baseline_AFF['data_source'] = baseline_AFF.apply(lambda row: 'BT' if bt_pattern.search(row.original_names) else 'CV' , axis = 1) + +baseline_AFF['score_type'] = baseline_AFF['original_names'].str.replace('bts_|test_', '', regex = True) + +score_type_uniqueN = set(baseline_AFF['score_type']) +cL1 = list(score_type_ordermapD.keys()) +cL2 = list(score_type_uniqueN) + +if set(cL1).issubset(cL2): + print('\nPASS: sorting df by score that is mapped onto the order I want') + baseline_AFF['score_order'] = baseline_AFF['score_type'].map(score_type_ordermapD) + baseline_AFF.sort_values(by = ['data_source', 'score_order'], ascending = [True, True], inplace = True) +else: + sys.exit('\nFAIL: could not sort df as score mapping for ordering failed') + +baseline_AFF['feature_group'] = feature_gp_nameAFF +baseline_AFF['sampling_type'] = sampling_type_name +baseline_AFF['tts_split'] = tts_split_name +baseline_AFF['n_features'] = n_featuresAFF +############################################################################### +#================ +# Residue level +# X_resprop_FN +#================ +feature_gp_nameRES = 'residue_prop' +n_featuresRES = len(X_resprop_FN) + +scores_mmRES = MultModelsCl_dissected(input_df = X[X_resprop_FN] + , target = y + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts[X_resprop_FN] + , blind_test_target = y_bts + , add_cm = True + , add_yn = True) + +baseline_allRES = pd.DataFrame(scores_mmRES) + +baseline_RES = baseline_allRES.filter(regex = 'bts_.*|test_.*|.*_time|TN|FP|FN|TP|.*_neg|.*_pos', axis = 0) +baseline_RES = baseline_RES.reset_index() +baseline_RES.rename(columns = {'index': 'original_names'}, inplace = True) + +# Indicate whether BT or CT +bt_pattern = re.compile(r'bts_.*') +baseline_RES['data_source'] = baseline_RES.apply(lambda row: 'BT' if bt_pattern.search(row.original_names) else 'CV' , axis = 1) + +baseline_RES['score_type'] = baseline_RES['original_names'].str.replace('bts_|test_', '', regex = True) + +score_type_uniqueN = set(baseline_RES['score_type']) +cL1 = list(score_type_ordermapD.keys()) +cL2 = list(score_type_uniqueN) + +if set(cL1).issubset(cL2): + print('\nPASS: sorting df by score that is mapped onto the order I want') + baseline_RES['score_order'] = baseline_RES['score_type'].map(score_type_ordermapD) + baseline_RES.sort_values(by = ['data_source', 'score_order'], ascending = [True, True], inplace = True) +else: + sys.exit('\nFAIL: could not sort df as score mapping for ordering failed') + +baseline_RES['feature_group'] = feature_gp_nameRES +baseline_RES['sampling_type'] = sampling_type_name +baseline_RES['tts_split'] = tts_split_name +baseline_RES['n_features'] = n_featuresRES +############################################################################### +#================ +# Residue level-AAindex +#X_resprop_FN - X_aaindex_Fnum +#================ +X_respropNOaaFN = list(set(X_resprop_FN) - set(X_aaindex_Fnum)) + +feature_gp_nameRNAA = 'ResPropNoAA' +n_featuresRNAA = len(X_respropNOaaFN) + +scores_mmRNAA = MultModelsCl_dissected(input_df = X[X_respropNOaaFN] + , target = y + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts[X_respropNOaaFN] + , blind_test_target = y_bts + , add_cm = True + , add_yn = True) + +baseline_allRNAA = pd.DataFrame(scores_mmRNAA) + +baseline_RNAA = baseline_allRNAA.filter(regex = 'bts_.*|test_.*|.*_time|TN|FP|FN|TP|.*_neg|.*_pos', axis = 0) +baseline_RNAA = baseline_RNAA.reset_index() +baseline_RNAA.rename(columns = {'index': 'original_names'}, inplace = True) + +# Indicate whether BT or CT +bt_pattern = re.compile(r'bts_.*') +baseline_RNAA['data_source'] = baseline_RNAA.apply(lambda row: 'BT' if bt_pattern.search(row.original_names) else 'CV' , axis = 1) + +baseline_RNAA['score_type'] = baseline_RNAA['original_names'].str.replace('bts_|test_', '', regex = True) + +score_type_uniqueN = set(baseline_RNAA['score_type']) +cL1 = list(score_type_ordermapD.keys()) +cL2 = list(score_type_uniqueN) + +if set(cL1).issubset(cL2): + print('\nPASS: sorting df by score that is mapped onto the order I want') + baseline_RNAA['score_order'] = baseline_RNAA['score_type'].map(score_type_ordermapD) + baseline_RNAA.sort_values(by = ['data_source', 'score_order'], ascending = [True, True], inplace = True) +else: + sys.exit('\nFAIL: could not sort df as score mapping for ordering failed') + +baseline_RNAA['feature_group'] = feature_gp_nameRNAA +baseline_RNAA['sampling_type'] = sampling_type_name +baseline_RNAA['tts_split'] = tts_split_name +baseline_RNAA['n_features'] = n_featuresRNAA +############################################################################### +#================ +# Structural cols-AAindex +#X_structural_FN - X_aaindex_Fnum +#================ +X_strNOaaFN = list(set(X_structural_FN) - set(X_aaindex_Fnum)) + +feature_gp_nameSNAA = 'StrNoAA' +n_featuresSNAA = len(X_strNOaaFN) + +scores_mmSNAA = MultModelsCl_dissected(input_df = X[X_strNOaaFN] + , target = y + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts[X_strNOaaFN] + , blind_test_target = y_bts + , add_cm = True + , add_yn = True) + +baseline_allSNAA = pd.DataFrame(scores_mmSNAA) + +baseline_SNAA = baseline_allSNAA.filter(regex = 'bts_.*|test_.*|.*_time|TN|FP|FN|TP|.*_neg|.*_pos', axis = 0) +baseline_SNAA = baseline_SNAA.reset_index() +baseline_SNAA.rename(columns = {'index': 'original_names'}, inplace = True) + +# Indicate whether BT or CT +bt_pattern = re.compile(r'bts_.*') +baseline_SNAA['data_source'] = baseline_SNAA.apply(lambda row: 'BT' if bt_pattern.search(row.original_names) else 'CV' , axis = 1) + +baseline_SNAA['score_type'] = baseline_SNAA['original_names'].str.replace('bts_|test_', '', regex = True) + +score_type_uniqueN = set(baseline_SNAA['score_type']) +cL1 = list(score_type_ordermapD.keys()) +cL2 = list(score_type_uniqueN) + +if set(cL1).issubset(cL2): + print('\nPASS: sorting df by score that is mapped onto the order I want') + baseline_SNAA['score_order'] = baseline_SNAA['score_type'].map(score_type_ordermapD) + baseline_SNAA.sort_values(by = ['data_source', 'score_order'], ascending = [True, True], inplace = True) +else: + sys.exit('\nFAIL: could not sort df as score mapping for ordering failed') + +baseline_SNAA['feature_group'] = feature_gp_nameSNAA +baseline_SNAA['sampling_type'] = sampling_type_name +baseline_SNAA['tts_split'] = tts_split_name +baseline_SNAA['n_features'] = n_featuresSNAA +############################################################################### +#%% COMBINING all FG dfs +#================ +# Combine all +# https://stackoverflow.com/questions/39862654/pandas-concat-of-multiple-data-frames-using-only-common-columns +#================ +dfs_combine = [baseline_EV, baseline_GN, baseline_STR, baseline_STB, baseline_AFF, baseline_RES , baseline_RNAA , baseline_SNAA] + +dfs_nrows = [] +for df in dfs_combine: + dfs_nrows = dfs_nrows + [len(df)] +dfs_nrows = max(dfs_nrows) + +dfs_ncols = [] +for df in dfs_combine: + dfs_ncols = dfs_ncols + [len(df.columns)] +dfs_ncols = max(dfs_ncols) + +# dfs_ncols = [] +# dfs_ncols2 = mode(dfs_ncols.append(len(df.columns) for df in dfs_combine) +# dfs_ncols2 + +expected_nrows = len(dfs_combine) * dfs_nrows +expected_ncols = dfs_ncols + +common_cols = list(set.intersection(*(set(df.columns) for df in dfs_combine))) + +if len(common_cols) == dfs_ncols : + combined_FG_baseline = pd.concat([df[common_cols] for df in dfs_combine], ignore_index=True) + fgs = combined_FG_baseline[['feature_group', 'n_features']] + fgs = fgs.drop_duplicates() + print('\nConcatenating dfs with feature groups after ML analysis (sampling type):' + , '\nNo. of dfs combining:', len(dfs_combine) + , '\nSampling type:', sampling_type + , '\nThe feature groups are:' + , '\n', fgs) + if len(combined_FG_baseline) == expected_nrows and len(combined_FG_baseline.columns) == expected_ncols: + print('\nPASS:', len(dfs_combine), 'dfs successfully combined' + , '\nnrows in combined_df:', len(combined_FG_baseline) + , '\nncols in combined_df:', len(combined_FG_baseline.columns)) + else: + print('\nFAIL: concatenating failed' + , '\nExpected nrows:', expected_nrows + , '\nGot:', len(combined_FG_baseline) + , '\nExpected ncols:', expected_ncols + , '\nGot:', len(combined_FG_baseline.columns)) + sys.exit() +else: + sys.exit('\nConcatenting dfs not possible,check numbers ') + +# # rpow bind +# if all(ll((baseline_EV.columns == baseline_GN.columns == baseline_STR.columns)): +# print('\nPASS:colnames match, proceeding to rowbind') +# comb_df = pd.concat()], axis = 0, ignore_index = True ) +############################################################################### +#==================== +# Write output file +#==================== + +combined_FG_baseline.to_csv(outFile) +print('\nFile successfully written:', outFile) +############################################################################### \ No newline at end of file diff --git a/scripts/ml/ml_data_dissected.py b/scripts/ml/ml_data_dissected.py index a589449..218cd30 100644 --- a/scripts/ml/ml_data_dissected.py +++ b/scripts/ml/ml_data_dissected.py @@ -5,785 +5,787 @@ 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)) +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 -#--------------- -# 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') + from imblearn.over_sampling import SMOTENC + from imblearn.under_sampling import EditedNearestNeighbours + from imblearn.under_sampling import RepeatedEditedNearestNeighbours -#--------------------------- -# 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)) + 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 -# 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!') + from sklearn.model_selection import train_test_split, cross_validate, cross_val_score + from sklearn.model_selection import StratifiedKFold,RepeatedStratifiedKFold, RepeatedKFold -#-------- -# 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'] ) - + from sklearn.pipeline import Pipeline, make_pipeline + import argparse + import re + #%% GLOBALS + rs = {'random_state': 42} + njobs = {'n_jobs': 10} -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 a common one - 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' - , 'mmcsm_lig' - , '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'] - -# 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' + 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) + }) - , '\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 + 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 a common one + 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' + , 'mmcsm_lig' + , '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'] + + # 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 diff --git a/scripts/ml/pnca_config_dissected.py b/scripts/ml/pnca_config_dissected.py deleted file mode 100644 index dafaff2..0000000 --- a/scripts/ml/pnca_config_dissected.py +++ /dev/null @@ -1,207 +0,0 @@ -#!/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('\n================================================================\n') - - -X_structural_FN = X_stability_FN + X_affinityFN + X_resprop_FN - X_aaindex_Fnum + X_str_Fnum + X_aap_Fcat -all_featuresN = X_evolFN + X_structural_FN + X_genomicFN - -print('\n================================================================' - - , '\nTotal Evolutionary features (n):' , len(X_evolFN) - , '\n--------------Evol. feature colnames:', X_evolFN - - , '\n================================================================' - - , '\n\nTotal structural features (n):', len(X_structural_FN) - - , '\n--------Stability ncols:' , len(X_stability_FN) - , '\n--------------Common stability colnames:' , X_common_stability_Fnum - , '\n--------------Foldx colnames:' , X_foldX_Fnum - - , '\n--------Affinity ncols:' , len(X_affinityFN) - , '\n--------------Common affinity colnames:' , common_affinity_Fnum - , '\n--------------Gene specific affinity colnames:', gene_affinity_colnames - - , '\n--------Residue prop ncols:' , len(X_resprop_FN) - , '\n--------------Residue Prop cols:' , X_str_Fnum - , '\n--------------AA change Prop cols:' , X_aap_Fcat - , '\n--------------AA index cols:' , X_aaindex_Fnum - - , '\n================================================================' - - , '\n\nTotal Genomic features (n):' , len(X_genomicFN) - , '\n--------MAF+OR cols:' , len(X_gn_mafor_Fnum) - , '\n--------------MAF+OR colnames:' , X_gn_mafor_Fnum - - , '\n--------Lineage cols:' , len(X_gn_linegae_Fnum) - , '\n--------------Lineage cols:' , X_gn_linegae_Fnum - - , '\n--------Other cols:' , len(X_gn_Fcat) - , '\n--------------Other cols:' , X_gn_Fcat - - , '\n================================================================') - -# Sanity check -if ( len(X.columns) == len(X_evolFN) + len(X_structural_FN) + len(X_genomicFN)): - print('\nPass: No. of features match') -else: - print('\nFail: Count of feature mismatch' - , '\nExpected:', len(X_evolFN) + len(X_structural_FN) + len(X_genomicFN) - , '\nGot:', len(X.columns)) - sys.exit() - -print('\n#####################################################################\n') - -# ############################################################################### -# #================== -# # Baseline models -# #================== -# mm_skf_scoresD = MultModelsCl_dissected(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')