diff --git a/scripts/ml/alr_8020.py b/scripts/ml/alr_8020.py new file mode 100755 index 0000000..3309f04 --- /dev/null +++ b/scripts/ml/alr_8020.py @@ -0,0 +1,203 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat May 28 05:25:30 2022 + +@author: tanu +""" + +import os + +gene = 'alr' +drug = 'cycloserine' +#total_mtblineage_uc = 8 + +homedir = os.path.expanduser("~") +os.chdir( homedir + '/git/LSHTM_analysis/scripts/ml/') + +from ml_data_8020 import * +setvars(gene,drug) +from ml_data_8020 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 import MultModelsCl + +################################################################################ +print('\n#####################################################################\n' + , '\nRunning ML analysis: 80/20 split' + , '\nGene name:', gene + , '\nDrug name:', drug) + +#================== +# Specify outdir +#================== + +outdir_ml = outdir + 'ml/tts_8020/' + +print('\nOutput directory:', outdir_ml) + +#%%########################################################################### +print('\nSanity checks:' + #, '\nML source data size:', x_features.shape + , '\nTotal input features:', len(X.columns) + , '\n' + , '\nTraining data size:', X.shape + , '\nTest data size:', X_bts.shape + , '\n' + , '\nTarget feature numbers (training data):', Counter(y) + , '\nTarget features ratio (training data:', yc1_ratio + , '\n' + , '\nTarget feature numbers (test data):', Counter(y_bts) + , '\nTarget features ratio (test data):', yc2_ratio + + , '\n\n#####################################################################\n') + +print('\n================================================================\n') + +print('Strucutral features (n):' + , len(X_ssFN) + , '\nThese are:' + , '\nCommon stablity features:', X_stabilityN + , '\nFoldX columns:', X_foldX_cols + , '\nOther struc columns:', X_str + , '\n================================================================\n') + +print('AAindex features (n):' + , len(X_aaindexFN) + , '\nThese are:\n' + , X_aaindexFN + , '\n================================================================\n') + +print('Evolutionary features (n):' + , len(X_evolFN) + , '\nThese are:\n' + , X_evolFN + , '\n================================================================\n') + +print('Genomic features (n):' + , len(X_genomicFN) + , '\nThese are:\n' + , X_genomic_mafor, '\n' + , X_genomic_linegae + , '\n================================================================\n') + +print('Categorical features (n):' + , len(categorical_FN) + , '\nThese are:\n' + , categorical_FN + , '\n================================================================\n') + +if ( len(X.columns) == len(X_ssFN) + len(X_aaindexFN) + len(X_evolFN) + len(X_genomicFN) + len(categorical_FN) ): + print('\nPass: No. of features match') +else: + sys.exit('\nFail: Count of feature mismatch') + +print('\n#####################################################################\n') + +################################################################################ +#================== +# Baseline models +#================== +mm_skf_scoresD = MultModelsCl(input_df = X + , target = y + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) + +baseline_all = pd.DataFrame(mm_skf_scoresD) +baseline_all = baseline_all.T +#baseline_train = baseline_all.filter(like='train_', axis=1) +baseline_CT = baseline_all.filter(like='test_', axis=1) +baseline_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +baseline_BT = baseline_all.filter(like='bts_', axis=1) +baseline_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +baseline_CT.to_csv(outdir_ml + gene.lower() + '_baseline_CT_allF.csv') +baseline_BT.to_csv(outdir_ml + gene.lower() + '_baseline_BT_allF.csv') + + +#%% SMOTE NC: Oversampling [Numerical + categorical] +mm_skf_scoresD7 = MultModelsCl(input_df = X_smnc + , target = y_smnc + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +smnc_all = pd.DataFrame(mm_skf_scoresD7) +smnc_all = smnc_all.T + +smnc_CT = smnc_all.filter(like='test_', axis=1) +smnc_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +smnc_BT = smnc_all.filter(like='bts_', axis=1) +smnc_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +smnc_CT.to_csv(outdir_ml + gene.lower() + '_smnc_CT_allF.csv') +smnc_BT.to_csv(outdir_ml + gene.lower() + '_smnc_BT_allF.csv') + +#%% ROS: Numerical + categorical +mm_skf_scoresD3 = MultModelsCl(input_df = X_ros + , target = y_ros + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +ros_all = pd.DataFrame(mm_skf_scoresD3) +ros_all = ros_all.T + +ros_CT = ros_all.filter(like='test_', axis=1) +ros_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +ros_BT = ros_all.filter(like='bts_', axis=1) +ros_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +ros_CT.to_csv(outdir_ml + gene.lower() + '_ros_CT_allF.csv') +ros_BT.to_csv(outdir_ml + gene.lower() + '_ros_BT_allF.csv') + +#%% RUS: Numerical + categorical +mm_skf_scoresD4 = MultModelsCl(input_df = X_rus + , target = y_rus + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +rus_all = pd.DataFrame(mm_skf_scoresD4) +rus_all = rus_all.T + +rus_CT = rus_all.filter(like='test_', axis=1) +rus_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +rus_BT = rus_all.filter(like='bts_' , axis=1) +rus_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +rus_CT.to_csv(outdir_ml + gene.lower() + '_rus_CT_allF.csv') +rus_BT.to_csv(outdir_ml + gene.lower() + '_rus_BT_allF.csv') + +#%% ROS + RUS Combined: Numerical + categorical +mm_skf_scoresD8 = MultModelsCl(input_df = X_rouC + , target = y_rouC + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +rouC_all = pd.DataFrame(mm_skf_scoresD8) +rouC_all = rouC_all.T + +rouC_CT = rouC_all.filter(like='test_', axis=1) +rouC_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +rouC_BT = rouC_all.filter(like='bts_', axis=1) +rouC_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +rouC_CT.to_csv(outdir_ml + gene.lower() + '_rouC_CT_allF.csv') +rouC_BT.to_csv(outdir_ml + gene.lower() + '_rouC_BT_allF.csv') diff --git a/scripts/ml/embb_8020.py b/scripts/ml/embb_8020.py new file mode 100755 index 0000000..5e8dd31 --- /dev/null +++ b/scripts/ml/embb_8020.py @@ -0,0 +1,203 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat May 28 05:25:30 2022 + +@author: tanu +""" + +import os + +gene = 'embB' +drug = 'ethambutol' +#total_mtblineage_uc = 8 + +homedir = os.path.expanduser("~") +os.chdir( homedir + '/git/LSHTM_analysis/scripts/ml/') + +from ml_data_8020 import * +setvars(gene,drug) +from ml_data_8020 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 import MultModelsCl + +################################################################################ +print('\n#####################################################################\n' + , '\nRunning ML analysis: 80/20 split' + , '\nGene name:', gene + , '\nDrug name:', drug) + +#================== +# Specify outdir +#================== + +outdir_ml = outdir + 'ml/tts_8020/' + +print('\nOutput directory:', outdir_ml) + +#%%########################################################################### +print('\nSanity checks:' + #, '\nML source data size:', x_features.shape + , '\nTotal input features:', len(X.columns) + , '\n' + , '\nTraining data size:', X.shape + , '\nTest data size:', X_bts.shape + , '\n' + , '\nTarget feature numbers (training data):', Counter(y) + , '\nTarget features ratio (training data:', yc1_ratio + , '\n' + , '\nTarget feature numbers (test data):', Counter(y_bts) + , '\nTarget features ratio (test data):', yc2_ratio + + , '\n\n#####################################################################\n') + +print('\n================================================================\n') + +print('Strucutral features (n):' + , len(X_ssFN) + , '\nThese are:' + , '\nCommon stablity features:', X_stabilityN + , '\nFoldX columns:', X_foldX_cols + , '\nOther struc columns:', X_str + , '\n================================================================\n') + +print('AAindex features (n):' + , len(X_aaindexFN) + , '\nThese are:\n' + , X_aaindexFN + , '\n================================================================\n') + +print('Evolutionary features (n):' + , len(X_evolFN) + , '\nThese are:\n' + , X_evolFN + , '\n================================================================\n') + +print('Genomic features (n):' + , len(X_genomicFN) + , '\nThese are:\n' + , X_genomic_mafor, '\n' + , X_genomic_linegae + , '\n================================================================\n') + +print('Categorical features (n):' + , len(categorical_FN) + , '\nThese are:\n' + , categorical_FN + , '\n================================================================\n') + +if ( len(X.columns) == len(X_ssFN) + len(X_aaindexFN) + len(X_evolFN) + len(X_genomicFN) + len(categorical_FN) ): + print('\nPass: No. of features match') +else: + sys.exit('\nFail: Count of feature mismatch') + +print('\n#####################################################################\n') + +################################################################################ +#================== +# Baseline models +#================== +mm_skf_scoresD = MultModelsCl(input_df = X + , target = y + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) + +baseline_all = pd.DataFrame(mm_skf_scoresD) +baseline_all = baseline_all.T +#baseline_train = baseline_all.filter(like='train_', axis=1) +baseline_CT = baseline_all.filter(like='test_', axis=1) +baseline_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +baseline_BT = baseline_all.filter(like='bts_', axis=1) +baseline_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +baseline_CT.to_csv(outdir_ml + gene.lower() + '_baseline_CT_allF.csv') +baseline_BT.to_csv(outdir_ml + gene.lower() + '_baseline_BT_allF.csv') + + +#%% SMOTE NC: Oversampling [Numerical + categorical] +mm_skf_scoresD7 = MultModelsCl(input_df = X_smnc + , target = y_smnc + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +smnc_all = pd.DataFrame(mm_skf_scoresD7) +smnc_all = smnc_all.T + +smnc_CT = smnc_all.filter(like='test_', axis=1) +smnc_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +smnc_BT = smnc_all.filter(like='bts_', axis=1) +smnc_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +smnc_CT.to_csv(outdir_ml + gene.lower() + '_smnc_CT_allF.csv') +smnc_BT.to_csv(outdir_ml + gene.lower() + '_smnc_BT_allF.csv') + +#%% ROS: Numerical + categorical +mm_skf_scoresD3 = MultModelsCl(input_df = X_ros + , target = y_ros + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +ros_all = pd.DataFrame(mm_skf_scoresD3) +ros_all = ros_all.T + +ros_CT = ros_all.filter(like='test_', axis=1) +ros_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +ros_BT = ros_all.filter(like='bts_', axis=1) +ros_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +ros_CT.to_csv(outdir_ml + gene.lower() + '_ros_CT_allF.csv') +ros_BT.to_csv(outdir_ml + gene.lower() + '_ros_BT_allF.csv') + +#%% RUS: Numerical + categorical +mm_skf_scoresD4 = MultModelsCl(input_df = X_rus + , target = y_rus + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +rus_all = pd.DataFrame(mm_skf_scoresD4) +rus_all = rus_all.T + +rus_CT = rus_all.filter(like='test_', axis=1) +rus_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +rus_BT = rus_all.filter(like='bts_' , axis=1) +rus_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +rus_CT.to_csv(outdir_ml + gene.lower() + '_rus_CT_allF.csv') +rus_BT.to_csv(outdir_ml + gene.lower() + '_rus_BT_allF.csv') + +#%% ROS + RUS Combined: Numerical + categorical +mm_skf_scoresD8 = MultModelsCl(input_df = X_rouC + , target = y_rouC + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +rouC_all = pd.DataFrame(mm_skf_scoresD8) +rouC_all = rouC_all.T + +rouC_CT = rouC_all.filter(like='test_', axis=1) +rouC_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +rouC_BT = rouC_all.filter(like='bts_', axis=1) +rouC_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +rouC_CT.to_csv(outdir_ml + gene.lower() + '_rouC_CT_allF.csv') +rouC_BT.to_csv(outdir_ml + gene.lower() + '_rouC_BT_allF.csv') diff --git a/scripts/ml/gid_8020.py b/scripts/ml/gid_8020.py new file mode 100755 index 0000000..9366c50 --- /dev/null +++ b/scripts/ml/gid_8020.py @@ -0,0 +1,203 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat May 28 05:25:30 2022 + +@author: tanu +""" + +import os + +gene = 'gid' +drug = 'streptomycin' +#total_mtblineage_uc = 8 + +homedir = os.path.expanduser("~") +os.chdir( homedir + '/git/LSHTM_analysis/scripts/ml/') + +from ml_data_8020 import * +setvars(gene,drug) +from ml_data_8020 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 import MultModelsCl + +################################################################################ +print('\n#####################################################################\n' + , '\nRunning ML analysis: 80/20 split' + , '\nGene name:', gene + , '\nDrug name:', drug) + +#================== +# Specify outdir +#================== + +outdir_ml = outdir + 'ml/tts_8020/' + +print('\nOutput directory:', outdir_ml) + +#%%########################################################################### +print('\nSanity checks:' + #, '\nML source data size:', x_features.shape + , '\nTotal input features:', len(X.columns) + , '\n' + , '\nTraining data size:', X.shape + , '\nTest data size:', X_bts.shape + , '\n' + , '\nTarget feature numbers (training data):', Counter(y) + , '\nTarget features ratio (training data:', yc1_ratio + , '\n' + , '\nTarget feature numbers (test data):', Counter(y_bts) + , '\nTarget features ratio (test data):', yc2_ratio + + , '\n\n#####################################################################\n') + +print('\n================================================================\n') + +print('Strucutral features (n):' + , len(X_ssFN) + , '\nThese are:' + , '\nCommon stablity features:', X_stabilityN + , '\nFoldX columns:', X_foldX_cols + , '\nOther struc columns:', X_str + , '\n================================================================\n') + +print('AAindex features (n):' + , len(X_aaindexFN) + , '\nThese are:\n' + , X_aaindexFN + , '\n================================================================\n') + +print('Evolutionary features (n):' + , len(X_evolFN) + , '\nThese are:\n' + , X_evolFN + , '\n================================================================\n') + +print('Genomic features (n):' + , len(X_genomicFN) + , '\nThese are:\n' + , X_genomic_mafor, '\n' + , X_genomic_linegae + , '\n================================================================\n') + +print('Categorical features (n):' + , len(categorical_FN) + , '\nThese are:\n' + , categorical_FN + , '\n================================================================\n') + +if ( len(X.columns) == len(X_ssFN) + len(X_aaindexFN) + len(X_evolFN) + len(X_genomicFN) + len(categorical_FN) ): + print('\nPass: No. of features match') +else: + sys.exit('\nFail: Count of feature mismatch') + +print('\n#####################################################################\n') + +################################################################################ +#================== +# Baseline models +#================== +mm_skf_scoresD = MultModelsCl(input_df = X + , target = y + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) + +baseline_all = pd.DataFrame(mm_skf_scoresD) +baseline_all = baseline_all.T +#baseline_train = baseline_all.filter(like='train_', axis=1) +baseline_CT = baseline_all.filter(like='test_', axis=1) +baseline_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +baseline_BT = baseline_all.filter(like='bts_', axis=1) +baseline_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +baseline_CT.to_csv(outdir_ml + gene.lower() + '_baseline_CT_allF.csv') +baseline_BT.to_csv(outdir_ml + gene.lower() + '_baseline_BT_allF.csv') + + +#%% SMOTE NC: Oversampling [Numerical + categorical] +mm_skf_scoresD7 = MultModelsCl(input_df = X_smnc + , target = y_smnc + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +smnc_all = pd.DataFrame(mm_skf_scoresD7) +smnc_all = smnc_all.T + +smnc_CT = smnc_all.filter(like='test_', axis=1) +smnc_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +smnc_BT = smnc_all.filter(like='bts_', axis=1) +smnc_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +smnc_CT.to_csv(outdir_ml + gene.lower() + '_smnc_CT_allF.csv') +smnc_BT.to_csv(outdir_ml + gene.lower() + '_smnc_BT_allF.csv') + +#%% ROS: Numerical + categorical +mm_skf_scoresD3 = MultModelsCl(input_df = X_ros + , target = y_ros + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +ros_all = pd.DataFrame(mm_skf_scoresD3) +ros_all = ros_all.T + +ros_CT = ros_all.filter(like='test_', axis=1) +ros_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +ros_BT = ros_all.filter(like='bts_', axis=1) +ros_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +ros_CT.to_csv(outdir_ml + gene.lower() + '_ros_CT_allF.csv') +ros_BT.to_csv(outdir_ml + gene.lower() + '_ros_BT_allF.csv') + +#%% RUS: Numerical + categorical +mm_skf_scoresD4 = MultModelsCl(input_df = X_rus + , target = y_rus + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +rus_all = pd.DataFrame(mm_skf_scoresD4) +rus_all = rus_all.T + +rus_CT = rus_all.filter(like='test_', axis=1) +rus_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +rus_BT = rus_all.filter(like='bts_' , axis=1) +rus_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +rus_CT.to_csv(outdir_ml + gene.lower() + '_rus_CT_allF.csv') +rus_BT.to_csv(outdir_ml + gene.lower() + '_rus_BT_allF.csv') + +#%% ROS + RUS Combined: Numerical + categorical +mm_skf_scoresD8 = MultModelsCl(input_df = X_rouC + , target = y_rouC + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +rouC_all = pd.DataFrame(mm_skf_scoresD8) +rouC_all = rouC_all.T + +rouC_CT = rouC_all.filter(like='test_', axis=1) +rouC_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +rouC_BT = rouC_all.filter(like='bts_', axis=1) +rouC_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +rouC_CT.to_csv(outdir_ml + gene.lower() + '_rouC_CT_allF.csv') +rouC_BT.to_csv(outdir_ml + gene.lower() + '_rouC_BT_allF.csv') diff --git a/scripts/ml/katg_8020.py b/scripts/ml/katg_8020.py new file mode 100755 index 0000000..20be7d9 --- /dev/null +++ b/scripts/ml/katg_8020.py @@ -0,0 +1,203 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat May 28 05:25:30 2022 + +@author: tanu +""" + +import os + +gene = 'katG' +drug = 'isoniazid' +#total_mtblineage_uc = 8 + +homedir = os.path.expanduser("~") +os.chdir( homedir + '/git/LSHTM_analysis/scripts/ml/') + +from ml_data_8020 import * +setvars(gene,drug) +from ml_data_8020 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 import MultModelsCl + +################################################################################ +print('\n#####################################################################\n' + , '\nRunning ML analysis: 80/20 split' + , '\nGene name:', gene + , '\nDrug name:', drug) + +#================== +# Specify outdir +#================== + +outdir_ml = outdir + 'ml/tts_8020/' + +print('\nOutput directory:', outdir_ml) + +#%%########################################################################### +print('\nSanity checks:' + #, '\nML source data size:', x_features.shape + , '\nTotal input features:', len(X.columns) + , '\n' + , '\nTraining data size:', X.shape + , '\nTest data size:', X_bts.shape + , '\n' + , '\nTarget feature numbers (training data):', Counter(y) + , '\nTarget features ratio (training data:', yc1_ratio + , '\n' + , '\nTarget feature numbers (test data):', Counter(y_bts) + , '\nTarget features ratio (test data):', yc2_ratio + + , '\n\n#####################################################################\n') + +print('\n================================================================\n') + +print('Strucutral features (n):' + , len(X_ssFN) + , '\nThese are:' + , '\nCommon stablity features:', X_stabilityN + , '\nFoldX columns:', X_foldX_cols + , '\nOther struc columns:', X_str + , '\n================================================================\n') + +print('AAindex features (n):' + , len(X_aaindexFN) + , '\nThese are:\n' + , X_aaindexFN + , '\n================================================================\n') + +print('Evolutionary features (n):' + , len(X_evolFN) + , '\nThese are:\n' + , X_evolFN + , '\n================================================================\n') + +print('Genomic features (n):' + , len(X_genomicFN) + , '\nThese are:\n' + , X_genomic_mafor, '\n' + , X_genomic_linegae + , '\n================================================================\n') + +print('Categorical features (n):' + , len(categorical_FN) + , '\nThese are:\n' + , categorical_FN + , '\n================================================================\n') + +if ( len(X.columns) == len(X_ssFN) + len(X_aaindexFN) + len(X_evolFN) + len(X_genomicFN) + len(categorical_FN) ): + print('\nPass: No. of features match') +else: + sys.exit('\nFail: Count of feature mismatch') + +print('\n#####################################################################\n') + +################################################################################ +#================== +# Baseline models +#================== +mm_skf_scoresD = MultModelsCl(input_df = X + , target = y + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) + +baseline_all = pd.DataFrame(mm_skf_scoresD) +baseline_all = baseline_all.T +#baseline_train = baseline_all.filter(like='train_', axis=1) +baseline_CT = baseline_all.filter(like='test_', axis=1) +baseline_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +baseline_BT = baseline_all.filter(like='bts_', axis=1) +baseline_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +baseline_CT.to_csv(outdir_ml + gene.lower() + '_baseline_CT_allF.csv') +baseline_BT.to_csv(outdir_ml + gene.lower() + '_baseline_BT_allF.csv') + + +#%% SMOTE NC: Oversampling [Numerical + categorical] +mm_skf_scoresD7 = MultModelsCl(input_df = X_smnc + , target = y_smnc + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +smnc_all = pd.DataFrame(mm_skf_scoresD7) +smnc_all = smnc_all.T + +smnc_CT = smnc_all.filter(like='test_', axis=1) +smnc_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +smnc_BT = smnc_all.filter(like='bts_', axis=1) +smnc_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +smnc_CT.to_csv(outdir_ml + gene.lower() + '_smnc_CT_allF.csv') +smnc_BT.to_csv(outdir_ml + gene.lower() + '_smnc_BT_allF.csv') + +#%% ROS: Numerical + categorical +mm_skf_scoresD3 = MultModelsCl(input_df = X_ros + , target = y_ros + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +ros_all = pd.DataFrame(mm_skf_scoresD3) +ros_all = ros_all.T + +ros_CT = ros_all.filter(like='test_', axis=1) +ros_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +ros_BT = ros_all.filter(like='bts_', axis=1) +ros_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +ros_CT.to_csv(outdir_ml + gene.lower() + '_ros_CT_allF.csv') +ros_BT.to_csv(outdir_ml + gene.lower() + '_ros_BT_allF.csv') + +#%% RUS: Numerical + categorical +mm_skf_scoresD4 = MultModelsCl(input_df = X_rus + , target = y_rus + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +rus_all = pd.DataFrame(mm_skf_scoresD4) +rus_all = rus_all.T + +rus_CT = rus_all.filter(like='test_', axis=1) +rus_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +rus_BT = rus_all.filter(like='bts_' , axis=1) +rus_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +rus_CT.to_csv(outdir_ml + gene.lower() + '_rus_CT_allF.csv') +rus_BT.to_csv(outdir_ml + gene.lower() + '_rus_BT_allF.csv') + +#%% ROS + RUS Combined: Numerical + categorical +mm_skf_scoresD8 = MultModelsCl(input_df = X_rouC + , target = y_rouC + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +rouC_all = pd.DataFrame(mm_skf_scoresD8) +rouC_all = rouC_all.T + +rouC_CT = rouC_all.filter(like='test_', axis=1) +rouC_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +rouC_BT = rouC_all.filter(like='bts_', axis=1) +rouC_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +rouC_CT.to_csv(outdir_ml + gene.lower() + '_rouC_CT_allF.csv') +rouC_BT.to_csv(outdir_ml + gene.lower() + '_rouC_BT_allF.csv') diff --git a/scripts/ml/ml_data_8020.py b/scripts/ml/ml_data_8020.py new file mode 100644 index 0000000..9049345 --- /dev/null +++ b/scripts/ml/ml_data_8020.py @@ -0,0 +1,761 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sun Mar 6 13:41:54 2022 + +@author: tanu +""" +def setvars(gene,drug): + #https://stackoverflow.com/questions/51695322/compare-multiple-algorithms-with-sklearn-pipeline + import os, sys + import pandas as pd + import numpy as np + print(np.__version__) + print(pd.__version__) + import pprint as pp + from copy import deepcopy + from collections import Counter + from sklearn.impute import KNNImputer as KNN + from imblearn.over_sampling import RandomOverSampler + from imblearn.under_sampling import RandomUnderSampler + from imblearn.over_sampling import SMOTE + from sklearn.datasets import make_classification + from imblearn.combine import SMOTEENN + from imblearn.combine import SMOTETomek + + from imblearn.over_sampling import SMOTENC + from imblearn.under_sampling import EditedNearestNeighbours + from imblearn.under_sampling import RepeatedEditedNearestNeighbours + + from sklearn.metrics import make_scorer, confusion_matrix, accuracy_score, balanced_accuracy_score, precision_score, average_precision_score, recall_score + from sklearn.metrics import roc_auc_score, roc_curve, f1_score, matthews_corrcoef, jaccard_score, classification_report + + from sklearn.model_selection import train_test_split, cross_validate, cross_val_score + from sklearn.model_selection import StratifiedKFold,RepeatedStratifiedKFold, RepeatedKFold + + from sklearn.pipeline import Pipeline, make_pipeline + #%% GLOBALS + rs = {'random_state': 42} + njobs = {'n_jobs': 10} + + scoring_fn = ({ 'mcc' : make_scorer(matthews_corrcoef) + , 'accuracy' : make_scorer(accuracy_score) + , 'fscore' : make_scorer(f1_score) + , 'precision' : make_scorer(precision_score) + , 'recall' : make_scorer(recall_score) + , 'roc_auc' : make_scorer(roc_auc_score) + , 'jcc' : make_scorer(jaccard_score) + }) + + skf_cv = StratifiedKFold(n_splits = 10 + #, shuffle = False, random_state= None) + , shuffle = True,**rs) + + rskf_cv = RepeatedStratifiedKFold(n_splits = 10 + , n_repeats = 3 + , **rs) + + mcc_score_fn = {'mcc': make_scorer(matthews_corrcoef)} + jacc_score_fn = {'jcc': make_scorer(jaccard_score)} + + #%% FOR LATER: Combine ED logo data + ########################################################################### + rs = {'random_state': 42} + njobs = {'n_jobs': 10} + homedir = os.path.expanduser("~") + + geneL_basic = ['pnca'] + geneL_na = ['gid'] + geneL_na_ppi2 = ['rpob'] + geneL_ppi2 = ['alr', 'embb', 'katg'] + + #num_type = ['int64', 'float64'] + num_type = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64'] + cat_type = ['object', 'bool'] + + #============== + # directories + #============== + datadir = homedir + '/git/Data/' + indir = datadir + drug + '/input/' + outdir = datadir + drug + '/output/' + + #======= + # input + #======= + + #--------- + # File 1 + #--------- + infile_ml1 = outdir + gene.lower() + '_merged_df3.csv' + #infile_ml2 = outdir + gene.lower() + '_merged_df2.csv' + + my_features_df = pd.read_csv(infile_ml1, index_col = 0) + my_features_df = my_features_df .reset_index(drop = True) + my_features_df.index + + my_features_df.dtypes + mycols = my_features_df.columns + + #--------- + # File 2 + #--------- + infile_aaindex = outdir + 'aa_index/' + gene.lower() + '_aa.csv' + aaindex_df = pd.read_csv(infile_aaindex, index_col = 0) + aaindex_df.dtypes + + #----------- + # check for non-numerical columns + #----------- + if any(aaindex_df.dtypes==object): + print('\naaindex_df contains non-numerical data') + + aaindex_df_object = aaindex_df.select_dtypes(include = cat_type) + print('\nTotal no. of non-numerial columns:', len(aaindex_df_object.columns)) + + expected_aa_ncols = len(aaindex_df.columns) - len(aaindex_df_object.columns) + + #----------- + # Extract numerical data only + #----------- + print('\nSelecting numerical data only') + aaindex_df = aaindex_df.select_dtypes(include = num_type) + + #--------------------------- + # aaindex: sanity check 1 + #--------------------------- + if len(aaindex_df.columns) == expected_aa_ncols: + print('\nPASS: successfully selected numerical columns only for aaindex_df') + else: + print('\nFAIL: Numbers mismatch' + , '\nExpected ncols:', expected_aa_ncols + , '\nGot:', len(aaindex_df.columns)) + + #--------------- + # check for NA + #--------------- + print('\nNow checking for NA in the remaining aaindex_cols') + c1 = aaindex_df.isna().sum() + c2 = c1.sort_values(ascending=False) + print('\nCounting aaindex_df cols with NA' + , '\nncols with NA:', sum(c2>0), 'columns' + , '\nDropping these...' + , '\nOriginal ncols:', len(aaindex_df.columns) + ) + aa_df = aaindex_df.dropna(axis=1) + + print('\nRevised df ncols:', len(aa_df.columns)) + + c3 = aa_df.isna().sum() + c4 = c3.sort_values(ascending=False) + + print('\nChecking NA in revised df...') + + if sum(c4>0): + sys.exit('\nFAIL: aaindex_df still contains cols with NA, please check and drop these before proceeding...') + else: + print('\nPASS: cols with NA successfully dropped from aaindex_df' + , '\nProceeding with combining aa_df with other features_df') + + #--------------------------- + # aaindex: sanity check 2 + #--------------------------- + expected_aa_ncols2 = len(aaindex_df.columns) - sum(c2>0) + if len(aa_df.columns) == expected_aa_ncols2: + print('\nPASS: ncols match' + , '\nExpected ncols:', expected_aa_ncols2 + , '\nGot:', len(aa_df.columns)) + else: + print('\nFAIL: Numbers mismatch' + , '\nExpected ncols:', expected_aa_ncols2 + , '\nGot:', len(aa_df.columns)) + + # Important: need this to identify aaindex cols + aa_df_cols = aa_df.columns + print('\nTotal no. of columns in clean aa_df:', len(aa_df_cols)) + + ############################################################################### + #%% Combining my_features_df and aaindex_df + #=========================== + # Merge my_df + aaindex_df + #=========================== + + if aa_df.columns[aa_df.columns.isin(my_features_df.columns)] == my_features_df.columns[my_features_df.columns.isin(aa_df.columns)]: + print('\nMerging on column: mutationinformation') + + if len(my_features_df) == len(aa_df): + expected_nrows = len(my_features_df) + print('\nProceeding to merge, expected nrows in merged_df:', expected_nrows) + else: + sys.exit('\nNrows mismatch, cannot merge. Please check' + , '\nnrows my_df:', len(my_features_df) + , '\nnrows aa_df:', len(aa_df)) + + #----------------- + # Reset index: mutationinformation + # Very important for merging + #----------------- + aa_df = aa_df.reset_index() + + expected_ncols = len(my_features_df.columns) + len(aa_df.columns) - 1 # for the no. of merging col + + #----------------- + # Merge: my_features_df + aa_df + #----------------- + merged_df = pd.merge(my_features_df + , aa_df + , on = 'mutationinformation') + + #--------------------------- + # aaindex: sanity check 3 + #--------------------------- + if len(merged_df.columns) == expected_ncols: + print('\nPASS: my_features_df and aa_df successfully combined' + , '\nnrows:', len(merged_df) + , '\nncols:', len(merged_df.columns)) + else: + sys.exit('\nFAIL: could not combine my_features_df and aa_df' + , '\nCheck dims and merging cols!') + + #-------- + # Reassign so downstream code doesn't need to change + #-------- + my_df = merged_df.copy() + + #%% Data: my_df + # Check if non structural pos have crept in + # IDEALLY remove from source! But for rpoB do it here + # Drop NA where numerical cols have them + if gene.lower() in geneL_na_ppi2: + #D1148 get rid of + na_index = my_df['mutationinformation'].index[my_df['mcsm_na_affinity'].apply(np.isnan)] + my_df = my_df.drop(index=na_index) + + # FIXED: complete data for all muts inc L114M, F115L, V123L, V125I, V131M + # if gene.lower() in ['embb']: + # na_index = my_df['mutationinformation'].index[my_df['ligand_distance'].apply(np.isnan)] + # my_df = my_df.drop(index=na_index) + + # # Sanity check for non-structural positions + # print('\nChecking for non-structural postions') + # na_index = my_df['mutationinformation'].index[my_df['ligand_distance'].apply(np.isnan)] + # if len(na_index) > 0: + # print('\nNon-structural positions detected for gene:', gene.lower() + # , '\nTotal number of these detected:', len(na_index) + # , '\These are at index:', na_index + # , '\nOriginal nrows:', len(my_df) + # , '\nDropping these...') + # my_df = my_df.drop(index=na_index) + # print('\nRevised nrows:', len(my_df)) + # else: + # print('\nNo non-structural positions detected for gene:', gene.lower() + # , '\nnrows:', len(my_df)) + + + ########################################################################### + #%% Add lineage calculation columns + #FIXME: Check if this can be imported from config? + total_mtblineage_uc = 8 + lineage_colnames = ['lineage_list_all', 'lineage_count_all', 'lineage_count_unique', 'lineage_list_unique', 'lineage_multimode'] + #bar = my_df[lineage_colnames] + my_df['lineage_proportion'] = my_df['lineage_count_unique']/my_df['lineage_count_all'] + my_df['dist_lineage_proportion'] = my_df['lineage_count_unique']/total_mtblineage_uc + ########################################################################### + #%% Active site annotation column + # change from numberic to categorical + + if my_df['active_site'].dtype in num_type: + my_df['active_site'] = my_df['active_site'].astype(object) + my_df['active_site'].dtype + #%% AA property change + #-------------------- + # Water prop change + #-------------------- + my_df['water_change'] = my_df['wt_prop_water'] + str('_to_') + my_df['mut_prop_water'] + my_df['water_change'].value_counts() + + water_prop_changeD = { + 'hydrophobic_to_neutral' : 'change' + , 'hydrophobic_to_hydrophobic' : 'no_change' + , 'neutral_to_neutral' : 'no_change' + , 'neutral_to_hydrophobic' : 'change' + , 'hydrophobic_to_hydrophilic' : 'change' + , 'neutral_to_hydrophilic' : 'change' + , 'hydrophilic_to_neutral' : 'change' + , 'hydrophilic_to_hydrophobic' : 'change' + , 'hydrophilic_to_hydrophilic' : 'no_change' + } + + my_df['water_change'] = my_df['water_change'].map(water_prop_changeD) + my_df['water_change'].value_counts() + + #-------------------- + # Polarity change + #-------------------- + my_df['polarity_change'] = my_df['wt_prop_polarity'] + str('_to_') + my_df['mut_prop_polarity'] + my_df['polarity_change'].value_counts() + + polarity_prop_changeD = { + 'non-polar_to_non-polar' : 'no_change' + , 'non-polar_to_neutral' : 'change' + , 'neutral_to_non-polar' : 'change' + , 'neutral_to_neutral' : 'no_change' + , 'non-polar_to_basic' : 'change' + , 'acidic_to_neutral' : 'change' + , 'basic_to_neutral' : 'change' + , 'non-polar_to_acidic' : 'change' + , 'neutral_to_basic' : 'change' + , 'acidic_to_non-polar' : 'change' + , 'basic_to_non-polar' : 'change' + , 'neutral_to_acidic' : 'change' + , 'acidic_to_acidic' : 'no_change' + , 'basic_to_acidic' : 'change' + , 'basic_to_basic' : 'no_change' + , 'acidic_to_basic' : 'change'} + + my_df['polarity_change'] = my_df['polarity_change'].map(polarity_prop_changeD) + my_df['polarity_change'].value_counts() + + #-------------------- + # Electrostatics change + #-------------------- + my_df['electrostatics_change'] = my_df['wt_calcprop'] + str('_to_') + my_df['mut_calcprop'] + my_df['electrostatics_change'].value_counts() + + calc_prop_changeD = { + 'non-polar_to_non-polar' : 'no_change' + , 'non-polar_to_polar' : 'change' + , 'polar_to_non-polar' : 'change' + , 'non-polar_to_pos' : 'change' + , 'neg_to_non-polar' : 'change' + , 'non-polar_to_neg' : 'change' + , 'pos_to_polar' : 'change' + , 'pos_to_non-polar' : 'change' + , 'polar_to_polar' : 'no_change' + , 'neg_to_neg' : 'no_change' + , 'polar_to_neg' : 'change' + , 'pos_to_neg' : 'change' + , 'pos_to_pos' : 'no_change' + , 'polar_to_pos' : 'change' + , 'neg_to_polar' : 'change' + , 'neg_to_pos' : 'change' + } + + my_df['electrostatics_change'] = my_df['electrostatics_change'].map(calc_prop_changeD) + my_df['electrostatics_change'].value_counts() + + #-------------------- + # Summary change: Create a combined column summarising these three cols + #-------------------- + detect_change = 'change' + check_prop_cols = ['water_change', 'polarity_change', 'electrostatics_change'] + #my_df['aa_prop_change'] = (my_df.values == detect_change).any(1).astype(int) + my_df['aa_prop_change'] = (my_df[check_prop_cols].values == detect_change).any(1).astype(int) + my_df['aa_prop_change'].value_counts() + my_df['aa_prop_change'].dtype + + my_df['aa_prop_change'] = my_df['aa_prop_change'].map({1:'change' + , 0: 'no_change'}) + + my_df['aa_prop_change'].value_counts() + my_df['aa_prop_change'].dtype + + #%% IMPUTE values for OR [check script for exploration: UQ_or_imputer] + #-------------------- + # Impute OR values + #-------------------- + #or_cols = ['or_mychisq', 'log10_or_mychisq', 'or_fisher'] + sel_cols = ['mutationinformation', 'or_mychisq', 'log10_or_mychisq'] + or_cols = ['or_mychisq', 'log10_or_mychisq'] + + print("count of NULL values before imputation\n") + print(my_df[or_cols].isnull().sum()) + + my_dfI = pd.DataFrame(index = my_df['mutationinformation'] ) + + + my_dfI = pd.DataFrame(KNN(n_neighbors=3, weights="uniform").fit_transform(my_df[or_cols]) + , index = my_df['mutationinformation'] + , columns = or_cols ) + my_dfI.columns = ['or_rawI', 'logorI'] + my_dfI.columns + my_dfI = my_dfI.reset_index(drop = False) # prevents old index from being added as a column + my_dfI.head() + print("count of NULL values AFTER imputation\n") + print(my_dfI.isnull().sum()) + + #------------------------------------------- + # OR df Merge: with original based on index + #------------------------------------------- + #my_df['index_bm'] = my_df.index + mydf_imputed = pd.merge(my_df + , my_dfI + , on = 'mutationinformation') + #mydf_imputed = mydf_imputed.set_index(['index_bm']) + + my_df['log10_or_mychisq'].isna().sum() + mydf_imputed['log10_or_mychisq'].isna().sum() + mydf_imputed['logorI'].isna().sum() # should be 0 + + len(my_df.columns) + len(mydf_imputed.columns) + + #----------------------------------------- + # REASSIGN my_df after imputing OR values + #----------------------------------------- + my_df = mydf_imputed.copy() + + if my_df['logorI'].isna().sum() == 0: + print('\nPASS: OR values imputed, data ready for ML') + else: + sys.exit('\nFAIL: something went wrong, Data not ready for ML. Please check upstream!') + + #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + #--------------------------------------- + # TODO: try other imputation like MICE + #--------------------------------------- + #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + #%%######################################################################## + #========================== + # Data for ML + #========================== + my_df_ml = my_df.copy() + + #%% Build X: input for ML + common_cols_stabiltyN = ['ligand_distance' + , 'ligand_affinity_change' + , 'duet_stability_change' + , 'ddg_foldx' + , 'deepddg' + , 'ddg_dynamut2' + , 'mmcsm_lig' + , 'contacts'] + + # Build stability columns ~ gene + if gene.lower() in geneL_basic: + X_stabilityN = common_cols_stabiltyN + cols_to_mask = ['ligand_affinity_change'] + + if gene.lower() in geneL_ppi2: + # X_stabilityN = common_cols_stabiltyN + ['mcsm_ppi2_affinity' , 'interface_dist'] + geneL_ppi2_st_cols = ['mcsm_ppi2_affinity', 'interface_dist'] + X_stabilityN = common_cols_stabiltyN + geneL_ppi2_st_cols + cols_to_mask = ['ligand_affinity_change', 'mcsm_ppi2_affinity'] + + if gene.lower() in geneL_na: + # X_stabilityN = common_cols_stabiltyN + ['mcsm_na_affinity'] + geneL_na_st_cols = ['mcsm_na_affinity'] + X_stabilityN = common_cols_stabiltyN + geneL_na_st_cols + cols_to_mask = ['ligand_affinity_change', 'mcsm_na_affinity'] + + if gene.lower() in geneL_na_ppi2: + # X_stabilityN = common_cols_stabiltyN + ['mcsm_na_affinity'] + ['mcsm_ppi2_affinity', 'interface_dist'] + geneL_na_ppi2_st_cols = ['mcsm_na_affinity'] + ['mcsm_ppi2_affinity', 'interface_dist'] + X_stabilityN = common_cols_stabiltyN + geneL_na_ppi2_st_cols + cols_to_mask = ['ligand_affinity_change', 'mcsm_na_affinity', 'mcsm_ppi2_affinity'] + + + X_foldX_cols = [ 'electro_rr', 'electro_mm', 'electro_sm', 'electro_ss' + , 'disulfide_rr', 'disulfide_mm', 'disulfide_sm', 'disulfide_ss' + , 'hbonds_rr', 'hbonds_mm', 'hbonds_sm', 'hbonds_ss' + , 'partcov_rr', 'partcov_mm', 'partcov_sm', 'partcov_ss' + , 'vdwclashes_rr', 'vdwclashes_mm', 'vdwclashes_sm', 'vdwclashes_ss' + , 'volumetric_rr', 'volumetric_mm', 'volumetric_ss' + ] + + X_str = ['rsa' + #, 'asa' + , 'kd_values' + , 'rd_values'] + + X_ssFN = X_stabilityN + X_str + X_foldX_cols + + X_evolFN = ['consurf_score' + , 'snap2_score' + , 'provean_score'] + + X_genomic_mafor = ['maf' + , 'logorI' + # , 'or_rawI' + # , 'or_mychisq' + # , 'or_logistic' + # , 'or_fisher' + # , 'pval_fisher' + ] + + X_genomic_linegae = ['lineage_proportion' + , 'dist_lineage_proportion' + #, 'lineage' # could be included as a category but it has L2;L4 formatting + , 'lineage_count_all' + , 'lineage_count_unique' + ] + + X_genomicFN = X_genomic_mafor + X_genomic_linegae + + X_aaindexFN = list(aa_df_cols) + + print('\nTotal no. of features for aaindex:', len(X_aaindexFN)) + + # numerical feature names + numerical_FN = X_ssFN + X_evolFN + X_genomicFN + X_aaindexFN + + # categorical feature names + categorical_FN = ['ss_class' + # , 'wt_prop_water' + # , 'mut_prop_water' + # , 'wt_prop_polarity' + # , 'mut_prop_polarity' + # , 'wt_calcprop' + # , 'mut_calcprop' + , 'aa_prop_change' + , 'electrostatics_change' + , 'polarity_change' + , 'water_change' + , 'drtype_mode_labels' # beware then you can't use it to predict [USED it for uq_v1, not v2] + , 'active_site' #[didn't use it for uq_v1] + #, 'gene_name' # will be required for the combined stuff + ] + + #---------------------------------------------- + # count numerical and categorical features + #---------------------------------------------- + + print('\nNo. of numerical features:', len(numerical_FN) + , '\nNo. of categorical features:', len(categorical_FN)) + + #======================= + # Masking columns: + # (mCSM-lig, mCSM-NA, mCSM-ppi2) values for lig_dist >10 + #======================= + # my_df_ml['mutationinformation'][my_df['ligand_distance']>10].value_counts() + # my_df_ml.groupby('mutationinformation')['ligand_distance'].apply(lambda x: (x>10)).value_counts() + + # my_df_ml.loc[(my_df_ml['ligand_distance'] > 10), 'ligand_affinity_change'] = 0 + # (my_df_ml['ligand_affinity_change'] == 0).sum() + + my_df_ml['mutationinformation'][my_df_ml['ligand_distance']>10].value_counts() + my_df_ml.groupby('mutationinformation')['ligand_distance'].apply(lambda x: (x>10)).value_counts() + my_df_ml.loc[(my_df_ml['ligand_distance'] > 10), cols_to_mask].value_counts() + + # mask the mcsm affinity related columns where ligand distance > 10 + my_df_ml.loc[(my_df_ml['ligand_distance'] > 10), cols_to_mask] = 0 + (my_df_ml['ligand_affinity_change'] == 0).sum() + + mask_check = my_df_ml[['mutationinformation', 'ligand_distance'] + cols_to_mask] + + # write file for check + mask_check.sort_values(by = ['ligand_distance'], ascending = True, inplace = True) + mask_check.to_csv(outdir + 'ml/' + gene.lower() + '_mask_check.csv') + + ##################################################################### + #================================================================ + # Training and BLIND test set: 80/20 + + # Throw away previous blind_test_df, and call the 30% data as blind_test + # as these were imputed values and initial analysis shows that this + # is not very representative + #================================================================ + my_df_ml[drug].isna().sum() + # blind_test_df = my_df_ml[my_df_ml[drug].isna()] + # blind_test_df.shape + + training_df = my_df_ml[my_df_ml[drug].notna()] + training_df.shape + + # Target 1: dst_mode + training_df[drug].value_counts() + training_df['dst_mode'].value_counts() + + #################################################################### + +############################################################################### +############################################################################### + # #%% extracting dfs based on numerical, categorical column names + # #---------------------------------- + # # WITHOUT the target var included + # #---------------------------------- + # num_df = training_df[numerical_FN] + # num_df.shape + + # cat_df = training_df[categorical_FN] + # cat_df.shape + + # all_df = training_df[numerical_FN + categorical_FN] + # all_df.shape + + # #------------------------------ + # # WITH the target var included: + # #'wtgt': with target + # #------------------------------ + # # drug and dst_mode should be the same thing + # num_df_wtgt = training_df[numerical_FN + ['dst_mode']] + # num_df_wtgt.shape + + # cat_df_wtgt = training_df[categorical_FN + ['dst_mode']] + # cat_df_wtgt.shape + + # all_df_wtgt = training_df[numerical_FN + categorical_FN + ['dst_mode']] + # all_df_wtgt.shape + + #%%######################################################################## + # #============ + # # ML data: OLD + # #============ + # #------ + # # X: Training and Blind test (BTS) + # #------ + # X = all_df_wtgt[numerical_FN + categorical_FN] # training data ALL + # X_bts = blind_test_df[numerical_FN + categorical_FN] # blind test data ALL + # #X = all_df_wtgt[numerical_FN] # training numerical only + # #X_bts = blind_test_df[numerical_FN] # blind test data numerical + + # #------ + # # y + # #------ + # y = all_df_wtgt['dst_mode'] # training data y + # y_bts = blind_test_df['dst_mode'] # blind data test y + + # #X_bts_wt = blind_test_df[numerical_FN + ['dst_mode']] + +############################################################################### +############################################################################### + #====================================== + # ML data: Train test split: 80/20 + # with stratification + # 80% : training_data for CV + # 20% : blind test + #====================================== + + # features: all_df or + x_features = training_df[numerical_FN + categorical_FN] + y_target = training_df['dst_mode'] + + # sanity check + if not 'dst_mode' in x_features.columns: + print('\nPASS: x_features has no target variable') + x_ncols = len(x_features.columns) + print('\nNo. of columns for x_features:', x_ncols) + # NEED It for scaling law split + #https://towardsdatascience.com/finally-why-we-use-an-80-20-split-for-training-and-test-data-plus-an-alternative-method-oh-yes-edc77e96295d + else: + sys.exit('\nFAIL: x_features has target variable included. FIX it and rerun!') + + #x_train, x_test, y_train, y_test # traditional var_names + # so my downstream code doesn't need to change + X, X_bts, y, y_bts = train_test_split(x_features, y_target + , test_size = 0.2 + , **rs + , stratify = y_target) + yc1 = Counter(y) + yc1_ratio = yc1[0]/yc1[1] + + yc2 = Counter(y_bts) + yc2_ratio = yc2[0]/yc2[1] + + print('\n-------------------------------------------------------------' + , '\nSuccessfully split data with stratification: 80/20 ' + , '\nTrain data size:', X.shape + , '\nTest data size:', X_bts.shape + , '\ny_train numbers:', yc1 + , '\ny_train ratio:',yc1_ratio + , '\n' + , '\ny_test_numbers:', yc2 + , '\ny_test ratio:', yc2_ratio + , '\n-------------------------------------------------------------' + ) + ########################################################################## + # Quick check + #(X['ligand_affinity_change']==0).sum() == (X['ligand_distance']>10).sum() + for i in range(len(cols_to_mask)): + ind = i+1 + print('\nindex:', i, '\nind:', ind) + print('\nMask count check:' + , (my_df_ml[cols_to_mask[i]]==0).sum() == (my_df_ml['ligand_distance']>10).sum() + ) + + print('Original Data\n', Counter(y) + , 'Data dim:', X.shape) + ########################################################################### + #%% + ########################################################################### + # RESAMPLING + ########################################################################### + #------------------------------ + # Simple Random oversampling + # [Numerical + catgeorical] + #------------------------------ + oversample = RandomOverSampler(sampling_strategy='minority') + X_ros, y_ros = oversample.fit_resample(X, y) + print('\nSimple 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('\nSimple 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('\nSimple 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 default + sm_nc = SMOTENC(categorical_features=categorical_colind, k_neighbors = k_sm, **rs, **njobs) + X_smnc, y_smnc = sm_nc.fit_resample(X, y) + print('\nSMOTE_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('\nSMOTE 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('\nSMOTE Over+Under Sampling combined\n', Counter(y_enn)) + + ############################################################################### + # TODO: Find over and undersampling JUST for categorical data diff --git a/scripts/ml/pnca_8020.py b/scripts/ml/pnca_8020.py new file mode 100755 index 0000000..21f45af --- /dev/null +++ b/scripts/ml/pnca_8020.py @@ -0,0 +1,203 @@ +#!/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/') + +from ml_data_8020 import * +setvars(gene,drug) +from ml_data_8020 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 import MultModelsCl + +################################################################################ +print('\n#####################################################################\n' + , '\nRunning ML analysis: 80/20 split' + , '\nGene name:', gene + , '\nDrug name:', drug) + +#================== +# Specify outdir +#================== + +outdir_ml = outdir + 'ml/tts_8020/' + +print('\nOutput directory:', outdir_ml) + +#%%########################################################################### +print('\nSanity checks:' + #, '\nML source data size:', x_features.shape + , '\nTotal input features:', len(X.columns) + , '\n' + , '\nTraining data size:', X.shape + , '\nTest data size:', X_bts.shape + , '\n' + , '\nTarget feature numbers (training data):', Counter(y) + , '\nTarget features ratio (training data:', yc1_ratio + , '\n' + , '\nTarget feature numbers (test data):', Counter(y_bts) + , '\nTarget features ratio (test data):', yc2_ratio + + , '\n\n#####################################################################\n') + +print('\n================================================================\n') + +print('Strucutral features (n):' + , len(X_ssFN) + , '\nThese are:' + , '\nCommon stablity features:', X_stabilityN + , '\nFoldX columns:', X_foldX_cols + , '\nOther struc columns:', X_str + , '\n================================================================\n') + +print('AAindex features (n):' + , len(X_aaindexFN) + , '\nThese are:\n' + , X_aaindexFN + , '\n================================================================\n') + +print('Evolutionary features (n):' + , len(X_evolFN) + , '\nThese are:\n' + , X_evolFN + , '\n================================================================\n') + +print('Genomic features (n):' + , len(X_genomicFN) + , '\nThese are:\n' + , X_genomic_mafor, '\n' + , X_genomic_linegae + , '\n================================================================\n') + +print('Categorical features (n):' + , len(categorical_FN) + , '\nThese are:\n' + , categorical_FN + , '\n================================================================\n') + +if ( len(X.columns) == len(X_ssFN) + len(X_aaindexFN) + len(X_evolFN) + len(X_genomicFN) + len(categorical_FN) ): + print('\nPass: No. of features match') +else: + sys.exit('\nFail: Count of feature mismatch') + +print('\n#####################################################################\n') + +################################################################################ +#================== +# Baseline models +#================== +mm_skf_scoresD = MultModelsCl(input_df = X + , target = y + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) + +baseline_all = pd.DataFrame(mm_skf_scoresD) +baseline_all = baseline_all.T +#baseline_train = baseline_all.filter(like='train_', axis=1) +baseline_CT = baseline_all.filter(like='test_', axis=1) +baseline_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +baseline_BT = baseline_all.filter(like='bts_', axis=1) +baseline_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +baseline_CT.to_csv(outdir_ml + gene.lower() + '_baseline_CT_allF.csv') +baseline_BT.to_csv(outdir_ml + gene.lower() + '_baseline_BT_allF.csv') + + +#%% SMOTE NC: Oversampling [Numerical + categorical] +mm_skf_scoresD7 = MultModelsCl(input_df = X_smnc + , target = y_smnc + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +smnc_all = pd.DataFrame(mm_skf_scoresD7) +smnc_all = smnc_all.T + +smnc_CT = smnc_all.filter(like='test_', axis=1) +smnc_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +smnc_BT = smnc_all.filter(like='bts_', axis=1) +smnc_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +smnc_CT.to_csv(outdir_ml + gene.lower() + '_smnc_CT_allF.csv') +smnc_BT.to_csv(outdir_ml + gene.lower() + '_smnc_BT_allF.csv') + +#%% ROS: Numerical + categorical +mm_skf_scoresD3 = MultModelsCl(input_df = X_ros + , target = y_ros + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +ros_all = pd.DataFrame(mm_skf_scoresD3) +ros_all = ros_all.T + +ros_CT = ros_all.filter(like='test_', axis=1) +ros_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +ros_BT = ros_all.filter(like='bts_', axis=1) +ros_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +ros_CT.to_csv(outdir_ml + gene.lower() + '_ros_CT_allF.csv') +ros_BT.to_csv(outdir_ml + gene.lower() + '_ros_BT_allF.csv') + +#%% RUS: Numerical + categorical +mm_skf_scoresD4 = MultModelsCl(input_df = X_rus + , target = y_rus + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +rus_all = pd.DataFrame(mm_skf_scoresD4) +rus_all = rus_all.T + +rus_CT = rus_all.filter(like='test_', axis=1) +rus_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +rus_BT = rus_all.filter(like='bts_' , axis=1) +rus_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +rus_CT.to_csv(outdir_ml + gene.lower() + '_rus_CT_allF.csv') +rus_BT.to_csv(outdir_ml + gene.lower() + '_rus_BT_allF.csv') + +#%% ROS + RUS Combined: Numerical + categorical +mm_skf_scoresD8 = MultModelsCl(input_df = X_rouC + , target = y_rouC + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +rouC_all = pd.DataFrame(mm_skf_scoresD8) +rouC_all = rouC_all.T + +rouC_CT = rouC_all.filter(like='test_', axis=1) +rouC_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +rouC_BT = rouC_all.filter(like='bts_', axis=1) +rouC_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +rouC_CT.to_csv(outdir_ml + gene.lower() + '_rouC_CT_allF.csv') +rouC_BT.to_csv(outdir_ml + gene.lower() + '_rouC_BT_allF.csv') diff --git a/scripts/ml/rpob_8020.py b/scripts/ml/rpob_8020.py new file mode 100755 index 0000000..021c074 --- /dev/null +++ b/scripts/ml/rpob_8020.py @@ -0,0 +1,203 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat May 28 05:25:30 2022 + +@author: tanu +""" + +import os + +gene = 'rpoB' +drug = 'rifampicin' +#total_mtblineage_uc = 8 + +homedir = os.path.expanduser("~") +os.chdir( homedir + '/git/LSHTM_analysis/scripts/ml/') + +from ml_data_8020 import * +setvars(gene,drug) +from ml_data_8020 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 import MultModelsCl + +################################################################################ +print('\n#####################################################################\n' + , '\nRunning ML analysis: 80/20 split' + , '\nGene name:', gene + , '\nDrug name:', drug) + +#================== +# Specify outdir +#================== + +outdir_ml = outdir + 'ml/tts_8020/' + +print('\nOutput directory:', outdir_ml) + +#%%########################################################################### +print('\nSanity checks:' + #, '\nML source data size:', x_features.shape + , '\nTotal input features:', len(X.columns) + , '\n' + , '\nTraining data size:', X.shape + , '\nTest data size:', X_bts.shape + , '\n' + , '\nTarget feature numbers (training data):', Counter(y) + , '\nTarget features ratio (training data:', yc1_ratio + , '\n' + , '\nTarget feature numbers (test data):', Counter(y_bts) + , '\nTarget features ratio (test data):', yc2_ratio + + , '\n\n#####################################################################\n') + +print('\n================================================================\n') + +print('Strucutral features (n):' + , len(X_ssFN) + , '\nThese are:' + , '\nCommon stablity features:', X_stabilityN + , '\nFoldX columns:', X_foldX_cols + , '\nOther struc columns:', X_str + , '\n================================================================\n') + +print('AAindex features (n):' + , len(X_aaindexFN) + , '\nThese are:\n' + , X_aaindexFN + , '\n================================================================\n') + +print('Evolutionary features (n):' + , len(X_evolFN) + , '\nThese are:\n' + , X_evolFN + , '\n================================================================\n') + +print('Genomic features (n):' + , len(X_genomicFN) + , '\nThese are:\n' + , X_genomic_mafor, '\n' + , X_genomic_linegae + , '\n================================================================\n') + +print('Categorical features (n):' + , len(categorical_FN) + , '\nThese are:\n' + , categorical_FN + , '\n================================================================\n') + +if ( len(X.columns) == len(X_ssFN) + len(X_aaindexFN) + len(X_evolFN) + len(X_genomicFN) + len(categorical_FN) ): + print('\nPass: No. of features match') +else: + sys.exit('\nFail: Count of feature mismatch') + +print('\n#####################################################################\n') + +################################################################################ +#================== +# Baseline models +#================== +mm_skf_scoresD = MultModelsCl(input_df = X + , target = y + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) + +baseline_all = pd.DataFrame(mm_skf_scoresD) +baseline_all = baseline_all.T +#baseline_train = baseline_all.filter(like='train_', axis=1) +baseline_CT = baseline_all.filter(like='test_', axis=1) +baseline_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +baseline_BT = baseline_all.filter(like='bts_', axis=1) +baseline_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +baseline_CT.to_csv(outdir_ml + gene.lower() + '_baseline_CT_allF.csv') +baseline_BT.to_csv(outdir_ml + gene.lower() + '_baseline_BT_allF.csv') + + +#%% SMOTE NC: Oversampling [Numerical + categorical] +mm_skf_scoresD7 = MultModelsCl(input_df = X_smnc + , target = y_smnc + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +smnc_all = pd.DataFrame(mm_skf_scoresD7) +smnc_all = smnc_all.T + +smnc_CT = smnc_all.filter(like='test_', axis=1) +smnc_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +smnc_BT = smnc_all.filter(like='bts_', axis=1) +smnc_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +smnc_CT.to_csv(outdir_ml + gene.lower() + '_smnc_CT_allF.csv') +smnc_BT.to_csv(outdir_ml + gene.lower() + '_smnc_BT_allF.csv') + +#%% ROS: Numerical + categorical +mm_skf_scoresD3 = MultModelsCl(input_df = X_ros + , target = y_ros + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +ros_all = pd.DataFrame(mm_skf_scoresD3) +ros_all = ros_all.T + +ros_CT = ros_all.filter(like='test_', axis=1) +ros_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +ros_BT = ros_all.filter(like='bts_', axis=1) +ros_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +ros_CT.to_csv(outdir_ml + gene.lower() + '_ros_CT_allF.csv') +ros_BT.to_csv(outdir_ml + gene.lower() + '_ros_BT_allF.csv') + +#%% RUS: Numerical + categorical +mm_skf_scoresD4 = MultModelsCl(input_df = X_rus + , target = y_rus + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +rus_all = pd.DataFrame(mm_skf_scoresD4) +rus_all = rus_all.T + +rus_CT = rus_all.filter(like='test_', axis=1) +rus_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +rus_BT = rus_all.filter(like='bts_' , axis=1) +rus_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +rus_CT.to_csv(outdir_ml + gene.lower() + '_rus_CT_allF.csv') +rus_BT.to_csv(outdir_ml + gene.lower() + '_rus_BT_allF.csv') + +#%% ROS + RUS Combined: Numerical + categorical +mm_skf_scoresD8 = MultModelsCl(input_df = X_rouC + , target = y_rouC + , var_type = 'mixed' + , skf_cv = skf_cv + , blind_test_input_df = X_bts + , blind_test_target = y_bts) +rouC_all = pd.DataFrame(mm_skf_scoresD8) +rouC_all = rouC_all.T + +rouC_CT = rouC_all.filter(like='test_', axis=1) +rouC_CT.sort_values(by = ['test_mcc'], ascending = False, inplace = True) + +rouC_BT = rouC_all.filter(like='bts_', axis=1) +rouC_BT.sort_values(by = ['bts_mcc'], ascending = False, inplace = True) + +# Write csv +rouC_CT.to_csv(outdir_ml + gene.lower() + '_rouC_CT_allF.csv') +rouC_BT.to_csv(outdir_ml + gene.lower() + '_rouC_BT_allF.csv')