added foldx_scaled and deepddg_scaled values added to combine_df.py and also used that script to merge all the dfs so that merged_df2 and merged_df3 are infact what we need for downstream processing

This commit is contained in:
Tanushree Tunstall 2021-09-10 16:58:36 +01:00
parent dda5d1ea93
commit 4ba4ff602e
5 changed files with 354 additions and 977 deletions

View file

@ -41,6 +41,7 @@ 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('~')
@ -92,19 +93,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
@ -122,49 +110,52 @@ 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
in_filename_dssp = gene.lower() + '_dssp.csv'
in_filename_kd = gene.lower() + '_kd.csv'
in_filename_rd = gene.lower() + '_rd.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
in_filename_dssp = gene.lower() + '_dssp.csv'
in_filename_kd = gene.lower() + '_kd.csv'
in_filename_rd = gene.lower() + '_rd.csv'
#in_filename_snpinfo = 'ns' + gene.lower() + '_snp_info_f.csv' # gwas f info
in_filename_afor = gene.lower() + '_af_or.csv'
in_filename_afor = gene.lower() + '_af_or.csv'
#in_filename_afor_kin = gene.lower() + '_af_or_kinship.csv'
infilename_dynamut = gene.lower() + '_complex_dynamut_norm.csv'
infilename_dynamut2 = gene.lower() + '_complex_dynamut2_norm.csv'
infilename_mcsm_na = gene.lower() + '_complex_mcsm_na_norm.csv'
infilename_mcsm_f_snps = gene.lower() + '_mcsm_formatted_snps.csv'
infile_mcsm = outdir + in_filename_mcsm
infile_foldx = outdir + in_filename_foldx
infile_mcsm = outdir + in_filename_mcsm
infile_foldx = outdir + in_filename_foldx
infile_deepddg = outdir + in_filename_deepddg
infile_dssp = outdir + in_filename_dssp
infile_kd = outdir + in_filename_kd
infile_rd = outdir + in_filename_rd
#infile_snpinfo = outdir + in_filename_snpinfo
infile_afor = outdir + in_filename_afor
#infile_afor_kin = outdir + in_filename_afor_kin
infile_dynamut = outdir + 'dynamut_results/' + infilename_dynamut
infile_dynamut2 = outdir + 'dynamut_results/dynamut2/' + infilename_dynamut2
infile_mcsm_na = outdir + 'mcsm_na_results/' + infilename_mcsm_na
infile_mcsm_f_snps = outdir + infilename_mcsm_f_snps
infile_dssp = outdir + in_filename_dssp
infile_kd = outdir + in_filename_kd
infile_rd = outdir + in_filename_rd
#infile_snpinfo = outdir + '/' + in_filename_snpinfo
infile_afor = outdir + '/' + in_filename_afor
#infile_afor_kin = outdir + '/' + in_filename_afor_kin
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============================================================')
# read csv
mcsm_df = pd.read_csv(infile_mcsm, sep = ',')
foldx_df = pd.read_csv(infile_foldx , sep = ',')
deepddg_df = pd.read_csv(infile_deepddg, sep = ',')
dssp_df = pd.read_csv(infile_dssp, sep = ',')
kd_df = pd.read_csv(infile_kd, sep = ',')
rd_df = pd.read_csv(infile_rd, sep = ',')
afor_df = pd.read_csv(infile_afor, sep = ',')
dynamut_df = pd.read_csv(infile_dynamut, sep = ',')
dynamut2_df = pd.read_csv(infile_dynamut2, sep = ',')
mcsm_na_df = pd.read_csv(infile_mcsm_na, sep = ',')
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===================================================================')
@ -174,12 +165,101 @@ r_join = 'right'
i_join = 'inner'
# end of variable assignment for input and output files
#%%============================================================================
#%%############################################################################
#=====================
# some preprocessing
#=====================
#-------------
# FoldX
#-------------
foldx_df.shape
#=======================
# scale foldx values
#=======================
# 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'].min()
foldx_max = foldx_df['ddg'].max()
foldx_min
foldx_max
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'].apply(foldx_scale)
print('Raw foldx scores:\n', foldx_df['ddg']
, '\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']>=0].count()
foldx_pos = c.get(key = 'ddg')
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======================================================')
# rename ddg column to ddg_foldx
foldx_df['ddg']
foldx_df = foldx_df.rename(columns = {'ddg':'ddg_foldx'})
foldx_df['ddg_foldx']
#-------------
# 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======================================================')
#%%=============================================================================
# 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,7 +269,7 @@ 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)
merging_cols_m1 = detect_common_cols(mcsm_df, foldx_df)
@ -205,8 +285,8 @@ print('==================================='
, '\nSecond merge: mcsm_foldx_dfs + deepddg'
, '\n===================================')
deepddg_df = pd.read_csv(infile_deepddg, sep = ',')
deepddg_df.columns
#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)
@ -218,9 +298,9 @@ 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)
merging_cols_m2 = detect_common_cols(dssp_df, kd_df)
@ -308,8 +388,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 <gene_match>_abcXXXrst
#----------------------
@ -360,16 +440,60 @@ 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
, '\n===================================================================')
# write csv
# # 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])
#%% end of script
#%%============================================================================
# combine dynamut, dynamut2, and mcsm_na
dfs_list = [dynamut_df, dynamut2_df, mcsm_na_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 = i_join)
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