From 478df927ccc1b53a06f3b3c67adfcd2c41015a72 Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Tue, 28 Jun 2022 21:51:02 +0100 Subject: [PATCH] horrible lineage analysis hell --- scripts/count_vars_ML.R | 194 +++-- scripts/data_extraction.py | 66 +- scripts/functions/combining_dfs_plotting.R | 33 +- .../COPY_ml_data_combined_7030.py | 806 ++++++++++++++++++ scripts/ml/combined_model/ml_data_combined | 73 ++ scripts/ml/combined_model/untitled0.py | 221 +++++ scripts/ml/functions/GetMLData.py | 32 +- scripts/ml/ml_data_cd_sl.py | 1 + scripts/ml/test_MultClfs.py | 19 +- scripts/plotting/LINEAGE.R | 325 +++++++ 10 files changed, 1669 insertions(+), 101 deletions(-) create mode 100644 scripts/ml/combined_model/COPY_ml_data_combined_7030.py create mode 100644 scripts/ml/combined_model/ml_data_combined create mode 100644 scripts/ml/combined_model/untitled0.py create mode 100644 scripts/plotting/LINEAGE.R diff --git a/scripts/count_vars_ML.R b/scripts/count_vars_ML.R index 12bc5bd..228d462 100644 --- a/scripts/count_vars_ML.R +++ b/scripts/count_vars_ML.R @@ -2,12 +2,21 @@ #source("~/git/LSHTM_analysis/config/alr.R") #source("~/git/LSHTM_analysis/config/embb.R") -##source("~/git/LSHTM_analysis/config/gid.R") +#source("~/git/LSHTM_analysis/config/gid.R") #source("~/git/LSHTM_analysis/config/katg.R") #source("~/git/LSHTM_analysis/config/pnca.R") source("~/git/LSHTM_analysis/config/rpob.R") -source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") +############################# +# GET the actual merged dfs +############################# +source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") + +############################# +# Output files: merged data +############################# +outfile_merged_df3 = paste0(outdir, '/', tolower(gene), '_merged_df3.csv') +outfile_merged_df2 = paste0(outdir, '/', tolower(gene), '_merged_df2.csv') ################################################ # Add acticve site indication @@ -18,46 +27,90 @@ merged_df2_comp$active_site = as.integer(merged_df2_comp$position %in% active_aa merged_df3$active_site = as.integer(merged_df3$position %in% active_aa_pos) merged_df3_comp$active_site = as.integer(merged_df3_comp$position %in% active_aa_pos) -# sanity check +# check +cols_sel = c('mutationinformation', 'dst', 'mutation_info_labels', 'dm_om_numeric', 'dst_mode') + +check_mdf2 = merged_df2[, cols_sel] +check_mdf2T = table(check_mdf2$mutationinformation, check_mdf2$dst_mode) +ft_mdf2 = as.data.frame.matrix(check_mdf2T) + +#================== +# CHECK: dst mode +#=================== +dst_check = all((ft_mdf2[,1]==0)==(ft_mdf2[,2]!=0)); dst_check + +#======================= +# CHECK: dst mode labels +#======================= +table(merged_df2$mutation_info_labels_orig) +table(merged_df2$mutation_info_labels_v1) +table(merged_df2$mutation_info_labels) + +dst_check1 = table(merged_df2$dst_mode)[1] == table(merged_df2$mutation_info_labels)[2] +dst_check2 = table(merged_df2$dst_mode)[2] == table(merged_df2$mutation_info_labels)[1] + +check12 = all(dst_check && all(dst_check1 == dst_check2)) + +if (check12) { + cat('\nPASS: dst mode labels verified. merged_df3 CAN be trusted! ') +}else{ + stop('FAIL: Something is wrong with the dst_mode column. Quitting!') +``} + +#========================== +# CHECK: active site labels +#========================== table(merged_df2$active_site) table(merged_df3$active_site) - -if( all(table(merged_df2$active_site) == table(as.integer(merged_df2$position %in% active_aa_pos))) && - all(table(merged_df3$active_site) == table(as.integer(merged_df3$position %in% active_aa_pos))) -){ +aa_check1 = all( table(merged_df2$active_site) == table(as.integer(merged_df2$position %in% active_aa_pos)) ) +aa_check2 = all( table(merged_df3$active_site) == table(as.integer(merged_df3$position %in% active_aa_pos)) ) + +if ( all(aa_check1 && aa_check2) ){ cat('\nActive site indications successfully applied to merged_dfs for gene:', tolower(gene)) } - gene gene_match nrow(merged_df3) -############################################## -#============= -# mutation_info: revised labels -#============== -table(merged_df3$mutation_info) -sum(table(merged_df3$mutation_info)) -table(merged_df3$mutation_info_orig) -############################################## +########################################### +#======================== +# CHECK: drtype: revised labels [Merged_df2] +#========================= +table(merged_df2$drtype) #orig +table(merged_df2$drtype_mode) +# mapping 2.1: numeric +# drtype_map = {'XDR': 5 +# , 'Pre-XDR': 4 +# , 'MDR': 3 +# , 'Pre-MDR': 2 +# , 'Other': 1 +# , 'Sensitive': 0} -#============= -# , dst_mode: revised labels -#============== -table(merged_df3$dst) # orig -sum(table(merged_df3$dst)) +# create a labels col that is mapped based on drtype_mode +merged_df2$drtype_mode_labels = merged_df2$drtype_mode +merged_df2$drtype_mode_labels = as.factor(merged_df2$drtype_mode) +levels(merged_df2$drtype_mode_labels) +levels(merged_df2$drtype_mode_labels) <- c('Sensitive', 'Other' + , 'Pre-MDR', 'MDR' + , 'Pre-XDR', 'XDR') +levels(merged_df2$drtype_mode_labels) +# check +a1 = all(table(merged_df2$drtype_mode) == table(merged_df2$drtype_mode_labels)) +b1 = sum(table(merged_df2$drtype_mode_labels)) == nrow(merged_df2) -table(merged_df3$dst_mode) -#table(merged_df3[dr_muts_col]) -sum(table(merged_df3$drtype_mode)) +if (all(a1 && b1)){ + cat("\nPASS: added drtype mode labels to merged_df2") +}else{ + stop("FAIL: could not add drtype mode labels to merged_df2") + ##quit() +} + ################################################# -############################################## -#============= -# drtype: revised labels -#============== +#======================= +# CHECK: drtype: revised labels [merged_df3] +#======================= table(merged_df3$drtype) #orig - table(merged_df3$drtype_mode) # mapping 2.1: numeric # drtype_map = {'XDR': 5 @@ -70,35 +123,80 @@ table(merged_df3$drtype_mode) # create a labels col that is mapped based on drtype_mode merged_df3$drtype_mode_labels = merged_df3$drtype_mode merged_df3$drtype_mode_labels = as.factor(merged_df3$drtype_mode) - levels(merged_df3$drtype_mode_labels) - levels(merged_df3$drtype_mode_labels) <- c('Sensitive', 'Other' , 'Pre-MDR', 'MDR' , 'Pre-XDR', 'XDR') levels(merged_df3$drtype_mode_labels) - +a2 = all(table(merged_df3$drtype_mode) == table(merged_df3$drtype_mode_labels)) +b2 = sum(table(merged_df3$drtype_mode_labels)) == nrow(merged_df3) # check -#table(merged_df3$drtype) -table(merged_df3$drtype_mode) -table(merged_df3$drtype_mode_labels) -sum(table(merged_df3$drtype_mode_labels)) +if (all(a2 && b2)){ + cat("\nPASS: added drtype mode labels to merged_df3") +}else{ + stop("FAIL: could not add drtype mode labels to merged_df3") + ##quit() +} ############################################## -# lineage -table(merged_df3$lineage) -sum(table(merged_df3$lineage_labels)) +#=============== +# CHECK: lineage +#=============== +l1 = table(merged_df3$lineage) == table(merged_df3$lineage_labels) +l2 = table(merged_df2$lineage) == table(merged_df2$lineage_labels) +l3 = sum(table(merged_df2$lineage_labels)) == nrow(merged_df2) +l4 = sum(table(merged_df3$lineage_labels)) == nrow(merged_df3) -cat("\nWriting merged_df3 for:" +if (all(l1 && l2 && l3 && l4) ){ + cat("\nPASS: lineage and lineage labels are identical!") +}else{ + stop("FAIL: could not verify lineage labels") + ##quit() +} + +############################################### +# #============= +# # mutation_info: revised labels +# #============== +# table(merged_df3$mutation_info) +# sum(table(merged_df3$mutation_info)) +# table(merged_df3$mutation_info_orig) +############################################## + +# #============= +# # , dst_mode: revised labels +# #============== +# table(merged_df3$dst) # orig +# sum(table(merged_df3$dst)) +# +# table(merged_df3$dst_mode) +# #table(merged_df3[dr_muts_col]) +# sum(table(merged_df3$drtype_mode)) + +############################################## +if ( all( check12 && aa_check1 && aa_check2 && a1 && b1 && a2 && b2 && l1 && l2 && l3 && l4) ){ + cat("\nWriting merged_dfs for:" , "\nDrug:", drug , "\nGene:", gene) -# write file -outfile_merged_df3 = paste0(outdir, '/', tolower(gene), '_merged_df3.csv') -outfile_merged_df3 -write.csv(merged_df3, outfile_merged_df3) + + write.csv(merged_df3, outfile_merged_df3) + write.csv(merged_df2, outfile_merged_df2) + + cat(paste("\nmerged df3 filename:", outfile_merged_df3 + , "\nmerged df2 filename:", outfile_merged_df2)) + +} else{ + stop("FAIL: Not able to write merged dfs. Please check numbers!") + #quit() +} -outfile_merged_df2 = paste0(outdir, '/', tolower(gene), '_merged_df2.csv') -outfile_merged_df2 -write.csv(merged_df2, outfile_merged_df2) +# write file +# outfile_merged_df3 = paste0(outdir, '/', tolower(gene), '_merged_df3.csv') +# outfile_merged_df3 +# write.csv(merged_df3, outfile_merged_df3) +# +# outfile_merged_df2 = paste0(outdir, '/', tolower(gene), '_merged_df2.csv') +# outfile_merged_df2 +# write.csv(merged_df2, outfile_merged_df2) ################################################### ################################################### @@ -133,5 +231,3 @@ write.csv(merged_df2, outfile_merged_df2) # # drtype: MDR and XDR # #table(df3$drtype) orig i.e. incorrect ones! # table(df3$drtype_mode_labels) -# -# diff --git a/scripts/data_extraction.py b/scripts/data_extraction.py index 20de8ca..1e118ba 100755 --- a/scripts/data_extraction.py +++ b/scripts/data_extraction.py @@ -60,9 +60,12 @@ import collections #%% dir and local imports homedir = os.path.expanduser('~') # set working dir -os.getcwd() -os.chdir(homedir + '/git/LSHTM_analysis/scripts') -os.getcwd() +# os.getcwd() +# os.chdir(homedir + '/git/LSHTM_analysis/scripts') +# os.getcwd() + +sys.path.append(homedir + '/git/LSHTM_analysis/scripts/') + #======================================================================= #%% command line args arg_parser = argparse.ArgumentParser() @@ -1550,6 +1553,29 @@ gene_LF3['dst_multimode'].value_counts() #gene_LF3['dst_noNA'] = gene_LF3['dst_multimode'].apply(lambda x: np.nanmax(x)) gene_LF3['dst_mode'] = gene_LF3['dst_multimode'].apply(lambda x: np.nanmax(x)) #ML +#----------------------------------------------------------------------------- +#----------------------------------------------------------------------------- +# NOTE: unexpected weirdness with above, so redoing it! +mmdf = pd.DataFrame(gene_LF3.groupby('mutationinformation')['dst_mode'].agg(multimode)) +mmdf['dst2'] = mmdf['dst_mode'].apply(lambda x: int(max(x))) +mmdf=mmdf.reset_index() + +# rename cols to make sure merge will have the names you expect +mmdf2 = mmdf.rename(columns = {'dst_mode':'dst_multimode', 'dst2':'dst_mode'}) + +# IMPORTANT! +gene_LF3_copy = gene_LF3.copy() +gene_LF3_copy.drop(["dst_mode", "dst_multimode", "dst_multimode_all"], axis = 1, inplace = True) + +# Now merge gene_LF3.copy and mmdf2 +gene_LF3_merged = pd.merge(gene_LF3_copy, mmdf2, on='mutationinformation') +df_check4 = gene_LF3_merged[['mutationinformation', 'dst', 'dst_multimode', 'dst_mode', 'position' ]] + +# now reassign the merged df to gene_LF3 for integration with downstream +gene_LF3 = gene_LF3_merged.copy() +#----------------------------------------------------------------------------- +#----------------------------------------------------------------------------- + # sanity checks #gene_LF3['dst_noNA'].equals(gene_LF3['dst_mode']) gene_LF3[drug].value_counts() @@ -1700,10 +1726,24 @@ lf_lin_split['lineage_numeric'].value_counts() # Add lineage_list: ALL values: #-------------------------------- # Add all lineages for each mutation -lf_lin_split['lineage_corrupt_list'] = lf_lin_split['lineage_corrupt'] -lf_lin_split['lineage_corrupt_list'].value_counts() +#lf_lin_split['lineage_corrupt_list'] = lf_lin_split['lineage_corrupt'].copy() +#lf_lin_split['lineage_corrupt_list'].value_counts() +lf_lin_split['lineage_corrupt'].value_counts() + #lf_lin_split['lineage_corrupt_list'] = lf_lin_split['mutationinformation'].map(lf_lin_split.groupby('mutationinformati -lf_lin_split['lineage_corrupt_list'] = lf_lin_split.groupby('mutationinformation').lineage_corrupt_list.apply(list) +#lf_lin_split['lineage_corrupt_list'] = lf_lin_split.groupby('mutationinformation').lineage_corrupt_list.apply(list) +lf_lin_tmp =lf_lin_split.groupby('Mut').lineage_corrupt.apply(list) +lf_lin_tmp = lf_lin_tmp.reset_index() +lf_lin_tmp.rename(columns={'lineage_corrupt': 'lineage_corrupt_list' }, inplace=True) +#lf_lin_split['lineage_corrupt_list'] = lf_lin_split.groupby('Mut').lineage_corrupt_list.apply(list).copy() + +#lf_lin_split['lineage_corrupt_list'] = lf_lin_tmp +lf_lin_merged = pd.merge(lf_lin_split, lf_lin_tmp, on='Mut') +lf_lin_split.shape +lf_lin_merged.shape + +# REASSIGN merged +lf_lin_split = lf_lin_merged.copy() lf_lin_split['lineage_corrupt_list'].value_counts() #-------------------------------- @@ -1727,10 +1767,18 @@ lf_lin_split['lineage_ulist'] = lf_lin_split['lineage_set'].apply(lambda x : li #------------------------------------- # Lineage numeric mode: multimode #------------------------------------- -lf_lin_split['lineage_multimode'] = lf_lin_split.groupby('mutationinformation')['lineage_numeric'].agg(multimode) -lf_lin_split['lineage_multimode'].value_counts() +#lf_lin_split['lineage_multimode'] = lf_lin_split.groupby('mutationinformation')['lineage_numeric'].agg(multimode) +#lf_lin_split['lineage_multimode'] = lf_lin_split.groupby('Mut')['lineage_numeric'].agg(multimode) -# cant take max as it doesn't mean anyting! +lin_mm_tmp = pd.DataFrame(lf_lin_split.groupby('Mut')['lineage_numeric'].agg(multimode)) +lin_mm_tmp=lin_mm_tmp.reset_index() +lin_mm_tmp.rename(columns={'lineage_numeric':'lineage_multimode'}, inplace=True) + +lf_lin_split_merged = pd.merge(lf_lin_split, lin_mm_tmp, on='Mut') +#lf_lin_split['lineage_multimode'].value_counts() # cant take max as it doesn't mean anyting! + +#REASSIGN +lf_lin_split = lf_lin_split_merged.copy() ############################################################################### #%% Select only the columns you want to merge from lf_lin_split diff --git a/scripts/functions/combining_dfs_plotting.R b/scripts/functions/combining_dfs_plotting.R index 88f62b6..26b214f 100644 --- a/scripts/functions/combining_dfs_plotting.R +++ b/scripts/functions/combining_dfs_plotting.R @@ -185,9 +185,6 @@ combining_dfs_plotting <- function( my_df_u } - - - # Quick formatting: ordering df and pretty labels #------------------------------ @@ -198,36 +195,12 @@ combining_dfs_plotting <- function( my_df_u #----------------------- # mutation_info_labels #----------------------- - merged_df2$mutation_info_labels = ifelse(merged_df2$mutation_info == dr_muts_col - , "DM", "OM") - merged_df2$mutation_info_labels = factor(merged_df2$mutation_info_labels) + #merged_df2$mutation_info_labels = ifelse(merged_df2$mutation_info == dr_muts_col + # , "DM", "OM") + #merged_df2$mutation_info_labels = factor(merged_df2$mutation_info_labels) #----------------------- # lineage labels #----------------------- - # merged_df2$lineage_labels = gsub("lineage", "L", merged_df2$lineage) - - # Already solved upstream where lineage now contains L prefix - # merged_df2$lineage_labels = factor(merged_df2$lineage_labels, c("L1" - # , "L2" - # , "L3" - # , "L4" - # , "L5" - # , "L6" - # , "L7" - # , "LBOV" - # , "L1;L2" - # , "L1;L3" - # , "L1;L4" - # , "L2;L3" - # , "L2;L3;L4" - # , "L2;L4" - # , "L2;L6" - # , "L2;LBOV" - # , "L3;L4" - # , "L4;L6" - # , "L4;L7" - # , "")) - merged_df2$lineage_labels = merged_df2$lineage #merged_df2$lineage_labels = as.factor(merged_df2$lineage_labels) #merged_df2$lineage_labels = factor(merged_df2$lineage_labels) diff --git a/scripts/ml/combined_model/COPY_ml_data_combined_7030.py b/scripts/ml/combined_model/COPY_ml_data_combined_7030.py new file mode 100644 index 0000000..274aff6 --- /dev/null +++ b/scripts/ml/combined_model/COPY_ml_data_combined_7030.py @@ -0,0 +1,806 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sun Mar 6 13:41:54 2022 + +@author: tanu +""" +def setvars(gene,drug): + #https://stackoverflow.com/questions/51695322/compare-multiple-algorithms-with-sklearn-pipeline + import os, sys + import pandas as pd + import numpy as np + print(np.__version__) + print(pd.__version__) + import pprint as pp + from copy import deepcopy + from collections import Counter + from sklearn.impute import KNNImputer as KNN + from imblearn.over_sampling import RandomOverSampler + from imblearn.under_sampling import RandomUnderSampler + from imblearn.over_sampling import SMOTE + from sklearn.datasets import make_classification + from imblearn.combine import SMOTEENN + from imblearn.combine import SMOTETomek + + from imblearn.over_sampling import SMOTENC + from imblearn.under_sampling import EditedNearestNeighbours + from imblearn.under_sampling import RepeatedEditedNearestNeighbours + + from sklearn.metrics import make_scorer, confusion_matrix, accuracy_score, balanced_accuracy_score, precision_score, average_precision_score, recall_score + from sklearn.metrics import roc_auc_score, roc_curve, f1_score, matthews_corrcoef, jaccard_score, classification_report + + from sklearn.model_selection import train_test_split, cross_validate, cross_val_score + from sklearn.model_selection import StratifiedKFold,RepeatedStratifiedKFold, RepeatedKFold + + from sklearn.pipeline import Pipeline, make_pipeline + import argparse + import re + #%% GLOBALS + tts_split = "70_30" + + 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 + ########################################################################### + + 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 column names to mask for affinity chanhes + if gene.lower() in geneL_basic: + #X_stabilityN = common_cols_stabiltyN + gene_affinity_colnames = []# not needed as its the common ones + cols_to_mask = ['ligand_affinity_change'] + + if gene.lower() in geneL_ppi2: + gene_affinity_colnames = ['mcsm_ppi2_affinity', 'interface_dist'] + #X_stabilityN = common_cols_stabiltyN + geneL_ppi2_st_cols + cols_to_mask = ['ligand_affinity_change', 'mcsm_ppi2_affinity'] + + if gene.lower() in geneL_na: + gene_affinity_colnames = ['mcsm_na_affinity'] + #X_stabilityN = common_cols_stabiltyN + geneL_na_st_cols + cols_to_mask = ['ligand_affinity_change', 'mcsm_na_affinity'] + + if gene.lower() in geneL_na_ppi2: + gene_affinity_colnames = ['mcsm_na_affinity'] + ['mcsm_ppi2_affinity', 'interface_dist'] + #X_stabilityN = common_cols_stabiltyN + geneL_na_ppi2_st_cols + cols_to_mask = ['ligand_affinity_change', 'mcsm_na_affinity', 'mcsm_ppi2_affinity'] + + #======================= + # Masking columns: + # (mCSM-lig, mCSM-NA, mCSM-ppi2) values for lig_dist >10 + #======================= + my_df_ml['mutationinformation'][my_df_ml['ligand_distance']>10].value_counts() + my_df_ml.groupby('mutationinformation')['ligand_distance'].apply(lambda x: (x>10)).value_counts() + my_df_ml.loc[(my_df_ml['ligand_distance'] > 10), cols_to_mask].value_counts() + + # mask the mcsm affinity related columns where ligand distance > 10 + my_df_ml.loc[(my_df_ml['ligand_distance'] > 10), cols_to_mask] = 0 + (my_df_ml['ligand_affinity_change'] == 0).sum() + + mask_check = my_df_ml[['mutationinformation', 'ligand_distance'] + cols_to_mask] + + #=================================================== + # write file for check + mask_check.sort_values(by = ['ligand_distance'], ascending = True, inplace = True) + mask_check.to_csv(outdir + 'ml/' + gene.lower() + '_mask_check.csv') + #=================================================== + ############################################################################### + #%% Feature groups (FG): Build X for Input ML + ############################################################################ + #=========================== + # FG1: Evolutionary features + #=========================== + X_evolFN = ['consurf_score' + , 'snap2_score' + , 'provean_score'] + + ############################################################################### + #======================== + # FG2: Stability features + #======================== + #-------- + # common + #-------- + X_common_stability_Fnum = [ + 'duet_stability_change' + , 'ddg_foldx' + , 'deepddg' + , 'ddg_dynamut2' + , 'contacts'] + #-------- + # FoldX + #-------- + X_foldX_Fnum = [ 'electro_rr', 'electro_mm', 'electro_sm', 'electro_ss' + , 'disulfide_rr', 'disulfide_mm', 'disulfide_sm', 'disulfide_ss' + , 'hbonds_rr', 'hbonds_mm', 'hbonds_sm', 'hbonds_ss' + , 'partcov_rr', 'partcov_mm', 'partcov_sm', 'partcov_ss' + , 'vdwclashes_rr', 'vdwclashes_mm', 'vdwclashes_sm', 'vdwclashes_ss' + , 'volumetric_rr', 'volumetric_mm', 'volumetric_ss'] + + X_stability_FN = X_common_stability_Fnum + X_foldX_Fnum + + ############################################################################### + #=================== + # FG3: Affinity features + #=================== + common_affinity_Fnum = ['ligand_distance' + , 'ligand_affinity_change' + , 'mmcsm_lig'] + + # if gene.lower() in geneL_basic: + # X_affinityFN = common_affinity_Fnum + # else: + # X_affinityFN = common_affinity_Fnum + gene_affinity_colnames + + X_affinityFN = common_affinity_Fnum + gene_affinity_colnames + + ############################################################################### + #============================ + # FG4: Residue level features + #============================ + #----------- + # AA index + #----------- + X_aaindex_Fnum = list(aa_df_cols) + print('\nTotal no. of features for aaindex:', len(X_aaindex_Fnum)) + + #----------------- + # surface area + # depth + # hydrophobicity + #----------------- + X_str_Fnum = ['rsa' + #, 'asa' + , 'kd_values' + , 'rd_values'] + + #--------------------------- + # Other aa properties + # active site indication + #--------------------------- + X_aap_Fcat = ['ss_class' + # , 'wt_prop_water' + # , 'mut_prop_water' + # , 'wt_prop_polarity' + # , 'mut_prop_polarity' + # , 'wt_calcprop' + # , 'mut_calcprop' + , 'aa_prop_change' + , 'electrostatics_change' + , 'polarity_change' + , 'water_change' + , 'active_site'] + + X_resprop_FN = X_aaindex_Fnum + X_str_Fnum + X_aap_Fcat + ############################################################################### + #======================== + # FG5: Genomic features + #======================== + X_gn_mafor_Fnum = ['maf' + #, 'logorI' + # , 'or_rawI' + # , 'or_mychisq' + # , 'or_logistic' + # , 'or_fisher' + # , 'pval_fisher' + ] + + X_gn_linegae_Fnum = ['lineage_proportion' + , 'dist_lineage_proportion' + #, 'lineage' # could be included as a category but it has L2;L4 formatting + , 'lineage_count_all' + , 'lineage_count_unique' + ] + + # X_gn_Fcat = ['drtype_mode_labels' # beware then you can't use it to predict [USED it for uq_v1, not v2] + # #, 'gene_name' # will be required for the combined stuff + # ] + X_gn_Fcat = [] + + 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: 70/30 + # dst with actual values : training set + # dst with imputed values : THROW AWAY [unrepresentative] + #================================================================ + 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() + + #################################################################### + #==================================== + # ML data: Train test split: 70/30 + # with stratification + # 70% : training_data for CV + # 30% : blind test + #===================================== + x_features = training_df[all_featuresN] + 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!') + #------------------- + # train-test split + #------------------- + #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.33 + , **rs + , stratify = y_target) + 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' + , '\nSplit:', tts_split + #, '\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-------------------------------------------------------------' + ) + ########################################################################## + # 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 + ########################################################################### + + print('\n#################################################################' + , '\nDim of X for gene:', gene.lower(), '\n', X.shape + , '\n###############################################################') diff --git a/scripts/ml/combined_model/ml_data_combined b/scripts/ml/combined_model/ml_data_combined new file mode 100644 index 0000000..f2e8198 --- /dev/null +++ b/scripts/ml/combined_model/ml_data_combined @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat Jun 25 11:07:30 2022 + +@author: tanu +""" + +import sys, os +import pandas as pd +import numpy as np +import re + +############################################################################### +homedir = os.path.expanduser("~") +sys.path.append(homedir + '/git/LSHTM_analysis/scripts/ml/functions') +############################################################################### +#==================== +# Import ML functions +#==================== +#from MultClfs import * +from GetMLData import * +from SplitTTS import * +#%% Load all gene files ####################################################### +# param dict +combined_model_paramD = {'data_combined_model' : True + , 'use_or' : False + , 'omit_all_genomic_features': False + , 'write_maskfile' : False + , 'write_outfile' : False } + +pnca_df = getmldata('pncA', 'pyrazinamide' , **combined_model_paramD) +embb_df = getmldata('embB', 'ethambutol' , **combined_model_paramD) +katg_df = getmldata('katG', 'isoniazid' , **combined_model_paramD) +rpob_df = getmldata('rpoB', 'rifampicin' , **combined_model_paramD) +gid_df = getmldata('gid' , 'streptomycin' , **combined_model_paramD) +alr_df = getmldata('alr' , 'cycloserine' , **combined_model_paramD) + +# quick check +foo = pd.concat([alr_df, pnca_df]) +check1 = foo.filter(regex= '.*_affinity|gene_name|ligand_distance', axis = 1) +# So, pd.concat will join correctly but introduce NAs. +# TODO: discuss whether to make these 0 and use it or just omit +# For now I am omitting these i.e combining only on common columns + +expected_nrows = len(pnca_df) + len(embb_df) + len(katg_df) + len(rpob_df) + len(gid_df) + len(alr_df) + +# finding common columns +dfs_combine = [pnca_df, embb_df, katg_df, rpob_df, gid_df, alr_df] +common_cols = list(set.intersection(*(set(df.columns) for df in dfs_combine))) +expected_ncols = np.min([len(pnca_df.columns)] + [len(embb_df.columns)] + [len(katg_df.columns)] + [len(rpob_df.columns)] + [len(gid_df.columns)] + [len(alr_df.columns)]) +expected_ncols + +if len(common_cols) == expected_ncols: + print('\nProceeding to combine based on common cols (n):', len(common_cols)) + combined_df = pd.concat([df[common_cols] for df in dfs_combine], ignore_index = False) + print('\nSuccessfully combined dfs:' + , '\nNo. of dfs combined:', len(dfs_combine) + , '\nDim of combined df:', combined_df.shape) +else: + print('\nFAIL: could not combine dfs, length mismatch' + , '\nExpected ncols:', expected_ncols + , '\nGot:', len(common_cols)) +#%% split data into different data types +tts_7030_paramD = {'data_type' : 'actual' + , 'split_type' : '70_30' + , 'oversampling' : True} + +data_CM_7030D = split_tts(ml_input_data = combined_df + , **tts_7030_paramD + , dst_colname = 'dst' + , target_colname = 'dst_mode' + , include_gene_name = False) # when not doing leave one group out \ No newline at end of file diff --git a/scripts/ml/combined_model/untitled0.py b/scripts/ml/combined_model/untitled0.py new file mode 100644 index 0000000..00a8470 --- /dev/null +++ b/scripts/ml/combined_model/untitled0.py @@ -0,0 +1,221 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat Jun 25 11:07:30 2022 + +@author: tanu +""" + +import sys, os +import pandas as pd +import numpy as np +import os, sys +import pandas as pd +import numpy as np +print(np.__version__) +print(pd.__version__) +import pprint as pp +from copy import deepcopy +from collections import Counter +from sklearn.impute import KNNImputer as KNN +from imblearn.over_sampling import RandomOverSampler +from imblearn.under_sampling import RandomUnderSampler +from imblearn.over_sampling import SMOTE +from sklearn.datasets import make_classification +from imblearn.combine import SMOTEENN +from imblearn.combine import SMOTETomek + +from imblearn.over_sampling import SMOTENC +from imblearn.under_sampling import EditedNearestNeighbours +from imblearn.under_sampling import RepeatedEditedNearestNeighbours + +from sklearn.metrics import make_scorer, confusion_matrix, accuracy_score, balanced_accuracy_score, precision_score, average_precision_score, recall_score +from sklearn.metrics import roc_auc_score, roc_curve, f1_score, matthews_corrcoef, jaccard_score, classification_report + +from sklearn.model_selection import train_test_split, cross_validate, cross_val_score +from sklearn.model_selection import StratifiedKFold,RepeatedStratifiedKFold, RepeatedKFold + +from sklearn.pipeline import Pipeline, make_pipeline +import argparse +import re +homedir = os.path.expanduser("~") +#%% Globals +rs = {'random_state': 42} +njobs = {'n_jobs': 10} +#%% Define split_tts function ################################################# +def split_tts(ml_input_data + , data_type = ['actual', 'complete', 'reverse'] + , split_type = ['70_30', '80_20', 'sl'] + , oversampling = True + , dst_colname = 'dst'# determine how to subset the actual vs reverse data + , target_colname = 'dst_mode'): + + print('\nInput params:' + , '\nDim of input df:' , ml_input_data.shape + , '\nData type to split:', data_type + , '\nSplit type:' , split_type + , '\ntarget colname:' , target_colname) + + if oversampling: + print('\noversampling enabled') + else: + print('\nNot generating oversampled or undersampled data') + + #==================================== + # evaluating use_data_type + #==================================== + if data_type == 'actual': + ml_data = ml_input_data[ml_input_data[dst_colname].notna()] + if data_type == 'complete': + ml_data = ml_input_data.copy() + if data_type == 'reverse': + ml_data = ml_input_data[ml_input_data[dst_colname].isna()] + #if_data_type == none + + #==================================== + # separate features and target + #==================================== + x_features = ml_data.drop([target_colname, dst_colname], axis = 1) + y_target = ml_data[target_colname] + + # 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) + else: + sys.exit('\nFAIL: x_features has target variable included. FIX it and rerun!') + + #==================================== + # Train test split + # with stratification + #===================================== + if split_type == '70_30': + tts_test_size = 0.33 + if split_type == '80_20': + tts_test_size = 0.2 + if split_type == 'sl': + tts_test_size = 1/np.sqrt(x_ncols) + train_sl = 1 - tts_test_size + + #------------------------- + # TTS split ~ split_type + #------------------------- + #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 = tts_test_size + , **rs + , stratify = y_target) + 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 + ############################################################################### + print('\n-------------------------------------------------------------' + , '\nSuccessfully generated training and test data:' + , '\nData used:' , data_type + , '\nSplit type:', split_type + + , '\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 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-------------------------------------------------------------' + ) + + if oversampling: + + ####################################################################### + # 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 # 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) + + print('\nGenerated resampled data as below:' + , '\n===========================' + , '\nRandom oversampling:' + , '\n===========================' + + , '\n\nTrain data size:', X_ros.shape + + , '\ny_train numbers:', y_ros + , '\n\ny_train ratio:', Counter(y_ros)[0]/Counter(y_ros)[0] + + , '\ny_test ratio:' , yc2_ratio + + , '\n-------------------------------------------------------------' + ) + + + # globals().update(locals()) # TROLOLOLOLOLOLS + + #return() \ No newline at end of file diff --git a/scripts/ml/functions/GetMLData.py b/scripts/ml/functions/GetMLData.py index f8155d9..37ecc1c 100755 --- a/scripts/ml/functions/GetMLData.py +++ b/scripts/ml/functions/GetMLData.py @@ -603,19 +603,20 @@ def getmldata(gene, drug # training_df[drug].value_counts() # training_df['dst_mode'].value_counts() - all_training_df = my_df_ml[all_featuresN] + #all_training_df = my_df_ml[all_featuresN] + # Getting the dst column as this will be required for tts_split() + if 'dst' in my_df_ml: + print('\ndst column exists') + if my_df_ml['dst'].equals(my_df_ml[drug]): + print('\nand this is identical to drug column:', drug) + + all_featuresN2 = all_featuresN + ['dst', 'dst_mode'] + all_training_df = my_df_ml[all_featuresN2] + + print('\nAll feature names:', all_featuresN2) #################################################################### - print('\n#################################################################' - , '\nSUCCESS: Extacted training data for gene:', gene.lower() - , '\nDim of training_df:', all_training_df.shape) - if use_or: - print('\nThis includes Odds Ratio') - else: - print('\nThis EXCLUDES Odds Ratio' - , '\n###############################################################') - #========================================================================== if write_maskfile: print('\nPASS: and now writing file to check masked columns and values:', outFile_mask_ml ) @@ -630,4 +631,15 @@ def getmldata(gene, drug else: print('\nPASS: But NOT writing processed file') #========================================================================== + + print('\n#################################################################' + , '\nSUCCESS: Extacted training data for gene:', gene.lower() + , '\nDim of training_df:', all_training_df.shape) + if use_or: + print('\nThis includes Odds Ratio' + , '\n###########################################################') + else: + print('\nThis EXCLUDES Odds Ratio' + , '\n############################################################') + return(all_training_df) \ No newline at end of file diff --git a/scripts/ml/ml_data_cd_sl.py b/scripts/ml/ml_data_cd_sl.py index 5941c30..2c08ef1 100755 --- a/scripts/ml/ml_data_cd_sl.py +++ b/scripts/ml/ml_data_cd_sl.py @@ -753,6 +753,7 @@ def setvars(gene,drug): 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) diff --git a/scripts/ml/test_MultClfs.py b/scripts/ml/test_MultClfs.py index f337470..f053c62 100755 --- a/scripts/ml/test_MultClfs.py +++ b/scripts/ml/test_MultClfs.py @@ -67,7 +67,7 @@ print('\n#####################################################################\n #================ # MultModelsCl: without formatted output #================ -mmD = MultModelsCl(input_df = X_smnc +mmD = MultModelsCl_noBT(input_df = X_smnc , target = y_smnc , var_type = 'mixed' , tts_split_type = tts_split_7030 @@ -77,12 +77,13 @@ mmD = MultModelsCl(input_df = X_smnc , blind_test_target = y_bts , add_cm = True , add_yn = True + , run_blind_test = True , return_formatted_output = False) #================ # MultModelsCl: WITH formatted output #================ -mmDF3 = MultModelsCl(input_df = X_smnc +mmDF3 = MultModelsCl_noBT(input_df = X_smnc , target = y_smnc , var_type = 'mixed' , tts_split_type = tts_split_7030 @@ -92,9 +93,21 @@ mmDF3 = MultModelsCl(input_df = X_smnc , blind_test_target = y_bts , add_cm = True , add_yn = True + , run_blind_test = True , return_formatted_output= True ) - +mmDF9= MultModelsCl_noBT(input_df = X + , target = y + , var_type = 'mixed' + , tts_split_type = tts_split_7030 + , resampling_type = 'none' + , skf_cv = None + , blind_test_df = X_bts + , blind_test_target = y_bts + , add_cm = True + , add_yn = True + , run_blind_test = True + , return_formatted_output= True ) #================= # test function #================= diff --git a/scripts/plotting/LINEAGE.R b/scripts/plotting/LINEAGE.R new file mode 100644 index 0000000..f5b1795 --- /dev/null +++ b/scripts/plotting/LINEAGE.R @@ -0,0 +1,325 @@ +library(tidyverse) +#install.packages("ggforce") +library("ggforce") +#install.packages("gginference") +library(gginference) +library(ggpubr) + +#%% read data +df = read.csv("/home/tanu/git/Data/pyrazinamide/output/pnca_merged_df2.csv") +#df = read.csv("/home/tanu/git/Data/pyrazinamide/output/pnca_merged_df3.csv") + +foo = as.data.frame(colnames(df)) + +my_df = df[ ,c('mutationinformation' + , 'snp_frequency' + , 'pos_count' + , 'lineage' + , 'lineage_multimode' + , 'dst' + , 'dst_mode')] + +#%% create sensitivity column ~ dst_mode +my_df$sensitivity = ifelse(my_df$dst_mode == 1, "R", "S") +table(my_df$dst_mode) +table(my_df$sensitivity) + +test = my_df[my_df$mutationinformation=="A102P",] + + + + +# fix the lineage_multimode labels +my_df$lineage_multimode +my_df$lineage_mm <- gsub("\\.0", "", my_df$lineage_multimode) +my_df$lineage_mm + +my_df$lineage_mm <- gsub("\\[|||]", "", my_df$lineage_mm) +str(my_df$lineage_mm) +table(my_df$lineage_mm) + +my_dfF = separate_rows(my_df, lineage_mm, sep = ",") +my_dfF = as.data.frame(my_dfF) + +table(my_dfF$lineage_mm) +my_dfF$lineage_mm <- gsub(" ", "", my_dfF$lineage_mm) +table(my_dfF$lineage_mm) + +# addd prefix L +my_dfF$lineage_mm = paste0("L", my_dfF$lineage_mm) +table(my_dfF$lineage_mm) + +if (class(my_df) == class(my_dfF)){ + cat('\nPASS: separated lineage multimode label column') + my_df = my_dfF +} else{ + cat('\nFAIL: could not split lineage multimode column') +} + +# select only L1-L4 and LBOV +sel_lineages = c("L1", "L2", "L3", "L4") +table(my_df$lineage_mm) +my_df2 = my_df[my_df$lineage_mm%in%sel_lineages,] +table(my_df2$lineage) +sum(table(my_df2$lineage_mm)) == nrow(my_df2) + + +dup_rows = my_df2[duplicated(my_df2[c('mutationinformation')]), ] +expected_nrows = nrow(my_df2) - nrow(dup_rows) +my_df3 = my_df2[!duplicated(my_df2[c('mutationinformation')]), ] + +if ( nrow(my_df3) == expected_nrows ) { + cat('\nPASS: duplicated rows removed') +}else{ + cat('\nFAIL: duplicated rows could not be removed') +} + +table(my_df3$lineage_mm) +str(my_df3$lineage_mm) + +# convert to factor +str(my_df3) +my_df3$lineage = as.factor(my_df3$lineage) +my_df3$lineage_mm = as.factor(my_df3$lineage_mm) +my_df3$sensitivity = as.factor(my_df3$sensitivity) + +str(my_df3$lineage_mm) + +#df2 = my_df2[1:100,] +df2 = my_df3 +sum(table(df2$mutationinformation)) + +table(df2$lineage_mm) +str(df2$lineage_mm) + +#df3 = df2[na.omit(df2$dst)] +#sum(is.na(df2$dst)) +df3 = df2[!is.na(df2$dst), ] +nrow(df3) + +#%% plot +#============ +# facet wrap +#============ +plot_data = df2 +plot_data = df3 +table(plot_data$mutationinformation, plot_data$lineage_mm, plot_data$dst) + +test2 = my_df[1:500, ] +test2 = my_df +test2 = test2[test2$lineage%in%sel_lineages,] +nrow(test2) + +# stats +f2 = test2[test2$mutationinformation == "Y95D",] +h = table(f2$lineage, f2$dst); h +h2 = table(f2$lineage, f2$dst_mode); h2 +length(h) +length(h2) + + +f2 = test2[test2$mutationinformation == "Y95D",] +h = table(f2$lineage, f2$dst); h +h2 = table(f2$lineage, f2$sensitivity); h2 +length(h) +length(h2) + +tm = "G97A" # 1 +tm = "L117R" +tm = "D63G" +tm = "A102P" +tm = "F13L" +tm = "E174G" +tm = "L182S" +tm = "L4S" + +f3 = test2[test2$mutationinformation == tm,] +h3 = table(f3$lineage, f3$sensitivity); h3 +print(h3) +print(class(h3)) +print(dim(h3)) +dim(h3)[1] # >1 +dim(h3)[2] #>1 +#h3 = table(f3$lineage); h3 +length(h3) + +h3v2 = table(f3$lineage, f3$sensitivity); h3v2 +length(h3v2) + +#if length is > 2, then get these +chisq.test(h3) +chisq.test(h3)$p.value + +#ggchisqtest(chisq.test(h3)) + +fisher.test(h3) +fisher.test(h3)$p.value + +######################### +muts = unique(my_df2$mutationinformation) +my_df = my_df2 + +# step1 : get muts with more than one lineage +lin_muts = NULL +for (i in muts) { + print (i) + s_mut = my_df[my_df$mutationinformation == i,] + s_tab = table(s_mut$lineage, s_mut$sensitivity) + #s_tab = table(s_mut$lineage) + #print(s_tab) + + #if (length(s_tab) > 1 ){ + # if (dim(s_tab)[1] > 1 ){ + # lin_muts = c(lin_muts, i) + if (dim(s_tab)[1] > 1 && dim(s_tab)[2] > 1){ + lin_muts = c(lin_muts, i) + + } +} + + +# # now from the above list, get only the ones that have both R and S +# muts_var = NULL +# for (i in lin_muts) { +# print (i) +# s_mut = my_df[my_df$mutationinformation == i,] +# s_tab = table(s_mut$lineage, s_mut$sensitivity) +# print(s_tab) +# print(dim(s_tab)[2]) # if this is one, we are uninterested +# if ( dim(s_tab)[2] > 1 ){ +# muts_var = c(muts_var, i) +# } +# } + + +# now final check +for (i in lin_muts) { + print (i) + s_mut = my_df[my_df$mutationinformation == i,] + s_tab = table(s_mut$lineage, s_mut$sensitivity) + print(s_tab) + print(c(i, "FT:", fisher.test(s_tab)$p.value)) + # print(dim(s_tab)[2]) # if this is one, we are uninterested + # if ( dim(s_tab)[2] > 1 ){ + # muts_var = c(muts_var, i) + # } + +} + +plot_df = my_df[my_df$mutationinformation%in%lin_muts,] + +#plot_df2 = plot_df[plot_df$lineage%in%sel_lineages,] + + + +table(plot_df$lineage) +length(unique(plot_df2$mutationinformation)) == length(lin_muts) + +#muts_var +lin_mutsL = plot_df$mutationinformation[plot_df$mutationinformation%in%lin_muts] + + +plot_df$p.value = NULL + +for (i in lin_muts) { + print (i) + s_mut = plot_df[plot_df$mutationinformation == i,] + print(s_mut) + s_tab = table(s_mut$lineage, s_mut$sensitivity) + print(s_tab) + ft_pvalue_i = round(fisher.test(s_tab)$p.value, 2) + + print(ft_pvalue_i) + + # #my_df[my_df['mutationinformation']==i,]['ft_pvalue']= ft_pvalue_i + #plot_df[plot_df['mutationinformation']==i,]['p.value']= ft_pvalue_i + + plot_df$p.value[plot_df$mutationinformation == i] <- ft_pvalue_i + #print(s_tab) + } + + + +plot_df2 = my_df[my_df$mutationinformation == c("A102P"),] +#https://stackoverflow.com/questions/72618364/how-to-use-geom-signif-from-ggpubr-with-a-chi-square-test + +######################### +library(grid) +#sp2 + annotation_custom(grob)+facet_wrap(~cyl, scales="free") +grob <- grobTree(textGrob("Scatter plot", x=0.1, y=0.95, hjust=0, + gp=gpar(col="red", fontsize=5, fontface="italic"))) + +############# +chi.test <- function(a, b) { + return(chisq.test(cbind(a, b))) +} + +ggplot(plot_df, aes(x = lineage + #, y = snp_frequency + , fill = factor(sensitivity))) + + geom_bar( + stat = 'count' + #stat = 'identity' + , position = 'dodge') + + facet_wrap(~mutationinformation + , scales = 'free_y') + + #coord_flip() + + stat_count(aes(y=..count../sum(..count..), label=p.value), geom="text", hjust=0) + + #geom_text(aes(label = p.value, x = -0.5, y = 1)) + + #geom_text(data = data.frame(lineage = c("L1", "L2", "L3", "L4"), p.value = "p.value" )) + #geom_text(aes(label = p.value), stat = "count") + + + #geom_text(aes(label=after_stat(count)), vjust=0, stat = "count") # shows numbers + + #geom_signif(comparisons = list(c("L1", "L2", "L3", "L4")), test = "fisher.test", y = 1) + + # geom_signif(data = data.frame(lineage = c("L1", "L2", "L3", "L4"),sensitivity = c("R", "S") ) + # , test = "fisher.test" ) + # , aes(y_position=c(5.3, 8.3), xmin=c(0.8, 0.8), xmax=c(1.2, 1.2)) + # ) + + + #geom_label(p.value) + #coord_flip() + # ggforce::facet_wrap_paginate(~mutationinformation + # , ncol = 5 + # , nrow = 5 + # , page = 10 + # ) + + + + +# with coord flip +ggplot(plot_data, aes(x = lineage_mm, fill = sensitivity)) + + geom_bar(position = 'dodge') + + facet_wrap(~mutationinformation) + coord_flip() + +#============ +# facet grid +#============ +ggplot(plot_data, aes(x = mutationinformation, fill = sensitivity)) + + geom_bar(position = 'dodge') + + facet_grid(~lineage_mm) + +# with coord flip +ggplot(plot_data, aes(x = mutationinformation, fill = sensitivity)) + + geom_bar(position = 'dodge') + + facet_grid(~lineage_mm)+ coord_flip() + +########################################## +#%% useful info +# https://stackoverflow.com/questions/13773770/split-comma-separated-strings-in-a-column-into-separate-rows +bardf = as.data.frame(bar) +class(bardf) == class(my_df) + +baz = my_df +baz = baz %>% + mutate(col2 = strsplit(as.character(col2), ",")) %>% + unnest(col2) +baz = as.data.frame(baz) +class(baz) == class(bar) +