#!/usr/bin/env python3 # -*- coding: utf-8 -*- ''' Created on Tue Aug 6 12:56:03 2019 @author: tanu ''' #======================================================================= # Task: combining all dfs to a single one # Input: 8 dfs #1) .lower()'_complex_mcsm_norm.csv' #2) .lower()_foldx.csv' #3) .lower()_dssp.csv' #4) .lower()_kd.csv' #5) .lower()_rd.csv' #6) 'ns' + .lower()_snp_info.csv' #7) .lower()_af_or.csv' #8) .lower() _af_or_kinship.csv # combining order #Merge1 = 1 + 2 #Merge2 = 3 + 4 #Merge3 = Merge2 + 5 #Merge4 = Merge1 + Merge3 #Merge5 = 6 + 7 #Merge6 = Merge5 + 8 #Merge7 = Merge4 + Merge6 # Output: single csv of all 8 dfs combined # useful link # https://stackoverflow.com/questions/23668427/pandas-three-way-joining-multiple-dataframes-on-columns #%% FIXME: let the script proceed even if files don't exist! # i.e example below # '/home/tanu/git/Data/ethambutol/output/dynamut_results/embb_complex_dynamut_norm.csv' #======================================================================= #%% load packages import sys, os import pandas as pd from pandas import DataFrame import numpy as np import argparse from functools import reduce #======================================================================= #%% specify input and curr dir homedir = os.path.expanduser('~') # set working dir os.getcwd() #os.chdir(homedir + '/git/LSHTM_analysis/scripts') os.getcwd() # FIXME: local imports #from combining import combine_dfs_with_checks from combining_FIXME import detect_common_cols from reference_dict import oneletter_aa_dict from reference_dict import low_3letter_dict from aa_code import get_aa_3lower from aa_code import get_aa_1upper # REGEX: as required # mcsm_regex = r'^([A-Za-z]{1})([0-9]+)([A-Za-z]{1})$' # mcsm_wt = mcsm_df['mutationinformation'].str.extract(mcsm_regex)[0] # mcsm_mut = mcsm_df['mutationinformation'].str.extract(mcsm_regex)[2] # gwas_regex = r'^([A-Za-z]{3})([0-9]+)([A-Za-z]{3})$' # gwas_wt = mcsm_df['mutation'].str.extract(gwas_regex)[0] # gwas_pos = mcsm_df['mutation'].str.extract(gwas_regex)[1] # gwas_mut = mcsm_df['mutation'].str.extract(gwas_regex)[2] #======================================================================= #%% command line args: case sensitive arg_parser = argparse.ArgumentParser() arg_parser.add_argument('-d', '--drug', help = 'drug name', default = '') arg_parser.add_argument('-g', '--gene', help = 'gene name', default = '') arg_parser.add_argument('--datadir', help = 'Data Directory. By default, it assmumes homedir + git/Data') arg_parser.add_argument('-i', '--input_dir', help = 'Input dir containing pdb files. By default, it assmumes homedir + + input') arg_parser.add_argument('-o', '--output_dir', help = 'Output dir for results. By default, it assmes homedir + + output') arg_parser.add_argument('--debug', action ='store_true', help = 'Debug Mode') args = arg_parser.parse_args() #======================================================================= #%% variable assignment: input and output drug = args.drug gene = args.gene datadir = args.datadir indir = args.input_dir outdir = args.output_dir gene_match = gene + '_p.' print('mut pattern for gene', gene, ':', gene_match) #%%======================================================================= #============== # directories #============== if not datadir: datadir = homedir + '/git/Data/' if not indir: indir = datadir + drug + '/input/' if not outdir: outdir = datadir + drug + '/output/' #======= # input #======= gene_list_normal = ["pnca", "katg", "rpob", "alr"] #FIXME: for gid, this should be SRY as this is the drug...please check!!!! if gene.lower() == "gid": print("\nReading mCSM file for gene:", gene) in_filename_mcsm = gene.lower() + '_complex_mcsm_norm_SRY.csv' # was incorrectly SAM previously if gene.lower() == "embb": print("\nReading mCSM file for gene:", gene) #in_filename_mcsm = gene.lower() + '_complex_mcsm_norm1.csv' #798 in_filename_mcsm = gene.lower() + '_complex_mcsm_norm2.csv' #844 if gene.lower() in gene_list_normal: print("\nReading mCSM file for gene:", gene) in_filename_mcsm = gene.lower() + '_complex_mcsm_norm.csv' infile_mcsm = outdir + in_filename_mcsm mcsm_df = pd.read_csv(infile_mcsm, sep = ',') in_filename_foldx = gene.lower() + '_foldx.csv' infile_foldx = outdir + in_filename_foldx foldx_df = pd.read_csv(infile_foldx , sep = ',') in_filename_deepddg = gene.lower() + '_ni_deepddg.csv' # change to decent filename and put it in the correct dir infile_deepddg = outdir + in_filename_deepddg deepddg_df = pd.read_csv(infile_deepddg, sep = ',') in_filename_dssp = gene.lower() + '_dssp.csv' infile_dssp = outdir + in_filename_dssp dssp_df_raw = pd.read_csv(infile_dssp, sep = ',') in_filename_kd = gene.lower() + '_kd.csv' infile_kd = outdir + in_filename_kd kd_df = pd.read_csv(infile_kd, sep = ',') in_filename_rd = gene.lower() + '_rd.csv' infile_rd = outdir + in_filename_rd rd_df = pd.read_csv(infile_rd, sep = ',') #in_filename_snpinfo = 'ns' + gene.lower() + '_snp_info_f.csv' # gwas f info #infile_snpinfo = outdir + in_filename_snpinfo in_filename_afor = gene.lower() + '_af_or.csv' infile_afor = outdir + in_filename_afor afor_df = pd.read_csv(infile_afor, sep = ',') #in_filename_afor_kin = gene.lower() + '_af_or_kinship.csv' #infile_afor_kin = outdir + in_filename_afor_kin infilename_dynamut2 = gene.lower() + '_dynamut2_norm.csv' infile_dynamut2 = outdir + 'dynamut_results/dynamut2/' + infilename_dynamut2 dynamut2_df = pd.read_csv(infile_dynamut2, sep = ',') infilename_mcsm_f_snps = gene.lower() + '_mcsm_formatted_snps.csv' infile_mcsm_f_snps = outdir + infilename_mcsm_f_snps mcsm_f_snps = pd.read_csv(infile_mcsm_f_snps, sep = ',', names = ['mutationinformation'], header = None) #------------------------------------------------------------------------------ # ONLY:for gene pnca and gid: End logic should pick this up! geneL_dy_na = ['gid'] if gene.lower() in geneL_dy_na : print("\nGene:", gene.lower() , "\nReading Dynamut and mCSM_na files") infilename_dynamut = gene.lower() + '_dynamut_norm.csv' # gid infile_dynamut = outdir + 'dynamut_results/' + infilename_dynamut dynamut_df = pd.read_csv(infile_dynamut, sep = ',') infilename_mcsm_na = gene.lower() + '_complex_mcsm_na_norm.csv' # gid infile_mcsm_na = outdir + 'mcsm_na_results/' + infilename_mcsm_na mcsm_na_df = pd.read_csv(infile_mcsm_na, sep = ',') # ONLY:for gene embb and alr: End logic should pick this up! geneL_ppi2 = ['embb', 'alr'] #if gene.lower() == "embb" or "alr": if gene.lower() in geneL_ppi2: infilename_mcsm_ppi2 = gene.lower() + '_complex_mcsm_ppi2_norm.csv' infile_mcsm_ppi2 = outdir + 'mcsm_ppi2/' + infilename_mcsm_ppi2 mcsm_ppi2_df = pd.read_csv(infile_mcsm_ppi2, sep = ',') if gene.lower() == "embb": sel_chain = "B" else: sel_chain = "A" #------------------------------------------------------------------------------ #======= # output #======= out_filename_comb = gene.lower() + '_all_params.csv' outfile_comb = outdir + out_filename_comb print('\nOutput filename:', outfile_comb , '\n===================================================================') # end of variable assignment for input and output files #%%############################################################################ #===================== # some preprocessing #===================== #=========== # FoldX #=========== foldx_df.shape #---------------------- # scale foldx values #---------------------- # rename ddg column to ddg_foldx foldx_df['ddg'] foldx_df = foldx_df.rename(columns = {'ddg':'ddg_foldx'}) foldx_df['ddg_foldx'] # Rescale values in Foldx_change col b/w -1 and 1 so negative numbers # stay neg and pos numbers stay positive foldx_min = foldx_df['ddg_foldx'].min() foldx_max = foldx_df['ddg_foldx'].max() foldx_min foldx_max # quick check len(foldx_df.loc[foldx_df['ddg_foldx'] >= 0]) len(foldx_df.loc[foldx_df['ddg_foldx'] < 0]) foldx_scale = lambda x : x/abs(foldx_min) if x < 0 else (x/foldx_max if x >= 0 else 'failed') foldx_df['foldx_scaled'] = foldx_df['ddg_foldx'].apply(foldx_scale) print('\nRaw foldx scores:\n', foldx_df['ddg_foldx'] , '\n---------------------------------------------------------------' , '\nScaled foldx scores:\n', foldx_df['foldx_scaled']) # additional check added fsmi = foldx_df['foldx_scaled'].min() fsma = foldx_df['foldx_scaled'].max() c = foldx_df[foldx_df['ddg_foldx']>=0].count() foldx_pos = c.get(key = 'ddg_foldx') c2 = foldx_df[foldx_df['foldx_scaled']>=0].count() foldx_pos2 = c2.get(key = 'foldx_scaled') if foldx_pos == foldx_pos2 and fsmi == -1 and fsma == 1: print('\nPASS: Foldx values scaled correctly b/w -1 and 1') else: print('\nFAIL: Foldx values scaled numbers MISmatch' , '\nExpected number:', foldx_pos , '\nGot:', foldx_pos2 , '\n======================================================') #------------------------- # foldx outcome category: # Remember, its inverse # +ve: Destabilising # -ve: Stabilising #-------------------------- foldx_df['foldx_outcome'] = foldx_df['ddg_foldx'].apply(lambda x: 'Destabilising' if x >= 0 else 'Stabilising') foldx_df[foldx_df['ddg_foldx']>=0].count() foc = foldx_df['foldx_outcome'].value_counts() if foc['Destabilising'] == foldx_pos and foc['Destabilising'] == foldx_pos2: print('\nPASS: Foldx outcome category created') else: print('\nFAIL: Foldx outcome category could NOT be created' , '\nExpected number:', foldx_pos , '\nGot:', foc[0] , '\n======================================================') sys.exit() #======================= # Deepddg # TODO: RERUN 'gid' #======================= deepddg_df.shape #-------------------------- # check if >1 chain #-------------------------- deepddg_df.loc[:,'chain_id'].value_counts() if len(deepddg_df.loc[:,'chain_id'].value_counts()) > 1: print("\nChains detected: >1" , "\nGene:", gene , "\nChains:", deepddg_df.loc[:,'chain_id'].value_counts().index) print('\nSelecting chain:', sel_chain, 'for gene:', gene) deepddg_df = deepddg_df[deepddg_df['chain_id'] == sel_chain] #-------------------------- # Check for duplicates #-------------------------- if len(deepddg_df['mutationinformation'].duplicated().value_counts())> 1: print("\nFAIL: Duplicates detected in DeepDDG infile" , "\nNo. of duplicates:" , deepddg_df['mutationinformation'].duplicated().value_counts()[1] , "\nformat deepDDG infile before proceeding") sys.exit() else: print("\nPASS: No duplicates detected in DeepDDG infile") #-------------------------- # Drop chain id col as other targets don't have it.Check for duplicates #-------------------------- col_to_drop = ['chain_id'] deepddg_df = deepddg_df.drop(col_to_drop, axis = 1) #------------------------- # scale Deepddg values #------------------------- # Rescale values in deepddg_change col b/w -1 and 1 so negative numbers # stay neg and pos numbers stay positive deepddg_min = deepddg_df['deepddg'].min() deepddg_max = deepddg_df['deepddg'].max() deepddg_scale = lambda x : x/abs(deepddg_min) if x < 0 else (x/deepddg_max if x >= 0 else 'failed') deepddg_df['deepddg_scaled'] = deepddg_df['deepddg'].apply(deepddg_scale) print('\nRaw deepddg scores:\n', deepddg_df['deepddg'] , '\n---------------------------------------------------------------' , '\nScaled deepddg scores:\n', deepddg_df['deepddg_scaled']) # additional check added dsmi = deepddg_df['deepddg_scaled'].min() dsma = deepddg_df['deepddg_scaled'].max() c = deepddg_df[deepddg_df['deepddg']>=0].count() deepddg_pos = c.get(key = 'deepddg') c2 = deepddg_df[deepddg_df['deepddg_scaled']>=0].count() deepddg_pos2 = c2.get(key = 'deepddg_scaled') if deepddg_pos == deepddg_pos2 and dsmi == -1 and dsma == 1: print('\nPASS: deepddg values scaled correctly b/w -1 and 1') else: print('\nFAIL: deepddg values scaled numbers MISmatch' , '\nExpected number:', deepddg_pos , '\nGot:', deepddg_pos2 , '\n======================================================') #-------------------------- # Deepddg outcome category #-------------------------- deepddg_df['deepddg_outcome'] = deepddg_df['deepddg'].apply(lambda x: 'Stabilising' if x >= 0 else 'Destabilising') deepddg_df[deepddg_df['deepddg']>=0].count() doc = deepddg_df['deepddg_outcome'].value_counts() if doc['Stabilising'] == deepddg_pos and doc['Stabilising'] == deepddg_pos2: print('\nPASS: Deepddg outcome category created') else: print('\nFAIL: Deepddg outcome category could NOT be created' , '\nExpected number:', deepddg_pos , '\nGot:', doc[0] , '\n======================================================') sys.exit() if deepddg_df['deepddg_scaled'].min() == -1 and deepddg_df['deepddg_scaled'].max() == 1: print('\nPASS: Deepddg data is scaled between -1 and 1', '\nproceeding with merge') #%%============================================================================ # Now merges begin print('===================================' , '\nFirst merge: mcsm + foldx' , '\n===================================') mcsm_df.shape # add 3 lowercase aa code for wt and mutant get_aa_3lower(df = mcsm_df , wt_colname = 'wild_type' , mut_colname = 'mutant_type' , col_wt = 'wt_aa_3lower' , col_mut = 'mut_aa_3lower') #mcsm_df.columns = mcsm_df.columns.str.lower() # foldx_df.shape #mcsm_foldx_dfs = combine_dfs_with_checks(mcsm_df, foldx_df, my_join = "outer") merging_cols_m1 = detect_common_cols(mcsm_df, foldx_df) mcsm_foldx_dfs = pd.merge(mcsm_df , foldx_df , on = merging_cols_m1 , how = "outer") ncols_m1 = len(mcsm_foldx_dfs.columns) print('\n\nResult of first merge:', mcsm_foldx_dfs.shape , '\n===================================================================') mcsm_foldx_dfs[merging_cols_m1].apply(len) mcsm_foldx_dfs[merging_cols_m1].apply(len) == len(mcsm_foldx_dfs) #%% for embB and any other targets where mCSM-lig hasn't run for ALL nsSNPs. # Get the empty cells to be full of meaningful info if mcsm_foldx_dfs.loc[:,'wild_type': 'mut_aa_3lower'].isnull().values.any(): print ('\nNAs detected in mcsm cols after merge.' , '\nCleaning data before merging deepddg_df') ############################## # Extract relevant col values # code to one ############################## # wt_reg = r'(^[A-Z]{1})' # print('wild_type:', wt_reg) # mut_reg = r'[0-9]+(\w{1})$' # print('mut type:', mut_reg) mcsm_foldx_dfs['wild_type'] = mcsm_foldx_dfs.loc[:,'mutationinformation'].str.extract(r'(^[A-Z]{1})') mcsm_foldx_dfs['position'] = mcsm_foldx_dfs.loc[:,'mutationinformation'].str.extract(r'([0-9]+)') mcsm_foldx_dfs['mutant_type'] = mcsm_foldx_dfs.loc[:,'mutationinformation'].str.extract(r'[0-9]+([A-Z]{1})$') # BEWARE: Bit of logic trap i.e if nan comes first # in chain column, then nan will be populated! #df['foo'] = df['chain'].unique()[0] mcsm_foldx_dfs['chain'] = np.where(mcsm_foldx_dfs[['chain']].isnull().all(axis=1) , mcsm_foldx_dfs['chain'].unique()[0] , mcsm_foldx_dfs['chain']) mcsm_foldx_dfs['ligand_id'] = np.where(mcsm_foldx_dfs[['ligand_id']].isnull().all(axis=1) , mcsm_foldx_dfs['ligand_id'].unique()[0] , mcsm_foldx_dfs['ligand_id']) #-------------------------------------------------------------------------- mcsm_foldx_dfs['wild_pos'] = mcsm_foldx_dfs.loc[:,'wild_type'] + mcsm_foldx_dfs.loc[:,'position'].astype(int).apply(str) mcsm_foldx_dfs['wild_chain_pos'] = mcsm_foldx_dfs.loc[:,'wild_type'] + mcsm_foldx_dfs.loc[:,'chain'] + mcsm_foldx_dfs.loc[:,'position'].astype(int).apply(str) ############# # Map 1 letter # code to 3Upper ############# # initialise a sub dict that is lookup dict for # 3-LETTER aa code to 1-LETTER aa code lookup_dict = dict() for k, v in oneletter_aa_dict.items(): lookup_dict[k] = v['three_letter_code_lower'] wt = mcsm_foldx_dfs['wild_type'].squeeze() # converts to a series that map works on mcsm_foldx_dfs['wt_aa_3lower'] = wt.map(lookup_dict) mut = mcsm_foldx_dfs['mutant_type'].squeeze() mcsm_foldx_dfs['mut_aa_3lower'] = mut.map(lookup_dict) else: print('\nNo NAs detected in mcsm_fold_dfs. Proceeding to merge deepddg_df') #%% print('===================================' , '\nSecond merge: mcsm_foldx_dfs + deepddg' , '\n===================================') # merge with mcsm_foldx_dfs and deepddg_df mcsm_foldx_deepddg_dfs = pd.merge(mcsm_foldx_dfs , deepddg_df , on = 'mutationinformation' , how = "left") mcsm_foldx_deepddg_dfs['deepddg_outcome'].value_counts() ncols_deepddg_merge = len(mcsm_foldx_deepddg_dfs.columns) mcsm_foldx_deepddg_dfs['position'] = mcsm_foldx_deepddg_dfs['position'].astype('int64') #%%============================================================================ #FIXME: select df with 'chain' to allow corret dim merging! print('===================================' , '\nThird merge: dssp + kd' , '\n===================================') dssp_df_raw.shape kd_df.shape rd_df.shape dssp_df = dssp_df_raw[dssp_df_raw['chain_id'] == sel_chain] dssp_df['chain_id'].value_counts() #dssp_kd_dfs = combine_dfs_with_checks(dssp_df, kd_df, my_join = "outer") merging_cols_m2 = detect_common_cols(dssp_df, kd_df) dssp_kd_dfs = pd.merge(dssp_df , kd_df , on = merging_cols_m2 , how = "outer") print('\n\nResult of third merge:', dssp_kd_dfs.shape , '\n===================================================================') #%%============================================================================ print('===================================' , '\nFourth merge: third merge + rd_df' , '\ndssp_kd_dfs + rd_df' , '\n===================================') #dssp_kd_rd_dfs = combine_dfs_with_checks(dssp_kd_dfs, rd_df, my_join = "outer") merging_cols_m3 = detect_common_cols(dssp_kd_dfs, rd_df) dssp_kd_rd_dfs = pd.merge(dssp_kd_dfs , rd_df , on = merging_cols_m3 , how = "outer") ncols_m3 = len(dssp_kd_rd_dfs.columns) print('\n\nResult of Third merge:', dssp_kd_rd_dfs.shape , '\n===================================================================') dssp_kd_rd_dfs[merging_cols_m3].apply(len) dssp_kd_rd_dfs[merging_cols_m3].apply(len) == len(dssp_kd_rd_dfs) #%%============================================================================ print('=======================================' , '\nFifth merge: Second merge + fourth merge' , '\nmcsm_foldx_dfs + dssp_kd_rd_dfs' , '\n=======================================') #combined_df = combine_dfs_with_checks(mcsm_foldx_dfs, dssp_kd_rd_dfs, my_join = "inner") #merging_cols_m4 = detect_common_cols(mcsm_foldx_dfs, dssp_kd_rd_dfs) #combined_df = pd.merge(mcsm_foldx_dfs, dssp_kd_rd_dfs, on = merging_cols_m4, how = "inner") #combined_df_expected_cols = ncols_m1 + ncols_m3 - len(merging_cols_m4) # with deepddg values merging_cols_m4 = detect_common_cols(mcsm_foldx_deepddg_dfs, dssp_kd_rd_dfs) combined_df = pd.merge(mcsm_foldx_deepddg_dfs , dssp_kd_rd_dfs , on = merging_cols_m4 , how = "inner") combined_df_expected_cols = ncols_deepddg_merge + ncols_m3 - len(merging_cols_m4) # FIXME: check logic, doesn't effect anything else! if not gene == "embB": print("\nGene is:", gene) if len(combined_df) == len(mcsm_df) and len(combined_df.columns) == combined_df_expected_cols: print('PASS: successfully combined 5 dfs' , '\nNo. of rows combined_df:', len(combined_df) , '\nNo. of cols combined_df:', len(combined_df.columns)) else: #sys.exit('FAIL: check individual df merges') print("\nGene is:", gene , "\ncombined_df length:", len(combined_df) , "\nmcsm_df_length:", len(mcsm_df) ) if len(combined_df.columns) == combined_df_expected_cols: print('PASS: successfully combined 5 dfs' , '\nNo. of rows combined_df:', len(combined_df) , '\nNo. of cols combined_df:', len(combined_df.columns)) else: sys.exit('FAIL: check individual merges') print('\nResult of Fourth merge:', combined_df.shape , '\n===================================================================') combined_df[merging_cols_m4].apply(len) combined_df[merging_cols_m4].apply(len) == len(combined_df) #%%============================================================================ # Format the combined df columns combined_df_colnames = combined_df.columns # check redundant columns combined_df['chain'].equals(combined_df['chain_id']) combined_df['wild_type'].equals(combined_df['wild_type_kd']) # has nan combined_df['wild_type'].equals(combined_df['wild_type_dssp']) # sanity check foo = combined_df[['wild_type', 'wild_type_kd', 'wt_3letter_caps', 'wt_aa_3lower', 'mut_aa_3lower']] # Drop cols cols_to_drop = ['chain_id', 'wild_type_kd', 'wild_type_dssp', 'wt_3letter_caps'] combined_df_clean = combined_df.drop(cols_to_drop, axis = 1) del(foo) #%%============================================================================ # Output columns out_filename_stab_struc = gene.lower() + '_comb_stab_struc_params.csv' outfile_stab_struc = outdir + out_filename_stab_struc print('Output filename:', outfile_stab_struc , '\n===================================================================') combined_df_clean # write csv print('\nWriting file: combined stability and structural parameters') combined_df_clean.to_csv(outfile_stab_struc, index = False) print('\nFinished writing file:' , '\nNo. of rows:', combined_df_clean.shape[0] , '\nNo. of cols:', combined_df_clean.shape[1]) #%%===================================================================== print('\n=======================================' , '\nFifth merge:' , '\ncombined_df_clean + afor_df ' , '\n=======================================') afor_cols = afor_df.columns afor_df.shape # create a mapping from the gwas mutation column i.e _abcXXXrst #---------------------- # call get_aa_upper(): # adds 3 more cols with one letter aa code #---------------------- get_aa_1upper(df = afor_df , gwas_mut_colname = 'mutation' , wt_colname = 'wild_type' , pos_colname = 'position' , mut_colname = 'mutant_type') afor_df['mutationinformation'] = afor_df['wild_type'] + afor_df['position'].map(str) + afor_df['mutant_type'] afor_cols = afor_df.columns merging_cols_m5 = detect_common_cols(combined_df_clean, afor_df) # remove position so that merging can take place without dtype conflicts merging_cols_m5.remove('position') # drop position column from afor_df afor_df = afor_df.drop(['position'], axis = 1) afor_cols = afor_df.columns # merge combined_stab_afor = pd.merge(combined_df_clean , afor_df , on = merging_cols_m5 , how = "left") comb_afor_df_cols = combined_stab_afor.columns comb_afor_expected_cols = len(combined_df_clean.columns) + len(afor_df.columns) - len(merging_cols_m5) if len(combined_stab_afor) == len(combined_df_clean) and len(combined_stab_afor.columns) == comb_afor_expected_cols: print('\nPASS: successfully combined 6 dfs' , '\nNo. of rows combined_stab_afor:', len(combined_stab_afor) , '\nNo. of cols combined_stab_afor:', len(combined_stab_afor.columns)) else: sys.exit('\nFAIL: check individual df merges') print('\n\nResult of Fifth merge:', combined_stab_afor.shape , '\n===================================================================') combined_stab_afor[merging_cols_m5].apply(len) combined_stab_afor[merging_cols_m5].apply(len) == len(combined_stab_afor) if (len(combined_stab_afor) - combined_stab_afor['mutation'].isna().sum()) == len(afor_df): print('\nPASS: Merge successful for af and or with matched numbers') if len(combined_stab_afor) - combined_stab_afor['mutation'].isna().sum() == len(afor_df)-len(afor_df[~afor_df['mutation'].isin(combined_stab_afor['mutation'])]): print("\nMismatched numbers, OR df has extra snps not found in mcsm df" , "\nNo. of nsSNPs with valid ORs:", len(afor_df) , "\nNo. of mcsm nsSNPs: ", len(combined_df_clean) , "\nNo. of OR nsSNPs not in mCSM df:" , len(afor_df[~afor_df['mutation'].isin(combined_stab_afor['mutation'])]) , "\nWriting these mutations to file:") orsnps_notmcsm = afor_df[~afor_df['mutation'].isin(combined_stab_afor['mutation'])] else: sys.exit('\nFAIL: merge unsuccessful for af and or') #%%============================================================================ # Output columns: when dynamut, dynamut2 and others weren't being combined out_filename_comb_afor = gene.lower() + '_comb_afor.csv' outfile_comb_afor = outdir + out_filename_comb_afor print('Output filename:', outfile_comb_afor , '\n===================================================================') # write csv print('Writing file: combined stability and afor') combined_stab_afor.to_csv(outfile_comb_afor, index = False) print('\nFinished writing file:' , '\nNo. of rows:', combined_stab_afor.shape[0] , '\nNo. of cols:', combined_stab_afor.shape[1]) #%%============================================================================ # combine dynamut, dynamut2, and mcsm_na #dfs_list = [dynamut_df, dynamut2_df, mcsm_na_df] # gid if gene.lower() == "pnca": dfs_list = [dynamut2_df] if gene.lower() == "gid": dfs_list = [dynamut_df, dynamut2_df, mcsm_na_df] if gene.lower() == "embb": dfs_list = [dynamut2_df, mcsm_ppi2_df] if gene.lower() == "katg": dfs_list = [dynamut2_df] if gene.lower() == "rpob": dfs_list = [dynamut2_df] if gene.lower() == "alr": dfs_list = [dynamut2_df, mcsm_ppi2_df] dfs_merged = reduce(lambda left,right: pd.merge(left , right , on = ['mutationinformation'] , how = 'inner') , dfs_list) # drop excess columns drop_cols = detect_common_cols(dfs_merged, combined_stab_afor) drop_cols.remove('mutationinformation') dfs_merged_clean = dfs_merged.drop(drop_cols, axis = 1) merging_cols_m6 = detect_common_cols(dfs_merged_clean, combined_stab_afor) len(dfs_merged_clean.columns) len(combined_stab_afor.columns) combined_all_params = pd.merge(combined_stab_afor , dfs_merged_clean , on = merging_cols_m6 , how = "inner") expected_ncols = len(dfs_merged_clean.columns) + len(combined_stab_afor.columns) - len(merging_cols_m6) expected_nrows = len(combined_stab_afor) if len(combined_all_params.columns) == expected_ncols and len(combined_all_params) == expected_nrows: print('\nPASS: All dfs combined') else: print('\nFAIL:lengths mismatch' , '\nExpected ncols:', expected_ncols , '\nGot:', len(dfs_merged_clean.columns) , '\nExpected nrows:', expected_nrows , '\nGot:', len(dfs_merged_clean) ) #%% Done for gid on 10/09/2021 # write csv print('Writing file: all params') combined_all_params.to_csv(outfile_comb, index = False) print('\nFinished writing file:' , '\nNo. of rows:', combined_all_params.shape[0] , '\nNo. of cols:', combined_all_params.shape[1]) #%% end of script