tanu/LSHTM_analysis
1
0
Fork 0
Code Issues Pull requests Releases Wiki Activity

defined method for formatting mcsm_results

Browse source
This commit is contained in:
Tanushree Tunstall 2020-04-14 11:30:36 +01:00
parent 45889990e7
commit 7aafa72e10
2 changed files with 245 additions and 206 deletions
Download patch file Download diff file Expand all files Collapse all files

413
mcsm/format_results.py
View file

@ -6,9 +6,9 @@
import os,sys import os,sys
import subprocess import subprocess
import argparse import argparse
import requests #import requests
import re import re
import time #import time
import pandas as pd import pandas as pd
from pandas.api.types import is_string_dtype from pandas.api.types import is_string_dtype
from pandas.api.types import is_numeric_dtype from pandas.api.types import is_numeric_dtype
@ -43,7 +43,7 @@ datadir = homedir + '/' + 'git/Data'
#======= #=======
# 1) result_urls (from outdir) # 1) result_urls (from outdir)
outdir = datadir + '/' + drug + '/' + 'output' outdir = datadir + '/' + drug + '/' + 'output'
in_filename = gene.lower() + '_mcsm_output.csv' #(outfile, from mcsm_results) in_filename = gene.lower() + '_mcsm_output.csv' #(outfile, from mcsm_results.py)
infile = outdir + '/' + in_filename infile = outdir + '/' + in_filename
print('Input filename:', in_filename print('Input filename:', in_filename
, '\nInput path(from output dir):', outdir , '\nInput path(from output dir):', outdir
@ -53,215 +53,256 @@ print('Input filename:', in_filename
# output # output
#======= #=======
outdir = datadir + '/' + drug + '/' + 'output' outdir = datadir + '/' + drug + '/' + 'output'
out_filename = gene.lower() + '_complex_mcsm_norm.csv' out_filename = gene.lower() + '_complex_mcsm_results.csv'
outfile = outdir + '/' + out_filename outfile = outdir + '/' + out_filename
print('Output filename:', out_filename print('Output filename:', out_filename
, '\nOutput path:', outdir , '\nOutput path:', outdir
, '\n=============================================================') , '\n=============================================================')
#======================================================================= #=======================================================================
print('Reading input file') def format_mcsm_output(mcsm_outputcsv):
mcsm_data = pd.read_csv(infile, sep = ',') """
@param mcsm_outputcsv: file containing mcsm results for all muts
(outfile from from mcsm_results.py)
@type string
mcsm_data.columns @return formatted mcsm output
# PredAffLog = affinity_change_log @type pandas df
# "DUETStability_Kcalpermol = DUET_change_kcalpermol
dforig_shape = mcsm_data.shape
print('dim of infile:', dforig_shape)
# change colnames to reflect units and no spaces, and replace '-' with '-' """
print('Assigning meaningful colnames i.e without spaces and hyphen and reflecting units' #############
, '\n===================================================================') # Read file
my_colnames_dict = {'Predicted Affinity Change': 'PredAffLog' #############
, 'Mutation information': 'Mutationinformation' mcsm_data = pd.read_csv(infile, sep = ',')
, 'Wild-type': 'Wild_type' dforig_shape = mcsm_data.shape
, 'Position': 'Position' print('dim of infile:', dforig_shape)
, 'Mutant-type': 'Mutant_type'
, 'Chain': 'Chain'
, 'Ligand ID': 'LigandID'
, 'Distance to ligand': 'Dis_lig_Ang'
, 'DUET stability change': 'DUET_change_kcalpermol'}
mcsm_data.rename(columns = my_colnames_dict, inplace = True) #############
mcsm_data.columns # rename cols
#%%=========================================================================== #############
# populate mutationinformation column:mcsm style muts {WT}<POS>{MUT} # format colnames: all lowercase, remove spaces and use '_' to join
print('Populating column : Mutationinformation which is currently empty\n', mcsm_data['Mutationinformation']) print('Assigning meaningful colnames i.e without spaces and hyphen and reflecting units'
mcsm_data['Mutationinformation'] = mcsm_data['Wild_type'] + mcsm_data['Position'].astype(str) + mcsm_data['Mutant_type'] , '\n===================================================================')
print('checking after populating:\n', mcsm_data['Mutationinformation'] my_colnames_dict = {'Predicted Affinity Change': 'PredAffLog' # relevant info from this col will be extracted and the column discarded
, '\n===================================================================') , 'Mutation information': 'mutation_information' # {wild_type}<position>{mutant_type}
, 'Wild-type': 'wild_type' # one letter amino acid code
, 'Position': 'position' # number
, 'Mutant-type': 'mutant_type' # one letter amino acid code
, 'Chain': 'chain' # single letter (caps)
, 'Ligand ID': 'ligand_id' # 3-letter code
, 'Distance to ligand': 'ligand_distance' # angstroms
, 'DUET stability change': 'duet_stability_change'} # in kcal/mol
# Remove spaces b/w pasted columns mcsm_data.rename(columns = my_colnames_dict, inplace = True)
print('removing white space within column: \Mutationinformation') #%%===========================================================================
mcsm_data['Mutationinformation'] = mcsm_data['Mutationinformation'].str.replace(' ', '') #################################
print('Correctly formatted column: Mutationinformation\n', mcsm_data['Mutationinformation'] # populate mutation_information
, '\n===================================================================') # col which is currently blank
#%%=========================================================================== #################################
# very important # populate mutation_information column:mcsm style muts {WT}<POS>{MUT}
print('Sanity check:' print('Populating column : mutation_information which is currently empty\n', mcsm_data['mutation_information'])
, '\nChecking duplicate mutations') mcsm_data['mutation_information'] = mcsm_data['wild_type'] + mcsm_data['position'].astype(str) + mcsm_data['mutant_type']
if mcsm_data['Mutationinformation'].duplicated().sum() == 0: print('checking after populating:\n', mcsm_data['mutation_information']
print('PASS: No duplicate mutations detected (as expected)' , '\n===================================================================')
, '\nDim of data:', mcsm_data.shape
, '\n===============================================================')
else:
print('FAIL (but not fatal): Duplicate mutations detected'
, '\nDim of df with duplicates:', mcsm_data.shape
, 'Removing duplicate entries')
mcsm_data = mcsm_data.drop_duplicates(['Mutationinformation'])
print('Dim of data after removing duplicate muts:', mcsm_data.shape
, '\n===============================================================')
#%%===========================================================================
# create DUET_outcome column: classification based on DUET stability values
print('Assigning col: DUET_outcome based on DUET stability values')
print('Sanity check:')
# count positive values in the DUET column
c = mcsm_data[mcsm_data['DUET_change_kcalpermol']>=0].count()
DUET_pos = c.get(key = 'DUET_change_kcalpermol')
# Assign category based on sign (+ve : Stabilising, -ve: Destabilising, Mind the spelling (British spelling))
mcsm_data['DUET_outcome'] = np.where(mcsm_data['DUET_change_kcalpermol']>=0, 'Stabilising', 'Destabilising')
mcsm_data['DUET_outcome'].value_counts()
if DUET_pos == mcsm_data['DUET_outcome'].value_counts()['Stabilising']:
print('PASS: DUET outcome assigned correctly')
else:
print('FAIL: DUET outcome assigned incorrectly'
, '\nExpected no. of stabilising mutations:', DUET_pos
, '\nGot no. of stabilising mutations', mcsm_data['DUET_outcome'].value_counts()['Stabilising']
, '\n===============================================================')
#%%===========================================================================
# Extract only the numeric part from col: Dis_lig_Ang
# number: '-?\d+\.?\d*'
mcsm_data['Dis_lig_Ang']
print('extracting numeric part of col: Dis_lig_Ang')
mcsm_data['Dis_lig_Ang'] = mcsm_data['Dis_lig_Ang'].str.extract('(\d+\.?\d*)')
mcsm_data['Dis_lig_Ang']
# changing dtype to numeric # Remove spaces b/w pasted columns
#if is_numeric_dtype(mcsm_data['Dis_lig_Ang']): print('removing white space within column: \mutation_information')
# print('Data type is already numeric, doing nothing') mcsm_data['mutation_information'] = mcsm_data['mutation_information'].str.replace(' ', '')
#else: print('Correctly formatted column: mutation_information\n', mcsm_data['mutation_information']
# print('Changing dtype in col: Dis_lig_Ang to numeric since Distance should be numeric') , '\n===================================================================')
## FIXME: either do it here, or in the end for all the required cols at once #%%===========================================================================
#%%=========================================================================== #############
# create Lig_outcome column: classification based on affinity change values # sanity check: drop dupliate muts
# the numerical and categorical parts need to be extracted from column: PredAffLog #############
# regex used # shouldn't exist as this should be eliminated at the time of running mcsm
# number: '-?\d+\.?\d*' print('Sanity check:'
# category: '\b(\w+ing)\b' , '\nChecking duplicate mutations')
print('Extracting numerical and categorical parts from the col: PredAffLog') if mcsm_data['mutation_information'].duplicated().sum() == 0:
print('to create two columns: affinity_change_log and Lig_outcome' print('PASS: No duplicate mutations detected (as expected)'
, '\n===================================================================') , '\nDim of data:', mcsm_data.shape
, '\n===============================================================')
else:
print('FAIL (but not fatal): Duplicate mutations detected'
, '\nDim of df with duplicates:', mcsm_data.shape
, 'Removing duplicate entries')
mcsm_data = mcsm_data.drop_duplicates(['mutation_information'])
print('Dim of data after removing duplicate muts:', mcsm_data.shape
, '\n===============================================================')
#%%===========================================================================
#############
# Create col: duet_outcome
#############
# classification based on DUET stability values
print('Assigning col: duet_outcome based on DUET stability values')
print('Sanity check:')
# count positive values in the DUET column
c = mcsm_data[mcsm_data['duet_stability_change']>=0].count()
DUET_pos = c.get(key = 'duet_stability_change')
# Assign category based on sign (+ve : Stabilising, -ve: Destabilising, Mind the spelling (British spelling))
mcsm_data['duet_outcome'] = np.where(mcsm_data['duet_stability_change']>=0, 'Stabilising', 'Destabilising')
mcsm_data['duet_outcome'].value_counts()
if DUET_pos == mcsm_data['duet_outcome'].value_counts()['Stabilising']:
print('PASS: DUET outcome assigned correctly')
else:
print('FAIL: DUET outcome assigned incorrectly'
, '\nExpected no. of stabilising mutations:', DUET_pos
, '\nGot no. of stabilising mutations', mcsm_data['duet_outcome'].value_counts()['Stabilising']
, '\n===============================================================')
#%%===========================================================================
#############
# Extract numeric
# part of ligand_distance col
#############
# Extract only the numeric part from col: ligand_distance
# number: '-?\d+\.?\d*'
mcsm_data['ligand_distance']
print('extracting numeric part of col: ligand_distance')
mcsm_data['ligand_distance'] = mcsm_data['ligand_distance'].str.extract('(\d+\.?\d*)')
mcsm_data['ligand_distance']
#%%===========================================================================
#############
# Create 2 columns:
# ligand_affinity_change and ligand_outcome
#############
# the numerical and categorical parts need to be extracted from column: PredAffLog
# regex used
# numerical part: '-?\d+\.?\d*'
# categorocal part: '\b(\w+ing)\b'
print('Extracting numerical and categorical parts from the col: PredAffLog')
print('to create two columns: ligand_affinity_change and ligand_outcome'
, '\n===================================================================')
# Extracting the predicted affinity change (numerical part) # 1) Extracting the predicted affinity change (numerical part)
mcsm_data['affinity_change_log'] = mcsm_data['PredAffLog'].str.extract('(-?\d+\.?\d*)', expand = True) mcsm_data['ligand_affinity_change'] = mcsm_data['PredAffLog'].str.extract('(-?\d+\.?\d*)', expand = True)
print(mcsm_data['affinity_change_log']) print(mcsm_data['ligand_affinity_change'])
# Extracting the categorical part (Destabillizing and Stabilizing) using word boundary ('ing')
#aff_regex = re.compile(r'\b(\w+ing)\b')
mcsm_data['Lig_outcome']= mcsm_data['PredAffLog'].str.extract(r'(\b\w+ing\b)', expand = True)
print(mcsm_data['Lig_outcome'])
print(mcsm_data['Lig_outcome'].value_counts())
american_spl = mcsm_data['Lig_outcome'].value_counts()
print('Changing to Bristish spellings for col: Lig_outcome')
mcsm_data['Lig_outcome'].replace({'Destabilizing': 'Destabilising', 'Stabilizing': 'Stabilising'}, inplace = True)
print(mcsm_data['Lig_outcome'].value_counts())
british_spl = mcsm_data['Lig_outcome'].value_counts()
# since series object will have different names on account of our spelling change # 2) Extracting the categorical part (Destabillizing and Stabilizing) using word boundary ('ing')
# use .equals #aff_regex = re.compile(r'\b(\w+ing)\b')
check = american_spl.values == british_spl.values mcsm_data['ligand_outcome']= mcsm_data['PredAffLog'].str.extract(r'(\b\w+ing\b)', expand = True)
if check.all(): print(mcsm_data['ligand_outcome'])
print('PASS: spelling change successfull' print(mcsm_data['ligand_outcome'].value_counts())
, '\nNo. of predicted affinity changes:\n', british_spl
, '\n===============================================================')
else:
print('FAIL: spelling change unsucessfull'
, '\nExpected:\n', american_spl
, '\nGot:\n', british_spl
, '\n===============================================================')
#%%===========================================================================
# check dtype in cols
print('Checking dtypes in all columns:\n', mcsm_data.dtypes
, '\n===================================================================')
print('Converting the following cols to numeric:'
, '\nDis_lig_Ang'
, '\nDUET_change_kcalpermol'
, '\naffinity_change_log'
, '\n===================================================================')
# using apply method to change stabilty and affinity values to numeric
numeric_cols = ['DUET_change_kcalpermol', 'affinity_change_log', 'Dis_lig_Ang']
mcsm_data[numeric_cols] = mcsm_data[numeric_cols].apply(pd.to_numeric)
# check dtype in cols
print('checking dtype after conversion')
cols_check = mcsm_data.select_dtypes(include='float64').columns.isin(numeric_cols)
if cols_check.all():
print('PASS: dtypes for selected cols:', numeric_cols
, '\nchanged to numeric'
, '\n===============================================================')
else:
print('FAIL:dtype change to numeric for selected cols unsuccessful'
, '\n===============================================================')
#mcsm_data['Dis_lig_Ang', 'affinity_change_log'].apply(is_numeric_dtype(mcsm_data['Dis_lig_Ang', 'affinity_change_log']))
print(mcsm_data.dtypes)
#%%=========================================================================== #############
# changing spelling: British
#############
# ensuring spellings are consistent
american_spl = mcsm_data['ligand_outcome'].value_counts()
print('Changing to Bristish spellings for col: ligand_outcome')
mcsm_data['ligand_outcome'].replace({'Destabilizing': 'Destabilising', 'Stabilizing': 'Stabilising'}, inplace = True)
print(mcsm_data['ligand_outcome'].value_counts())
british_spl = mcsm_data['ligand_outcome'].value_counts()
# compare series values since index will differ from spelling change
check = american_spl.values == british_spl.values
if check.all():
print('PASS: spelling change successfull'
, '\nNo. of predicted affinity changes:\n', british_spl
, '\n===============================================================')
else:
print('FAIL: spelling change unsucessfull'
, '\nExpected:\n', american_spl
, '\nGot:\n', british_spl
, '\n===============================================================')
#%%===========================================================================
#############
# ensuring corrrect dtype columns
#############
# check dtype in cols
print('Checking dtypes in all columns:\n', mcsm_data.dtypes
, '\n===================================================================')
print('Converting the following cols to numeric:'
, '\nligand_distance'
, '\nduet_stability_change'
, '\nligand_affinity_change'
, '\n===================================================================')
#%%=========================================================================== # using apply method to change stabilty and affinity values to numeric
# Normalise the DUET and affinity change cols numeric_cols = ['duet_stability_change', 'ligand_affinity_change', 'ligand_distance']
#converter = lambda x : x*2 if x < 10 else (x*3 if x < 20 else x) mcsm_data[numeric_cols] = mcsm_data[numeric_cols].apply(pd.to_numeric)
duet_min = mcsm_data['DUET_change_kcalpermol'].min() # check dtype in cols
duet_max = mcsm_data['DUET_change_kcalpermol'].max() print('checking dtype after conversion')
cols_check = mcsm_data.select_dtypes(include='float64').columns.isin(numeric_cols)
if cols_check.all():
print('PASS: dtypes for selected cols:', numeric_cols
, '\nchanged to numeric'
, '\n===============================================================')
else:
print('FAIL:dtype change to numeric for selected cols unsuccessful'
, '\n===============================================================')
print(mcsm_data.dtypes)
#%%===========================================================================
converter = lambda x : x/abs(duet_min) if x < 0 else (x/duet_max if x >= 0 else 'failed') #############
# scale duet values
#############
# Rescale values in DUET_change col b/w -1 and 1 so negative numbers
# stay neg and pos numbers stay positive
duet_min = mcsm_data['duet_stability_change'].min()
duet_max = mcsm_data['duet_stability_change'].max()
mcsm_data['DUET_change_kcalpermol'] duet_scale = lambda x : x/abs(duet_min) if x < 0 else (x/duet_max if x >= 0 else 'failed')
mcsm_data['ratioDUET'] = mcsm_data['DUET_change_kcalpermol'].apply(converter)
mcsm_data['ratioDUET']
#%%===========================================================================
# Normalise the affinity change cols
aff_min = mcsm_data['affinity_change_log'].min()
aff_max = mcsm_data['affinity_change_log'].max()
converter = lambda x : x/abs(aff_min) if x < 0 else (x/aff_max if x >= 0 else 'failed') mcsm_data['duet_scaled'] = mcsm_data['duet_stability_change'].apply(duet_scale)
#converter(mcsm_data['affinity_change_log']) print('Raw duet scores:\n', mcsm_data['duet_stability_change']
, '\n---------------------------------------------------------------'
, '\nScaled duet scores:\n', mcsm_data['duet_scaled'])
mcsm_data['affinity_change_log'] #%%===========================================================================
mcsm_data['ratioPredAff'] = mcsm_data['affinity_change_log'].apply(converter) #############
mcsm_data['ratioPredAff'] # scale affinity values
#============================================================================= #############
# Removing PredAff log column as it is not needed? # rescale values in affinity change col b/w -1 and 1 so negative numbers
print('Removing col: PredAffLog since relevant info has been extracted from it') # stay neg and pos numbers stay positive
mcsm_dataf = mcsm_data.drop(columns = ['PredAffLog']) aff_min = mcsm_data['ligand_affinity_change'].min()
#%%=========================================================================== aff_max = mcsm_data['ligand_affinity_change'].max()
expected_cols_toadd = 4
dforig_len = dforig_shape[1] aff_scale = lambda x : x/abs(aff_min) if x < 0 else (x/aff_max if x >= 0 else 'failed')
expected_cols = dforig_len + expected_cols_toadd
if len(mcsm_dataf.columns) == expected_cols: mcsm_data['affinity_scaled'] = mcsm_data['ligand_affinity_change'].apply(aff_scale)
print('PASS: formatting successful' print('Raw affinity scores:\n', mcsm_data['ligand_affinity_change']
, '\nformatted df has expected no. of cols:', expected_cols , '\n---------------------------------------------------------------'
, '\n---------------------------------------------------------------' , '\nScaled affinity scores:\n', mcsm_data['affinity_scaled'])
, '\ncolnames:', mcsm_dataf.columns #=============================================================================
, '\n----------------------------------------------------------------' # Removing PredAff log column as it is not needed?
, '\ndtypes in cols:', mcsm_dataf.dtypes print('Removing col: PredAffLog since relevant info has been extracted from it')
, '\n----------------------------------------------------------------' mcsm_dataf = mcsm_data.drop(columns = ['PredAffLog'])
, '\norig data shape:', dforig_shape #%%===========================================================================
, '\nformatted df shape:', mcsm_dataf.shape #############
, '\n===============================================================') # sanity check before writing file
else: #############
print('FAIL: something went wrong in formatting df' expected_cols_toadd = 4
, '\nExpected no. of cols:', expected_cols dforig_len = dforig_shape[1]
, '\nGot no. of cols:', len(mcsm_dataf.columns) expected_cols = dforig_len + expected_cols_toadd
, '\nCheck formatting' if len(mcsm_dataf.columns) == expected_cols:
, '\n===============================================================') print('PASS: formatting successful'
, '\nformatted df has expected no. of cols:', expected_cols
, '\n---------------------------------------------------------------'
, '\ncolnames:', mcsm_dataf.columns
, '\n----------------------------------------------------------------'
, '\ndtypes in cols:', mcsm_dataf.dtypes
, '\n----------------------------------------------------------------'
, '\norig data shape:', dforig_shape
, '\nformatted df shape:', mcsm_dataf.shape
, '\n===============================================================')
else:
print('FAIL: something went wrong in formatting df'
, '\nExpected no. of cols:', expected_cols
, '\nGot no. of cols:', len(mcsm_dataf.columns)
, '\nCheck formatting'
, '\n===============================================================')
return mcsm_dataf
#%%============================================================================ #%%============================================================================
# call function
mcsm_df_formatted = format_mcsm_output(infile)
# writing file # writing file
print('Writing formatted df to csv') print('Writing formatted df to csv')
mcsm_dataf.to_csv(outfile, index = False) mcsm_df_formatted.to_csv(outfile, index = False)
print('Finished writing file:' print('Finished writing file:'
, '\nFilename:', out_filename , '\nFilename:', out_filename
, '\nPath:', outdir , '\nPath:', outdir
, '\nExpected no. of rows:', len(mcsm_dataf) , '\nExpected no. of rows:', len(mcsm_df_formatted)
, '\nExpected no. of cols:', len(mcsm_dataf.columns) , '\nExpected no. of cols:', len(mcsm_df_formatted)
, '\n=============================================================') , '\n=============================================================')
#%% #%%
#End of script #End of script

26
mcsm/mcsm_results.py
View file

@ -58,13 +58,10 @@ print('Output filename:', out_filename
, '\nOutput path:', outdir , '\nOutput path:', outdir
, '\n=============================================================') , '\n=============================================================')
#%% global variables
#HOST = "http://biosig.unimelb.edu.au"
#PREDICTION_URL = f"{HOST}/mcsm_lig/prediction"
#======================================================================= #=======================================================================
def fetch_results(urltextfile): def fetch_results(urltextfile):
""" """
Extract results data from the results page Extract results data using the prediction url
@params result_page of request_results() @params result_page of request_results()
@type response object @type response object
@ -90,36 +87,37 @@ def fetch_results(urltextfile):
def build_result_dict(web_result_raw): def build_result_dict(web_result_raw):
""" """
Format web results which is inconveniently preformatted! Build dict of mcsm output for a single mutation
Format web results which is preformatted to enable building result dict
# preformatted string object: Problematic! # preformatted string object: Problematic!
# make format consistent
# roughly bring these to the same format as the other @params web_result_raw: directly from html parser extraction
@params web_result_raw directly from html parser extraction
@type string @type string
@returns result dict @returns result dict
@type {} @type {}
""" """
# remove blank lines from output # remove blank lines from output
mytext = os.linesep.join([s for s in web_result_raw.splitlines() if s]) mytext = os.linesep.join([s for s in web_result_raw.splitlines() if s])
# Predicted affintiy change and DUET stability change cols # affinity change and DUET stability change cols are are split over
# are are split over multiple lines and Mutation information is empty! # multiple lines and Mutation information is empty!
mytext = mytext.replace('ange:\n', 'ange: ') mytext = mytext.replace('ange:\n', 'ange: ')
#print(mytext) # print(mytext)
# initiliase result_dict # initiliase result_dict
result_dict = {} result_dict = {}
for line in mytext.split('\n'): for line in mytext.split('\n'):
fields = line.split(':') fields = line.split(':')
#print(fields) # print(fields)
if len(fields) > 1: # since Mutaton information is empty if len(fields) > 1: # since Mutaton information is empty
dict_entry = dict([(x, y) for x, y in zip(fields[::2], fields[1::2])]) dict_entry = dict([(x, y) for x, y in zip(fields[::2], fields[1::2])])
result_dict.update(dict_entry) result_dict.update(dict_entry)
return result_dict return result_dict
#======================================================================= #=======================================================================
#%% call function #%% call function
#request_results(infile_url) #request_results(infile_url)