#!/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: 12/13/14 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 (ONLY for pncA, but omitted for the final run) #9) .lower()'_dynamut2.csv' #10) .lower()'_dynamut.csv' #11) .lower()'_mcsm_na.csv' #12) .lower()'_mcsm_ppi2.csv' #13) .lower()'_consurf.csv' #14) .lower()'_snap2.csv' # combining order # 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') sys.path.append(homedir + '/git/LSHTM_analysis/scripts') os.getcwd() #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'] 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 in_filename_mcsm = gene.lower() + '_complex_mcsm_norm3.csv' #851 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) # more output added ## consurf [change colnames] infilename_consurf = gene.lower() + '_consurf_grades_f.csv' infile_consurf = outdir + 'consurf/'+ infilename_consurf consurf_df = pd.read_csv(infile_consurf, sep = ',') ## SNAP2 [add normalised score] infilename_snap2 = gene.lower() + '_snap2_output.csv' infile_snap2 = outdir + 'snap2/'+ infilename_snap2 snap2_df = pd.read_csv(infile_snap2, sep = ',') #------------------------------------------------------------------------------ # ONLY: for gene 'gid' and 'rpob': End logic should pick this up! geneL_na = ['gid', 'rpob'] if gene.lower() in geneL_na: print("\nGene:", gene.lower() , "\nReading mCSM_na files") 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 = ',') geneL_dy = ['gid'] if gene.lower() in geneL_dy: print("\nGene:", gene.lower() , "\nReading Dynamut files") infilename_dynamut = gene.lower() + '_dynamut_norm.csv' # gid infile_dynamut = outdir + 'dynamut_results/' + infilename_dynamut dynamut_df = pd.read_csv(infile_dynamut, sep = ',') # ONLY: for genes 'alr', 'embb', 'katg' and 'rpob': End logic should pick this up! geneL_ppi2 = ['alr', 'embb', 'katg', 'rpob'] 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 #======= # outfile 3 out_filename_comb = gene.lower() + '_all_params.csv' outfile_comb = outdir + out_filename_comb # outfile 2 out_filename_comb_afor = gene.lower() + '_comb_afor.csv' outfile_comb_afor = outdir + out_filename_comb_afor # outfile 1 out_filename_stab_struc = gene.lower() + '_comb_stab_struc_params.csv' outfile_stab_struc = outdir + out_filename_stab_struc # end of variable assignment for input and output files #%%############################################################################ #===================== # some preprocessing #===================== #=========== # KD #=========== kd_df.shape # geneL_kd = ['alr'] # if gene.lower() in geneL_kd: # print('\nRunning gene:', gene.lower() # ,'\nChecking start numbering') if kd_df['wild_type_kd'].str.contains('X').any(): print('\nDetected X in wild_type_kd' , '\nRunning gene:', gene.lower() , '\nChecking start numbering') kd_df = kd_df[~kd_df['wild_type_kd'].str.contains('X')] #=========== # 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 #======================= 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] #-------------------------- # Drop chain_id col as other # targets don't have it. #-------------------------- col_to_drop = ['chain_id'] deepddg_df = deepddg_df.drop(col_to_drop, axis = 1) #-------------------------- # 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") #------------------------- # 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======================================================') 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') #-------------------------- # Deepddg outcome category #-------------------------- if 'deepddg_outcome' not in deepddg_df.columns: print('\nCreating column: deepddg_outcome') 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() print(doc) else: print('\nColumn exists: deepddg_outcome') t1 = deepddg_df['deepddg_outcome'].value_counts() deepddg_df['deepddg_outcome2'] = deepddg_df['deepddg'].apply(lambda x: 'Stabilising' if x >= 0 else 'Destabilising') t2 = deepddg_df['deepddg_outcome2'].value_counts() print('\n', t1, '\n', t2) #-------------------------- # Drop deepddg_outcome2 col #-------------------------- col_to_drop2 = ['deepddg_outcome2'] deepddg_df = deepddg_df.drop(col_to_drop2, axis = 1) if all(t1 == t2): print('\nPASS: Deepddg_outcome category checked!') doc = deepddg_df['deepddg_outcome'].value_counts() else: print('\nMISmatch in deepddg_outcome counts' , '\n:', t1 , '\n:', t2) 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() #======================= # Consurf #======================= consurf_df.shape #---------------------- # rename colums #---------------------- consurf_df.columns print('\nRenaming cols and assigning pretty column names') geneL_consurf = ['alr', 'katg', 'rpob'] if gene.lower() in geneL_consurf: consurf_df = consurf_df.rename(columns={'POS' : 'position_consurf'}) #--------------------------- # Specify the offset #--------------------------- print('\nAdding offset value for gene:', gene.lower()) if gene.lower() == 'alr': offset_val = 34 print('\nUsing offset val:', offset_val) if gene.lower() == 'katg': offset_val = 23 print('\nUsing offset val:', offset_val) if gene.lower() == 'rpob': offset_val = 28 print('\nUsing offset val:', offset_val) consurf_df['position'] = consurf_df['position_consurf'] + offset_val else: consurf_df = consurf_df.rename(columns={'POS' : 'position'}) consurf_df = consurf_df.rename(columns={'SEQ' : 'wild_type' , '3LATOM': 'wt_3upper' , 'SCORE' : 'consurf_score' , 'COLOR' : 'consurf_colour_str' , 'CONFIDENCEINTERVAL' : 'consurf_ci' , 'CONFIDENCEINTERVALCOLORS' : 'consurf_ci_colour' , 'MSADATA' : 'consurf_msa_data' , 'RESIDUEVARIETY' : 'consurf_aa_variety'}) # quick check if len(consurf_df) == len(kd_df): print('\nPASS: length of consurf df is as expected' , '\nProceeding to format consurf df') else: print('\nFAIL: length mismatch' , '\nExpected nrows:', len(rd_df) , '\nGot:', len(consurf_df)) consurf_df.dtypes consurf_df['consurf_score'] = consurf_df['consurf_score'].astype(float) consurf_df['consurf_colour'] = consurf_df['consurf_colour_str'].str.extract(r'(\d).*') consurf_df['consurf_colour'] = consurf_df['consurf_colour'].astype(int) consurf_df['consurf_colour_rev'] = consurf_df['consurf_colour_str'].str.replace(r'.\*','0') # non struc position are assigned a *, replacing that with a 0 so its all integer consurf_df['consurf_colour_rev'] = consurf_df['consurf_colour_rev'].astype(int) consurf_df['consurf_ci_upper'] = consurf_df['consurf_ci'].str.extract(r'(.*):') consurf_df['consurf_ci_upper'] = consurf_df['consurf_ci_upper'].astype(float) consurf_df['consurf_ci_lower'] = consurf_df['consurf_ci'].str.extract(r':(.*)') consurf_df['consurf_ci_lower'] = consurf_df['consurf_ci_lower'].astype(float) #consurf_df['wt_3upper_f'] = consurf_df['wt_3upper'].str.extract(r'^\w{3}(\d+.*)') #consurf_df['wt_3upper_f'] consurf_df['chain'] = consurf_df['wt_3upper'].str.extract(r':(.*)') consurf_df['wt_3upper'] = consurf_df['wt_3upper'].str.replace(r'(\d+:.*)', '') #------------------------- # scale consurf values #------------------------- # Rescale values in consurf_score col b/w -1 and 1 so negative numbers # stay neg and pos numbers stay positive consurf_min = consurf_df['consurf_score'].min() consurf_max = consurf_df['consurf_score'].max() consurf_min consurf_max # quick check len(consurf_df.loc[consurf_df['consurf_score'] >= 0]) len(consurf_df.loc[consurf_df['consurf_score'] < 0]) consurf_scale = lambda x : x/abs(consurf_min) if x < 0 else (x/consurf_max if x >= 0 else 'failed') consurf_df['consurf_scaled'] = consurf_df['consurf_score'].apply(consurf_scale) print('\nRaw consurf scores:\n', consurf_df['consurf_score'] , '\n---------------------------------------------------------------' , '\nScaled consurf scores:\n', consurf_df['consurf_scaled']) # additional check added csmi = consurf_df['consurf_scaled'].min() csma = consurf_df['consurf_scaled'].max() c = consurf_df[consurf_df['consurf_score']>=0].count() consurf_pos = c.get(key = 'consurf_score') c2 = consurf_df[consurf_df['consurf_scaled']>=0].count() consurf_pos2 = c2.get(key = 'consurf_scaled') if consurf_pos == consurf_pos2 and csmi == -1 and csma == 1: print('\nPASS: Consurf values scaled correctly b/w -1 and 1') else: print('\nFAIL: Consurf values scaled numbers MISmatch' , '\nExpected number:', consurf_pos , '\nGot:', consurf_pos2 , '\n======================================================') consurf_df.dtypes consurf_df.columns #--------------------------- # select columns # (and also determine order) # this removes redundant cols: # consurf_colour_str # consurf_ci #--------------------------- consurf_col_order = ['position' , 'wild_type' , 'chain' , 'wt_3upper' , 'consurf_score' , 'consurf_scaled' , 'consurf_colour' , 'consurf_colour_rev' , 'consurf_ci_upper' , 'consurf_ci_lower' , 'consurf_ci_colour' , 'consurf_msa_data' , 'consurf_aa_variety'] consurf_df_f = consurf_df[consurf_col_order] # CHECK: whether a general rule or a gene specific rule! if consurf_df_f['chain'].isna().sum() > 0: print('\nNaN detected in column chain for consurf df') #if gene.lower() == 'embb': print('\nFurther consurf df processing for gene:', gene.lower()) print('\nDropping Nan from column name chain') consurf_df_f = consurf_df_f[consurf_df_f['chain'].notna()] #======================= # SNAP2 #======================= snap2_df.shape #---------------------- # rename colums #---------------------- geneL_snap2 = ['alr', 'katg', 'rpob'] if gene.lower() in geneL_snap2: print('\nReading SNAP2 for gene:', gene.lower() , '\nOffset column also being read' , '\nRenaming columns...' , '\nColumn mutationinformation exists. Renaming SNAP2 column variant --> mutationinformation') snap2_df = snap2_df.rename(columns = {'mutationinformation': 'mutationinformation' , 'Variant' : 'mutationinformation_snap2' , 'Predicted Effect' : 'snap2_outcome' , 'Score' : 'snap2_score' , 'Expected Accuracy': 'snap2_accuracy_pc'}) else: print('\nReading SNAP2 for gene:', gene.lower() , '\nNo offset column for SNAP2' , '\nRenaming columns...' , '\nRenaming SNAP2 column variant --> mutationinformation') snap2_df = snap2_df.rename(columns = {'Variant' : 'mutationinformation' , 'Predicted Effect' : 'snap2_outcome' , 'Score' : 'snap2_score' , 'Expected Accuracy': 'snap2_accuracy_pc'}) snap2_df.columns snap2_df.head() snap2_df.dtypes snap2_df['snap2_accuracy_pc'] = snap2_df['snap2_accuracy_pc'].str.replace('%','') snap2_df['snap2_accuracy_pc'] = snap2_df['snap2_accuracy_pc'].astype(int) #------------------------- # scale snap2 values #------------------------- # Rescale values in snap2_score col b/w -1 and 1 so negative numbers # stay neg and pos numbers stay positive snap2_min = snap2_df['snap2_score'].min() snap2_max = snap2_df['snap2_score'].max() snap2_min snap2_max # quick check len(snap2_df.loc[snap2_df['snap2_score'] >= 0]) len(snap2_df.loc[snap2_df['snap2_score'] < 0]) snap2_scale = lambda x : x/abs(snap2_min) if x < 0 else (x/snap2_max if x >= 0 else 'failed') snap2_df['snap2_scaled'] = snap2_df['snap2_score'].apply(snap2_scale) print('\nRaw snap2 scores:\n', snap2_df['snap2_score'] , '\n---------------------------------------------------------------' , '\nScaled snap2 scores:\n', snap2_df['snap2_scaled']) # additional check added ssmi = snap2_df['snap2_scaled'].min() ssma = snap2_df['snap2_scaled'].max() sn = snap2_df[snap2_df['snap2_score']>=0].count() snap2_pos = sn.get(key = 'snap2_score') sn2 = snap2_df[snap2_df['snap2_scaled']>=0].count() snap2_pos2 = sn2.get(key = 'snap2_scaled') if snap2_pos == snap2_pos2 and csmi == -1 and csma == 1: print('\nPASS: Snap2 values scaled correctly b/w -1 and 1') else: print('\nFAIL: snap2 values scaled numbers MISmatch' , '\nExpected number:', snap2_pos , '\nGot:', snap2_pos2 , '\n======================================================') #------------------------------------- # select columns # (and also determine order) # renumbering already done using # bash and corrected file is read in #------------------------------------- snap2_df.dtypes snap2_df.columns geneL_snap2 = ['alr', 'katg', 'rpob'] if gene.lower() in geneL_snap2: print('\nSelecting cols SNAP2 for gene:', gene.lower()) snap2_df_f = snap2_df[['mutationinformation' , 'mutationinformation_snap2' , 'snap2_score' , 'snap2_scaled' , 'snap2_accuracy_pc' , 'snap2_outcome']] else: print('\nSelecting cols SNAP2 for gene:', gene.lower()) snap2_df_f = snap2_df[['mutationinformation' , 'snap2_score' , 'snap2_scaled' , 'snap2_accuracy_pc' , 'snap2_outcome']] #%%============================================================================ # 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') #%%============================================================================ # 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") , how = "inner") 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('===================================' , '\nFourth merge*: fourth merge + consurf_df' , '\dssp_kd_rd_dfs + consurf_df' , '\n===================================') #dssp_kd_rd_dfs = combine_dfs_with_checks(dssp_kd_dfs, rd_df, my_join = "outer") merging_cols_m3_v2 = detect_common_cols(dssp_kd_rd_dfs, consurf_df) dssp_kd_rd_con_dfs = pd.merge(dssp_kd_rd_dfs , consurf_df , on = merging_cols_m3_v2 , how = "outer") ncols_m3_v2 = len(dssp_kd_rd_con_dfs.columns) print('\n\nResult of fourth merge*:', dssp_kd_rd_con_dfs.shape , '\n===================================================================') dssp_kd_rd_con_dfs[merging_cols_m3_v2].apply(len) dssp_kd_rd_con_dfs[merging_cols_m3_v2].apply(len) == len(dssp_kd_rd_con_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) # Check: whether logic effects 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) combined_df_clean.columns del(foo) #%%============================================================================ #--------------------- # Output 1: write csv #--------------------- print('\nWriting file: combined stability and structural parameters' , '\nOutput 1 filename:', outfile_stab_struc , '\n===================================================================\n') 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 2: write csv # when dynamut, dynamut2 and others weren't being combined #--------------------- print('\nWriting file: combined stability and afor' , '\nOutput 2 filename:', outfile_comb_afor , '\n===================================================================\n') 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, mcsm_ppi2_df] if gene.lower() == "rpob": dfs_list = [dynamut2_df, mcsm_na_df, mcsm_ppi2_df] if gene.lower() == "alr": dfs_list = [dynamut2_df, mcsm_ppi2_df] # noticed that with revised rpoB that mcsm-NA had one less position, # Hence this condition else the last check fails with discrepancy for expected_nrows if len(dfs_list) > 1: join_type = 'outer' else: join_type = 'inner' print('\nUsing join type: "', join_type, '" for the last but one merge') dfs_merged = reduce(lambda left,right: pd.merge(left , right , on = ['mutationinformation'] #, how = 'inner') , how = join_type) , 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) ) # Drop cols if combined_all_params.columns.str.contains(r'_x$|_y$', regex = True).any(): print('\nDuplicate column names detected...' , '\nDropping these before writing file') extra_cols_to_drop = list(combined_all_params.columns.str.extract(r'(.*_x$|.*_y$)', expand = True).dropna()[0]) print('\nTotal cols:', len(combined_all_params.columns) ,'\nDropping:', len(extra_cols_to_drop), 'columns') #extra_cols_to_drop = ['chain_x', 'chain_y'] combined_all_params = combined_all_params.drop(extra_cols_to_drop, axis = 1) else: print('\nNo duplicate column names detected, just writing file' , '\nTotal cols:', len(combined_all_params.columns) ) #%%============================================================================ #--------------------- # Output 3: write csv #--------------------- print('\nWriting file: all params') print('\nOutput 3 filename:', outfile_comb , '\n===================================================================\n') 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