bring in embb stuff which was in the wrong branch
This commit is contained in:
parent
1e3670f935
commit
9c37dbee31
6 changed files with 816 additions and 98 deletions
159
mcsm_ppi2/format_results_mcsm_ppi2.py
Executable file
159
mcsm_ppi2/format_results_mcsm_ppi2.py
Executable file
|
|
@ -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)
|
||||
#%%#####################################################################
|
||||
79
mcsm_ppi2/run_format_results_mcsm_ppi2.py
Executable file
79
mcsm_ppi2/run_format_results_mcsm_ppi2.py
Executable file
|
|
@ -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 + <drug> + input')
|
||||
arg_parser.add_argument('-o', '--output_dir', help = 'Output dir for results. By default, it assmes homedir + <drug> + 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=============================================================')
|
||||
|
||||
#%%#####################################################################
|
||||
|
|
@ -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'
|
||||
|
||||
#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'
|
||||
infile_mcsm = outdir + in_filename_mcsm
|
||||
mcsm_df = pd.read_csv(infile_mcsm, sep = ',')
|
||||
|
||||
infile_mcsm = outdir + in_filename_mcsm
|
||||
infile_foldx = outdir + in_filename_foldx
|
||||
infile_deepddg = outdir + in_filename_deepddg
|
||||
in_filename_foldx = gene.lower() + '_foldx.csv'
|
||||
infile_foldx = outdir + in_filename_foldx
|
||||
foldx_df = pd.read_csv(infile_foldx , sep = ',')
|
||||
|
||||
infile_dssp = outdir + in_filename_dssp
|
||||
infile_kd = outdir + in_filename_kd
|
||||
infile_rd = outdir + in_filename_rd
|
||||
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_snpinfo = outdir + '/' + in_filename_snpinfo
|
||||
infile_afor = outdir + '/' + in_filename_afor
|
||||
#infile_afor_kin = outdir + '/' + in_filename_afor_kin
|
||||
in_filename_dssp = gene.lower() + '_dssp.csv'
|
||||
infile_dssp = outdir + in_filename_dssp
|
||||
dssp_df = pd.read_csv(infile_dssp, 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
|
||||
in_filename_kd = gene.lower() + '_kd.csv'
|
||||
infile_kd = outdir + in_filename_kd
|
||||
kd_df = pd.read_csv(infile_kd, sep = ',')
|
||||
|
||||
#, '\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,23 +500,40 @@ 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 <gene_match>_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))
|
||||
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])
|
||||
#%%============================================================================
|
||||
# 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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
149
scripts/deepddg_format.py
Executable file
149
scripts/deepddg_format.py
Executable file
|
|
@ -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) <gene>.lower()'_mcsm_formatted_snps.csv'
|
||||
#2) <gene>.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 + <drug> + input')
|
||||
arg_parser.add_argument('-o', '--output_dir', help = 'Output dir for results. By default, it assmes homedir + <drug> + 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
|
||||
|
|
@ -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.'
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue