diff --git a/mcsm_ppi2/format_results_mcsm_ppi2.py b/mcsm_ppi2/format_results_mcsm_ppi2.py new file mode 100755 index 0000000..69d4835 --- /dev/null +++ b/mcsm_ppi2/format_results_mcsm_ppi2.py @@ -0,0 +1,159 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Aug 19 14:33:51 2020 + +@author: tanu +""" +#%% load packages +import os,sys +homedir = os.path.expanduser('~') +import subprocess +import argparse +import requests +import re +import time +from bs4 import BeautifulSoup +import pandas as pd +import numpy as np +from pandas.api.types import is_string_dtype +from pandas.api.types import is_numeric_dtype + +sys.path.append(homedir + '/git/LSHTM_analysis/scripts') +from reference_dict import up_3letter_aa_dict +from reference_dict import oneletter_aa_dict + +#%%##################################################################### + +def format_mcsm_ppi2_output(mcsm_ppi2_output_csv): + """ + @param mcsm_ppi2_output_csv: file containing mcsm_ppi2_results for all muts + which is the result of combining all mcsm_ppi2 batch results, and using + bash scripts to combine all the batch results into one file. + Formatting df to a pandas df and output as csv. + @type string + + @return (not true) formatted csv for mcsm_ppi2 output + @type pandas df + + """ + ############# + # Read file + ############# + mcsm_ppi2_data_raw = pd.read_csv(mcsm_ppi2_output_csv, sep = ',') + + # strip white space from both ends in all columns + mcsm_ppi2_data = mcsm_ppi2_data_raw.apply(lambda x: x.str.strip() if x.dtype == 'object' else x) + + dforig_shape = mcsm_ppi2_data.shape + print('dimensions of input file:', dforig_shape) + + ############# + # Map 3 letter + # code to one + ############# + # 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 up_3letter_aa_dict.items(): + lookup_dict[k] = v['one_letter_code'] + wt = mcsm_ppi2_data['wild-type'].squeeze() # converts to a series that map works on + mcsm_ppi2_data['w_type'] = wt.map(lookup_dict) + mut = mcsm_ppi2_data['mutant'].squeeze() + mcsm_ppi2_data['m_type'] = mut.map(lookup_dict) + + # ############# + # # CHECK + # # 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_upper'] + # wt = mcsm_ppi2_data['w_type'].squeeze() #converts to a series that map works on + # mcsm_ppi2_data['WILD'] = wt.map(lookup_dict) + # mut = mcsm_ppi2_data['m_type'].squeeze() + # mcsm_ppi2_data['MUT'] = mut.map(lookup_dict) + + # # check + # mcsm_ppi2_data['wild-type'].equals(mcsm_ppi2_data['WILD']) + # mcsm_ppi2_data['mutant'].equals(mcsm_ppi2_data['MUT']) +#%%============================================================================ + ############# + # rename cols + ############# + # format colnames: all lowercase and consistent colnames + mcsm_ppi2_data.columns + print('Assigning meaningful colnames' + , '\n=======================================================') + + my_colnames_dict = {'chain': 'chain' + , 'wild-type': 'wt_upper' + , 'res-number': 'position' + , 'mutant': 'mut_upper' + , 'distance-to-interface': 'interface_dist' + , 'mcsm-ppi2-prediction': 'mcsm_ppi2_affinity' + , 'affinity': 'mcsm_ppi2_outcome' + , 'w_type': 'wild_type' # one letter amino acid code + , 'm_type': 'mutant_type' # one letter amino acid code +} + + mcsm_ppi2_data.rename(columns = my_colnames_dict, inplace = True) + mcsm_ppi2_data.columns + + ############# + # create mutationinformation column + ############# + #mcsm_ppi2_data['mutationinformation'] = mcsm_ppi2_data['wild_type'] + mcsm_ppi2_data.position.map(str) + mcsm_ppi2_data['mutant_type'] + mcsm_ppi2_data['mutationinformation'] = mcsm_ppi2_data.loc[:,'wild_type'] + mcsm_ppi2_data.loc[:,'position'].astype(int).apply(str) + mcsm_ppi2_data.loc[:,'mutant_type'] + +#%%===================================================================== + ######################### + # scale mcsm_ppi2 values + ######################### + # Rescale values in mcsm_ppi2_affinity col b/w -1 and 1 so negative numbers + # stay neg and pos numbers stay positive + mcsm_ppi2_min = mcsm_ppi2_data['mcsm_ppi2_affinity'].min() + mcsm_ppi2_max = mcsm_ppi2_data['mcsm_ppi2_affinity'].max() + + mcsm_ppi2_scale = lambda x : x/abs(mcsm_ppi2_min) if x < 0 else (x/mcsm_ppi2_max if x >= 0 else 'failed') + + mcsm_ppi2_data['mcsm_ppi2_scaled'] = mcsm_ppi2_data['mcsm_ppi2_affinity'].apply(mcsm_ppi2_scale) + print('Raw mcsm_ppi2 scores:\n', mcsm_ppi2_data['mcsm_ppi2_affinity'] + , '\n---------------------------------------------------------------' + , '\nScaled mcsm_ppi2 scores:\n', mcsm_ppi2_data['mcsm_ppi2_scaled']) + + c = mcsm_ppi2_data[mcsm_ppi2_data['mcsm_ppi2_affinity']>=0].count() + mcsm_ppi2_pos = c.get(key = 'mcsm_ppi2_affinity') + + c2 = mcsm_ppi2_data[mcsm_ppi2_data['mcsm_ppi2_scaled']>=0].count() + mcsm_ppi2_pos2 = c2.get(key = 'mcsm_ppi2_scaled') + + if mcsm_ppi2_pos == mcsm_ppi2_pos2: + print('\nPASS: Affinity values scaled correctly') + else: + print('\nFAIL: Affinity values scaled numbers MISmatch' + , '\nExpected number:', mcsm_ppi2_pos + , '\nGot:', mcsm_ppi2_pos2 + , '\n======================================================') + +#%%===================================================================== + ############# + # reorder columns + ############# + mcsm_ppi2_data.columns + mcsm_ppi2_dataf = mcsm_ppi2_data[['mutationinformation' + , 'mcsm_ppi2_affinity' + , 'mcsm_ppi2_scaled' + , 'mcsm_ppi2_outcome' + , 'interface_dist' + , 'wild_type' + , 'position' + , 'mutant_type' + , 'wt_upper' + , 'mut_upper' + , 'chain']] + return(mcsm_ppi2_dataf) +#%%##################################################################### diff --git a/mcsm_ppi2/run_format_results_mcsm_ppi2.py b/mcsm_ppi2/run_format_results_mcsm_ppi2.py new file mode 100755 index 0000000..65752f5 --- /dev/null +++ b/mcsm_ppi2/run_format_results_mcsm_ppi2.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Feb 12 12:15:26 2021 + +@author: tanu +""" +#%% load packages +import sys, os +homedir = os.path.expanduser('~') +#sys.path.append(homedir + '/git/LSHTM_analysis/mcsm_ppi2') + +from format_results_mcsm_ppi2 import * +######################################################################## +# TODO: add cmd line args +#%% command line args +arg_parser = argparse.ArgumentParser() +arg_parser.add_argument('-d', '--drug' , help = 'drug name (case sensitive)', default = None) +arg_parser.add_argument('-g', '--gene' , help = 'gene name (case sensitive)', default = None) +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('--mkdir_name' , help = 'Output dir for processed results. This will be created if it does not exist') +arg_parser.add_argument('-m', '--make_dirs' , help = 'Make dir for input and output', action='store_true') + +arg_parser.add_argument('--debug' , action = 'store_true' , help = 'Debug Mode') + +args = arg_parser.parse_args() +#%%============================================================================ +# variable assignment: input and output paths & filenames +drug = args.drug +gene = args.gene +datadir = args.datadir +indir = args.input_dir +outdir = args.output_dir +#outdir_ppi2 = args.mkdir_name +make_dirs = args.make_dirs + +#======= +# dirs +#======= +if not datadir: + datadir = homedir + '/git/Data/' + +if not indir: + indir = datadir + drug + '/input/' + +if not outdir: + outdir = datadir + drug + '/output/' + +#if not mkdir_name: +# outdir_ppi2 = outdir + 'mcsm_ppi2/' + +outdir_ppi2 = outdir + 'mcsm_ppi2/' + +# Input file +infile_mcsm_ppi2 = outdir_ppi2 + gene.lower() + '_output_combined_clean.csv' + +# Formatted output file +outfile_mcsm_ppi2_f = outdir_ppi2 + gene.lower() + '_complex_mcsm_ppi2_norm.csv' + +#========================== +# CALL: format_results_mcsm_na() +# Data: gid+streptomycin +#========================== +print('Formatting results for:', infile_mcsm_ppi2) +mcsm_ppi2_df_f = format_mcsm_ppi2_output(mcsm_ppi2_output_csv = infile_mcsm_ppi2) + +# writing file +print('Writing formatted df to csv') +mcsm_ppi2_df_f.to_csv(outfile_mcsm_ppi2_f, index = False) + +print('Finished writing file:' + , '\nFile:', outfile_mcsm_ppi2_f + , '\nExpected no. of rows:', len(mcsm_ppi2_df_f) + , '\nExpected no. of cols:', len(mcsm_ppi2_df_f.columns) + , '\n=============================================================') + +#%%##################################################################### \ No newline at end of file diff --git a/scripts/combining_dfs.py b/scripts/combining_dfs.py index 634af18..53361c7 100755 --- a/scripts/combining_dfs.py +++ b/scripts/combining_dfs.py @@ -34,6 +34,11 @@ Created on Tue Aug 6 12:56:03 2019 # 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 @@ -41,13 +46,14 @@ 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.chdir(homedir + '/git/LSHTM_analysis/scripts') os.getcwd() # FIXME: local imports @@ -92,19 +98,6 @@ outdir = args.output_dir gene_match = gene + '_p.' print('mut pattern for gene', gene, ':', gene_match) -# !"Redundant, now that improvements have been made! -# See section "REGEX" -# nssnp_match = gene_match +'[A-Za-z]{3}[0-9]+[A-Za-z]{3}' -# print('nsSNP for gene', gene, ':', nssnp_match) - -# wt_regex = gene_match.lower()+'([A-Za-z]{3})' -# print('wt regex:', wt_regex) - -# mut_regex = r'[0-9]+(\w{3})$' -# print('mt regex:', mut_regex) - -# pos_regex = r'([0-9]+)' -# print('position regex:', pos_regex) #%%======================================================================= #============== # directories @@ -121,65 +114,263 @@ if not outdir: #======= # input #======= -#in_filename_mcsm = gene.lower() + '_complex_mcsm_norm.csv' -in_filename_mcsm = gene.lower() + '_complex_mcsm_norm_SAM.csv' # gidb -in_filename_foldx = gene.lower() + '_foldx.csv' -in_filename_deepddg = gene.lower() + '_ni_deepddg.csv' # change to decent filename and put it in the correct dir +gene_list_normal = ["pnca", "katg", "rpob", "alr"] -in_filename_dssp = gene.lower() + '_dssp.csv' -in_filename_kd = gene.lower() + '_kd.csv' -in_filename_rd = gene.lower() + '_rd.csv' +if gene.lower() == "gid": + print("\nReading mCSM file for gene:", gene) + in_filename_mcsm = gene.lower() + '_complex_mcsm_norm_SAM.csv' +if gene.lower() == "embb": + print("\nReading mCSM file for gene:", gene) + in_filename_mcsm = gene.lower() + '_complex_mcsm_norm1.csv' +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_snpinfo = 'ns' + gene.lower() + '_snp_info_f.csv' # gwas f info -in_filename_afor = gene.lower() + '_af_or.csv' -#in_filename_afor_kin = gene.lower() + '_af_or_kinship.csv' +in_filename_foldx = gene.lower() + '_foldx.csv' +infile_foldx = outdir + in_filename_foldx +foldx_df = pd.read_csv(infile_foldx , sep = ',') -infile_mcsm = outdir + in_filename_mcsm -infile_foldx = outdir + in_filename_foldx -infile_deepddg = outdir + in_filename_deepddg +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 = ',') -infile_dssp = outdir + in_filename_dssp -infile_kd = outdir + in_filename_kd -infile_rd = outdir + in_filename_rd +in_filename_dssp = gene.lower() + '_dssp.csv' +infile_dssp = outdir + in_filename_dssp +dssp_df = pd.read_csv(infile_dssp, sep = ',') -#infile_snpinfo = outdir + '/' + in_filename_snpinfo -infile_afor = outdir + '/' + in_filename_afor -#infile_afor_kin = outdir + '/' + in_filename_afor_kin +in_filename_kd = gene.lower() + '_kd.csv' +infile_kd = outdir + in_filename_kd +kd_df = pd.read_csv(infile_kd, sep = ',') -print('\nInput path:', indir - , '\nOutput path:', outdir, '\n' - , '\nInput filename mcsm:', infile_mcsm - , '\nInput filename foldx:', infile_foldx, '\n' - , '\nInput filename deepddg', infile_deepddg , '\n' - , '\nInput filename dssp:', infile_dssp - , '\nInput filename kd:', infile_kd - , '\nInput filename rd', infile_rd - - #, '\nInput filename snp info:', infile_snpinfo, '\n' - , '\nInput filename af or:', infile_afor - #, '\nInput filename afor kinship:', infile_afor_kin - , '\n============================================================') +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 = ',') + +#------------------------------------------------------------ +# ONLY:for gene pnca and gid: End logic should pick this up! +geneL_dy_na = ["pnca", "gid"] +#if gene.lower() == "pnca" or "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 "embb" or "alr": + 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 = ',') +#-------------------------------------------------------------- +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) #======= # output #======= out_filename_comb = gene.lower() + '_all_params.csv' -outfile_comb = outdir + '/' + out_filename_comb +outfile_comb = outdir + out_filename_comb print('Output filename:', outfile_comb , '\n===================================================================') - -o_join = 'outer' -l_join = 'left' -r_join = 'right' -i_join = 'inner' - # 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('Raw 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 + +#------------------------- +# 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('Raw 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() + +#-------------------------- +# 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) + +#-------------------------- +# subset chain +#-------------------------- +if gene.lower() == "embb": + sel_chain = "B" +else: + sel_chain = "A" + +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 itCheck for duplicates +#-------------------------- +col_to_drop = ['chain_id'] +deepddg_df = deepddg_df.drop(col_to_drop, axis = 1) + +#%%============================================================================= +# Now merges begin +#%%============================================================================= print('===================================' , '\nFirst merge: mcsm + foldx' , '\n===================================') -mcsm_df = pd.read_csv(infile_mcsm, sep = ',') +mcsm_df.shape # add 3 lowercase aa code for wt and mutant get_aa_3lower(df = mcsm_df @@ -189,42 +380,100 @@ get_aa_3lower(df = mcsm_df , col_mut = 'mut_aa_3lower') #mcsm_df.columns = mcsm_df.columns.str.lower() -foldx_df = pd.read_csv(infile_foldx , sep = ',') +# foldx_df.shape -#mcsm_foldx_dfs = combine_dfs_with_checks(mcsm_df, foldx_df, my_join = o_join) +#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 = o_join) +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 +# get the empty cells to be full of meaningful info +if mcsm_foldx_dfs.loc[:,'wild_type': 'mut_aa_3lower'].isnull().values.any(): + print ("NAs detected in mcsm cols after merge") + + ############################## + # 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) + #%% print('===================================' , '\nSecond merge: mcsm_foldx_dfs + deepddg' , '\n===================================') -deepddg_df = pd.read_csv(infile_deepddg, sep = ',') -deepddg_df.columns - # merge with mcsm_foldx_dfs and deepddg_df -mcsm_foldx_deepddg_dfs = pd.merge(mcsm_foldx_dfs, deepddg_df, on = 'mutationinformation', how = l_join) +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') + #%%============================================================================ print('===================================' , '\Third merge: dssp + kd' , '\n===================================') -dssp_df = pd.read_csv(infile_dssp, sep = ',') -kd_df = pd.read_csv(infile_kd, sep = ',') -rd_df = pd.read_csv(infile_rd, sep = ',') +dssp_df.shape +kd_df.shape +rd_df.shape -#dssp_kd_dfs = combine_dfs_with_checks(dssp_df, kd_df, my_join = o_join) +#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 = o_join) +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===================================================================') @@ -233,10 +482,12 @@ 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 = o_join) +#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 = o_join) +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) @@ -249,24 +500,41 @@ 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 = i_join) + +#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 = i_join) +#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 = i_join) +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) -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') - +# 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===================================================================') @@ -281,7 +549,7 @@ 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 +# sanity check foo = combined_df[['wild_type', 'wild_type_kd', 'wt_3letter_caps', 'wt_aa_3lower', 'mut_aa_3lower']] # Drop cols @@ -308,8 +576,8 @@ print('\n=======================================' , '\ncombined_df_clean + afor_df ' , '\n=======================================') -afor_df = pd.read_csv(infile_afor, sep = ',') afor_cols = afor_df.columns +afor_df.shape # create a mapping from the gwas mutation column i.e _abcXXXrst #---------------------- @@ -335,7 +603,11 @@ 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 = l_join) +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) @@ -347,20 +619,28 @@ if len(combined_stab_afor) == len(combined_df_clean) and len(combined_stab_afor. else: sys.exit('\nFAIL: check individual df merges') -print('\n\nResult of Fourth merge:', combined_stab_afor.shape +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' - , '\nNo. of nsSNPs with valid ORs: ', len(afor_df)) -else: - sys.exit('\nFAIL: merge unsuccessful for af and or') +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 +# 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 @@ -372,4 +652,61 @@ 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]) -#%% end of script +#%%============================================================================ +# combine dynamut, dynamut2, and mcsm_na +#dfs_list = [dynamut_df, dynamut2_df, mcsm_na_df] # gid + +if gene.lower() == "pnca": + dfs_list = [dynamut_df, 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 \ No newline at end of file diff --git a/scripts/data_extraction.py b/scripts/data_extraction.py index 5582632..aac7cdb 100755 --- a/scripts/data_extraction.py +++ b/scripts/data_extraction.py @@ -70,7 +70,6 @@ arg_parser.add_argument('-m', '--make_dirs', help = 'Make dir for input and outp arg_parser.add_argument('--debug', action ='store_true', help = 'Debug Mode') - args = arg_parser.parse_args() #======================================================================= #%% variable assignment: input and output paths & filenames @@ -81,9 +80,6 @@ indir = args.input_dir outdir = args.output_dir make_dirs = args.make_dirs -#drug = 'streptomycin' -#gene = 'gid' - #%% input and output dirs and files #======= # dirs @@ -1373,4 +1369,4 @@ if dr_muts.isin(other_muts).sum() & other_muts.isin(dr_muts).sum() > 0: print(u'\u2698' * 50, '\nEnd of script: Data extraction and writing files' '\n' + u'\u2698' * 50 ) -#%% end of script \ No newline at end of file +#%% end of script diff --git a/scripts/deepddg_format.py b/scripts/deepddg_format.py new file mode 100755 index 0000000..20b2dcb --- /dev/null +++ b/scripts/deepddg_format.py @@ -0,0 +1,149 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +''' +Created on Tue Aug 6 12:56:03 2019 + +@author: tanu +''' +#======================================================================= +# Task: format deep ddg df to allow easy merging + +# Input: 2 dfs +#1) .lower()'_mcsm_formatted_snps.csv' +#2) .lower()_complex_ddg_results.csv' +#======================================================================= +#%% load packages +import sys, os +import pandas as pd +from pandas import DataFrame +import numpy as np +#from varname import nameof +import argparse +#======================================================================= +#%% specify input and curr dir +homedir = os.path.expanduser('~') + +# set working dir +os.getcwd() +os.chdir(homedir + '/git/LSHTM_analysis/scripts') +os.getcwd() +#=======================================================================#%% 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 +#%%======================================================================= +#============== +# directories +#============== +if not datadir: + datadir = homedir + '/git/Data/' + +if not indir: + indir = datadir + drug + '/input/' + +if not outdir: + outdir = datadir + drug + '/output/' + +#======= +# input +#======= +in_filename_mcsm_snps = gene.lower() + '_mcsm_formatted_snps.csv' +infile_mcsm_snps = outdir + in_filename_mcsm_snps + +in_filename_deepddg = gene.lower() + '_complex_ddg_results.csv' +infile_deepddg = outdir + 'deep_ddg/' + in_filename_deepddg + +print('\nInput path:', indir + , '\nOutput path:', outdir, '\n' + , '\nInput filename mcsm snps', infile_mcsm_snps , '\n' + , '\nInput filename deepddg', infile_deepddg , '\n' + , '\n============================================================') + +#======= +# output +#======= +#out_filename_deepddg = gene.lower() + '_ni_deepddg.txt' +out_filename_deepddg = gene.lower() + '_ni_deepddg.csv' +outfile_deepddg_f = outdir + out_filename_deepddg + +print('Output filename:', outfile_deepddg_f + , '\n===================================================================') +# end of variable assignment for input and output files +#%%============================================================================ +print('===================================' + , '\nmcsm muts' + , '\n===================================') + +mcsm_muts_df = pd.read_csv(infile_mcsm_snps , header = None, sep = ',', names = ['mutationinformation']) +mcsm_muts_df.columns + +#%%============================================================================ +print('===================================' + , '\nDeep ddg' + , '\n===================================') + +deepddg_df = pd.read_csv(infile_deepddg, sep = ',') +deepddg_df.columns + +deepddg_df.rename(columns = {'#chain' : 'chain_id' + , 'WT' : 'wild_type_deepddg' + , 'ResID' : 'position' + , 'Mut' : 'mutant_type_deepddg'} + , inplace = True) +deepddg_df.columns +deepddg_df['mutationinformation'] = deepddg_df['wild_type_deepddg'] + deepddg_df['position'].map(str) + deepddg_df['mutant_type_deepddg'] +deepddg_df.columns + +# add deepddg outcome column: <0--> Destabilising, >0 --> Stabilising +deepddg_df['deepddg_outcome'] = np.where(deepddg_df['deepddg'] < 0, 'Destabilising', 'Stabilising') +deepddg_df['deepddg_outcome'].value_counts() + +# should be identical in count ot Destabilising and stabilising respectively +len(deepddg_df.loc[deepddg_df['deepddg'] < 0]) +len(deepddg_df.loc[deepddg_df['deepddg'] >= 0]) + +#---------------------------------------------- +# drop extra columns to allow clean merging +#---------------------------------------------- +#deepddg_short_df = deepddg_df.drop(['chain_id', 'wild_type_deepddg', 'position', 'mutant_type_deepddg'], axis = 1) + +#---------------------------------------------- +# embb (where gene-target has > 1 chain) +# include chain else the numbering will be messed up! +#---------------------------------------------- +deepddg_short_df = deepddg_df.drop(['wild_type_deepddg', 'position', 'mutant_type_deepddg'], axis = 1) + +# rearrange columns +deepddg_short_df.columns +deepddg_short_df = deepddg_short_df[["chain_id", "mutationinformation", "deepddg", "deepddg_outcome"]] + +#%% combine with mcsm snps +deepddg_mcsm_muts_dfs = pd.merge(deepddg_short_df + , mcsm_muts_df + , on = 'mutationinformation' + , how = 'right') +deepddg_mcsm_muts_dfs ['deepddg_outcome'].value_counts() + +#%%============================================================================ +# write csv +print('Writing file: formatted deepddg and only mcsm muts') +deepddg_mcsm_muts_dfs.to_csv(outfile_deepddg_f, index = False) +print('\nFinished writing file:' + , '\nNo. of rows:', deepddg_mcsm_muts_dfs.shape[0] + , '\nNo. of cols:', deepddg_mcsm_muts_dfs.shape[1]) +#%% end of script diff --git a/scripts/rd_df.py b/scripts/rd_df.py index 7eab903..102530d 100755 --- a/scripts/rd_df.py +++ b/scripts/rd_df.py @@ -45,8 +45,6 @@ arg_parser.add_argument('--debug', action='store_true', help = 'Debug Mode') args = arg_parser.parse_args() #======================================================================= #%% variable assignment: input and output -#drug = 'pyrazinamide' -#gene = 'pncA' drug = args.drug gene = args.gene gene_match = gene + '_p.'