From 05dd9698c4ac17dc204b18cbd2a1dca9338cf649 Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Fri, 17 Jun 2022 13:41:25 +0100 Subject: [PATCH] added aa_index data for running ml --- scripts/ml/ml_data.py | 275 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 226 insertions(+), 49 deletions(-) mode change 100755 => 100644 scripts/ml/ml_data.py diff --git a/scripts/ml/ml_data.py b/scripts/ml/ml_data.py old mode 100755 new mode 100644 index ec12f20..d5fbe11 --- a/scripts/ml/ml_data.py +++ b/scripts/ml/ml_data.py @@ -29,15 +29,15 @@ def setvars(gene,drug): 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) @@ -50,21 +50,30 @@ def setvars(gene,drug): skf_cv = StratifiedKFold(n_splits = 10 #, shuffle = False, random_state= None) , shuffle = True,**rs) - + rskf_cv = RepeatedStratifiedKFold(n_splits = 10 , n_repeats = 3 , **rs) - + mcc_score_fn = {'mcc': make_scorer(matthews_corrcoef)} jacc_score_fn = {'jcc': make_scorer(jaccard_score)} - + #%% FOR LATER: Combine ED logo data - #%% FOR LATER: active aa site annotations **DONE on 15/05/2022 as part of generating merged_dfs + #%% 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 #============== @@ -75,41 +84,175 @@ def setvars(gene,drug): #======= # input #======= + + #--------- + # File 1 + #--------- infile_ml1 = outdir + gene.lower() + '_merged_df3.csv' #infile_ml2 = outdir + gene.lower() + '_merged_df2.csv' - my_df = pd.read_csv(infile_ml1, index_col = 0) - my_df.dtypes - my_df_cols = my_df.columns + 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 - geneL_basic = ['pnca'] - geneL_na = ['gid'] - geneL_na_ppi2 = ['rpob'] - geneL_ppi2 = ['alr', 'embb', 'katg'] - #%% get cols - mycols = my_df.columns + my_features_df.dtypes + mycols = my_features_df.columns - # # change from numberic to - # num_type = ['int64', 'float64'] - # cat_type = ['object', 'bool'] + #--------- + # File 2 + #--------- + infile_aaindex = outdir + 'aa_index/' + gene.lower() + '_aa.csv' + aaindex_df = pd.read_csv(infile_aaindex, index_col = 0) + aaindex_df.dtypes - # if my_df['active_aa_pos'].dtype in num_type: - # my_df['active_aa_pos'] = my_df['active_aa_pos'].astype(object) - # my_df['active_aa_pos'].dtype + #----------- + # check for non-numerical columns + #----------- + if any(aaindex_df.dtypes==object): + print('\naaindex_df contains non-numerical data') - # FIXME: if this is not structural, remove from source.. + 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) - # FIXME: either impute or remove! - # for embb (L114M, F115L, V123L, V125I, V131M) delete for now - if gene.lower() in ['embb']: - na_index = my_df['mutationinformation'].index[my_df['ligand_distance'].apply(np.isnan)] - my_df = my_df.drop(index=na_index)# RERUN embb with the 5 values now present - + # 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? @@ -119,6 +262,12 @@ def setvars(gene,drug): 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 @@ -137,7 +286,7 @@ def setvars(gene,drug): , '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() @@ -146,7 +295,7 @@ def setvars(gene,drug): #-------------------- 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' @@ -164,7 +313,7 @@ def setvars(gene,drug): , '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() @@ -173,7 +322,7 @@ def setvars(gene,drug): #-------------------- 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' @@ -192,10 +341,10 @@ def setvars(gene,drug): , '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 #-------------------- @@ -208,11 +357,14 @@ def setvars(gene,drug): 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'] @@ -223,7 +375,7 @@ def setvars(gene,drug): my_dfI = pd.DataFrame(index = my_df['mutationinformation'] ) - my_dfI = pd.DataFrame(KNN(n_neighbors= 5, weights="uniform").fit_transform(my_df[or_cols]) + 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'] @@ -236,15 +388,15 @@ def setvars(gene,drug): #------------------------------------------- # OR df Merge: with original based on index #------------------------------------------- - my_df['index_bm'] = my_df.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']) + #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() + mydf_imputed['logorI'].isna().sum() # should be 0 len(my_df.columns) len(mydf_imputed.columns) @@ -253,13 +405,24 @@ def setvars(gene,drug): # 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() - + #========================== # BLIND test set #========================== @@ -309,7 +472,7 @@ def setvars(gene,drug): 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' @@ -347,12 +510,18 @@ def setvars(gene,drug): X_genomicFN = X_genomic_mafor + X_genomic_linegae + X_aaindexFN = list(aa_df_cols) + + print('\nTotal no. of features for aaindex:', len(X_aaindexFN)) + #%% Construct numerical and categorical column names # numerical feature names -# numerical_FN = common_cols_stabiltyN + foldX_cols + X_strFN + X_evolFN + X_genomicFN + # numerical_FN = common_cols_stabiltyN + foldX_cols + X_strFN + X_evolFN + X_genomicFN - numerical_FN = X_ssFN + X_evolFN + X_genomicFN + #numerical_FN = X_ssFN + X_evolFN + X_genomicFN + numerical_FN = X_ssFN + X_evolFN + X_genomicFN + X_aaindexFN + #categorical feature names categorical_FN = ['ss_class' # , 'wt_prop_water' @@ -365,9 +534,17 @@ def setvars(gene,drug): , 'electrostatics_change' , 'polarity_change' , 'water_change' - , 'drtype_mode_labels' # beware then you can use it to predict - #, 'active_aa_pos' # TODO? + #, 'drtype_mode_labels' # beware then you can't use it to predict [USED it for uq_v1] + , 'active_site' + #, '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: @@ -393,7 +570,7 @@ def setvars(gene,drug): # write file for check mask_check.sort_values(by = ['ligand_distance'], ascending = True, inplace = True) mask_check.to_csv(outdir + 'ml/' + gene.lower() + '_mask_check.csv') - + #%% extracting dfs based on numerical, categorical column names #---------------------------------- # WITHOUT the target var included @@ -437,7 +614,7 @@ def setvars(gene,drug): #------ y = all_df_wtgt['dst_mode'] # training data y y_bts = blind_test_df['dst_mode'] # blind data test y - + #X_bts_wt = blind_test_df[numerical_FN + ['dst_mode']] # Quick check