From d3d82623d27ad404a5143829eae480a0b7ba63f7 Mon Sep 17 00:00:00 2001
From: Tanushree Tunstall <tanu@tunstall.in>
Date: Thu, 9 Jul 2020 14:08:27 +0100
Subject: [PATCH] added consistent style scripts to format kd & rd values

---
 scripts/kd_df.py | 255 +++++++++++++++++++++++++++++++++++++++++++++++
 scripts/rd_df.py | 203 +++++++++++++++++++++++++++++++++++++
 2 files changed, 458 insertions(+)
 create mode 100755 scripts/kd_df.py
 create mode 100755 scripts/rd_df.py

diff --git a/scripts/kd_df.py b/scripts/kd_df.py
new file mode 100755
index 0000000..e36d660
--- /dev/null
+++ b/scripts/kd_df.py
@@ -0,0 +1,255 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+'''
+Created on Tue Aug  6 12:56:03 2019
+
+@author: tanu
+'''
+#=======================================================================
+# Task: Hydrophobicity (Kd) values for amino acid sequence using the
+# Kyt&-Doolittle.
+# Same output as using the expasy server (link below)
+# Input: fasta file
+    
+# Output: csv file with 
+
+# useful links
+# https://biopython.org/DIST/docs/api/Bio.SeqUtils.ProtParamData-pysrc.html
+# https://web.expasy.org/protscale/pscale/protscale_help.html
+#=======================================================================
+#%% load packages
+import sys, os
+import argparse
+import pandas as pd
+import numpy as np
+from pylab import *
+from Bio.SeqUtils import ProtParamData
+from Bio.SeqUtils.ProtParam import ProteinAnalysis
+from Bio import SeqIO
+#from Bio.Alphabet.IUPAC import IUPACProtein
+import pprint as pp
+#=======================================================================
+#%% specify homedir and curr dir
+homedir = os.path.expanduser('~')
+
+# set working dir
+os.getcwd()
+os.chdir(homedir + '/git/LSHTM_analysis/scripts')
+os.getcwd()
+#=======================================================================
+#%% command line args
+arg_parser = argparse.ArgumentParser()
+arg_parser.add_argument('-d', '--drug', help='drug name', default = None)
+arg_parser.add_argument('-g', '--gene', help='gene name', default = None)
+#arg_parser.add_argument('-p', '--plot', help='show plot',  action='store_true')
+
+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('-fasta','--fasta_file', help = 'fasta file. By default, it assmumes a file called <gene>.fasta.txt in input_dir') 
+
+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.'
+
+data_dir     = args.datadir
+indir        = args.input_dir
+outdir       = args.output_dir
+
+fasta_filename = args.fasta_file
+
+#plot = args.plot
+DEBUG        = args.debug
+
+#============
+# directories
+#============
+if data_dir:
+    datadir = data_dir
+else:
+    datadir = homedir + '/' + 'git/Data'
+
+if not indir:
+    indir = datadir + '/' + drug + '/' + 'input'
+
+if not outdir:
+    outdir = datadir + '/' + drug + '/' + 'output'
+
+#=======
+# input
+#=======
+if fasta_filename:
+    in_filename_fasta = fasta_filename
+else:
+    in_filename_fasta = gene.lower() + '.fasta.txt'
+    
+infile_fasta = indir + '/' + in_filename_fasta
+print('Input fasta file:', infile_fasta
+      , '\n============================================================')
+
+#=======
+# output 
+#=======
+out_filename_kd = gene.lower() + '_kd.csv'
+outfile_kd =  outdir + '/' + out_filename_kd
+print('Output file:', outfile_kd
+      , '\n=============================================================')
+#%% end of variable assignment for input and output files
+#=======================================================================
+#%% kd values from fasta file and output csv
+def kd_to_csv(inputfasta, outputkdcsv, windowsize = 3):
+    """
+    Calculate kd (hydropathy values) from input fasta file
+
+    @param inputfasta: fasta file
+    @type: string
+
+    @param outputkdcsv: csv file with kd values
+    @type: string
+
+    @param windowsize: windowsize to perform KD calcs on (Kyte&-Doolittle)
+    @type: numeric
+
+    @return: none, writes kd values df as csv
+    """
+    #========================
+    # read input fasta file
+    #========================
+    fh = open(inputfasta)
+
+    for record in SeqIO.parse(fh, 'fasta'):
+        id = record.id
+        seq = record.seq
+        num_residues = len(seq)
+    fh.close()
+
+    sequence = str(seq)
+    X = ProteinAnalysis(sequence)
+
+    #===================
+    # calculate KD values: same as the expasy server
+    #===================
+    my_window = windowsize
+    offset = round((my_window/2)-0.5)
+    # edge weight is set to  default (100%)
+    
+    kd_values = (X.protein_scale(ProtParamData.kd , window = my_window))
+    # sanity checks 
+    print('Sequence Length:', num_residues)
+    print('kd_values Length:',len(kd_values))
+    print('Window Length:', my_window)
+    print('Window Offset:', offset)
+    print('=================================================================')
+    print('Checking:len(kd values) is as expected for the given window size & offset...')
+    expected_length =  num_residues - (my_window - offset) 
+    if len(kd_values) == expected_length:
+        print('PASS: expected and actual length of kd values match')
+    else:
+        print('FAIL: length mismatch'
+              ,'\nExpected length:', expected_length
+              ,'\nActual length:', len(kd_values)
+              , '\n=========================================================')
+ 
+    #===================
+    # creating two dfs
+    #===================
+    # 1) aa sequence and 2) kd_values. Then reset index for each df 
+    # which will allow easy merging of the two dfs.
+
+    # df1: df of aa seq with index reset to start from 1 
+    # (reflective of the actual aa position in a sequence)
+    # Name column of wt as 'wild_type' to be the same name used 
+    # in the file required for merging later.
+    dfSeq = pd.DataFrame({'wild_type_kd':list(sequence)})
+    dfSeq.index = np.arange(1, len(dfSeq) + 1) # python is not inclusive
+
+    # df2: df of kd_values with index reset to start from offset + 1 and 
+    # subsequent matched length of the kd_values
+    dfVals = pd.DataFrame({'kd_values':kd_values})
+    dfVals.index = np.arange(offset + 1, len(dfVals) + 1 + offset)
+
+    # sanity checks
+    max(dfVals['kd_values'])
+    min(dfVals['kd_values'])
+    
+    #===================
+    # concatenating dfs
+    #===================
+    # Merge the two on index 
+    # (as these are now reflective of the aa position numbers): df1 and df2 
+    # This will introduce NaN where there is missing values. In our case this
+    # will be 2 (first and last ones based on window size and offset)
+
+    kd_df = pd.concat([dfSeq, dfVals], axis = 1)
+
+    #============================
+    # renaming index to position
+    #============================
+    kd_df = kd_df.rename_axis('position')
+    kd_df.head
+
+    print('Checking: position col i.e. index should be numeric')
+    if kd_df.index.dtype == 'int64':
+        print('PASS: position col is numeric'
+              , '\ndtype is:', kd_df.index.dtype)
+    else:
+        print('FAIL: position col is not numeric'
+              , '\nConverting to numeric')
+        kd_df.index.astype('int64')
+        print('Checking dtype for after conversion:\n'
+              , '\ndtype is:', kd_df.index.dtype
+              , '\n=========================================================')
+
+    # Ensuring lowercase column names for consistency
+    kd_df.columns = kd_df.columns.str.lower()
+    
+    #===============
+    # writing file
+    #===============
+    print('Writing file:'
+          , '\nFilename:', outputkdcsv
+          , '\nExpected no. of rows:', len(kd_df)
+          , '\nExpected no. of cols:', len(kd_df.columns)
+          , '\n=============================================================')
+
+    kd_df.to_csv(outputkdcsv, header = True, index = True)
+          
+    #===============
+    # plot: optional!
+    #===============
+    # http://www.dalkescientific.com/writings/NBN/plotting.html
+
+    # FIXME: save fig
+    # extract just pdb if from 'id' to pass to title of plot
+    # foo = re.match(r'(^[0-9]{1}\w{3})', id).groups(1)
+    #if doplot:
+    plot(kd_values, linewidth = 1.0)
+    #axis(xmin = 1, xmax = num_residues)
+    xlabel('Residue Number')
+    ylabel('Hydrophobicity')
+    title('K&D Hydrophobicity for ' + id)
+    show()
+        
+#%% end of function
+#=======================================================================
+def main():
+    print('Running hydropathy calcs with following params\n'
+        , '\nInput fasta file:', in_filename_fasta
+        , '\nOutput:', out_filename_kd)
+    kd_to_csv(infile_fasta, outfile_kd, 3)
+    print('Finished writing file:'
+        , '\nFile:', outfile_kd
+        , '\n=============================================================')
+    
+if __name__ == '__main__':
+    main()
+#%% end of script    
+#=======================================================================
diff --git a/scripts/rd_df.py b/scripts/rd_df.py
new file mode 100755
index 0000000..81863b5
--- /dev/null
+++ b/scripts/rd_df.py
@@ -0,0 +1,203 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+'''
+Created on Tue Aug  6 12:56:03 2019
+
+@author: tanu
+'''
+#=============================================================================
+# Task: Residue depth (rd) processing to generate a df with residue_depth(rd)
+# values
+
+# FIXME: source file is MANUALLY downloaded from the website
+# Input: '.tsv' i.e residue depth txt file (output from .zip file manually
+# downloaded from the website). 
+# This should be integrated into the pipeline
+
+# Output: .csv with 3 cols i.e position, rd_values & 3-letter wt aa code(caps)
+#=============================================================================
+#%% load packages
+import sys, os
+import argparse
+import pandas as pd
+#=============================================================================
+#%% specify input and curr dir
+homedir = os.path.expanduser('~')
+
+# set working dir
+os.getcwd()
+os.chdir(homedir + '/git/LSHTM_analysis/meta_data_analysis')
+os.getcwd()
+#=======================================================================
+#%% command line args
+arg_parser = argparse.ArgumentParser()
+arg_parser.add_argument('-d', '--drug', help='drug name', 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('-rd','--rd_file', help = 'residue depth file. By default, it assmumes a file called <gene>_rd.tsv in output_dir') 
+
+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.'
+
+data_dir     = args.datadir
+indir        = args.input_dir
+outdir       = args.output_dir
+
+rd_filename  = args.rd_file
+
+DEBUG        = args.debug
+
+#============
+# directories
+#============
+if data_dir:
+    datadir = data_dir
+else:
+    datadir = homedir + '/' + 'git/Data'
+
+if not indir:
+    indir = datadir + '/' + drug + '/' + 'input'
+
+if not outdir:
+    outdir = datadir + '/' + drug + '/' + 'output'
+
+#======
+# input
+#=======
+if rd_filename:
+    in_filename_rd = rd_filename
+else:
+    #in_filename_rd = '3pl1_rd.tsv'
+    in_filename_rd = gene.lower() + '_rd.tsv'
+ 
+infile_rd = outdir + '/' + in_filename_rd
+print('Input file:', infile_rd
+      , '\n=============================================================')
+      
+#=======
+# output 
+#=======
+out_filename_rd = gene.lower() + '_rd.csv'
+outfile_rd =  outdir + '/' + out_filename_rd
+print('Output file:', outfile_rd
+      , '\n=============================================================')
+      
+#%% end of variable assignment for input and output files
+#=======================================================================
+#%% rd values from <gene>_rd.tsv values
+def rd_to_csv(inputtsv, outputrdcsv):
+    """
+    formats residue depth values from input file
+
+    @param inputtsv: tsv file downloaded from {INSERT LINK}
+    @type inputtsv: string
+
+    @param outputrdsv: csv file with rd values
+    @type outfile_rd: string
+
+    @return: none, writes rd values df as csv
+    """
+    #========================
+    # read downloaded tsv file
+    #========================
+    #%%  Read input file
+    rd_data  = pd.read_csv(inputtsv, sep = '\t')  
+    print('Reading input file:', inputtsv
+          , '\nNo. of rows:', len(rd_data)
+          , '\nNo. of cols:', len(rd_data.columns))
+
+    print('Column names:', rd_data.columns
+        , '\n===========================================================')
+    #========================
+    # creating position col
+    #========================
+    # Extracting residue number from index and assigning 
+    # the values to a column [position]. Then convert the position col to numeric.
+    rd_data['position'] = rd_data.index.str.extract('([0-9]+)').values
+
+    # converting position to numeric      
+    rd_data['position'] = pd.to_numeric(rd_data['position'])
+    rd_data['position'].dtype
+
+    print('Extracted residue num from index and assigned as a column:'
+          , '\ncolumn name: position'
+          , '\ntotal no. of cols now:', len(rd_data.columns)
+          , '\n=========================================================')
+
+    #========================
+    # Renaming amino-acid 
+    # and all-atom cols
+    #========================
+    print('Renaming columns:'
+          , '\ncolname==> # chain:residue: wt_3letter_caps'
+          , '\nYES... the column name *actually* contains a # ..!'
+          , '\ncolname==> all-atom: rd_values'
+          , '\n=========================================================')
+
+    rd_data.rename(columns = {'# chain:residue':'wt_3letter_caps', 'all-atom':'rd_values'}, inplace = True)
+    print('Column names:', rd_data.columns)
+
+    #========================
+    # extracting df with the
+    # desired columns
+    #========================
+    print('Extracting relevant columns for writing df as csv')
+
+    rd_df = rd_data[['position','rd_values','wt_3letter_caps']]
+
+    if len(rd_df) == len(rd_data):
+        print('PASS: extracted df has expected no. of rows'
+              ,'\nExtracted df dim:'
+              ,'\nNo. of rows:', len(rd_df)
+              ,'\nNo. of cols:', len(rd_df.columns))
+    else:
+        print('FAIL: no. of rows mimatch'
+              , '\nExpected no. of rows:', len(rd_data)
+              , '\nGot no. of rows:', len(rd_df)
+              , '\n=====================================================')
+
+    # Ensuring lowercase column names for consistency
+    rd_df.columns = rd_df.columns.str.lower()
+
+    #===============
+    # writing file
+    #===============
+    print('Writing file:'
+          , '\nFilename:', outputrdcsv
+#          , '\nPath:',  outdir
+#          , '\nExpected no. of rows:', len(rd_df)
+#          , '\nExpected no. of cols:', len(rd_df.columns)
+          , '\n=========================================================')
+
+    rd_df.to_csv(outputrdcsv, header = True, index = False)
+
+#%% end of function
+#=======================================================================
+#%% call function
+#rd_to_csv(infile_rd, outfile_rd)
+#=======================================================================
+def main():
+    print('residue depth using the following params'
+        , '\nInput residue depth file:', in_filename_rd
+        , '\nOutput:', out_filename_rd)
+    rd_to_csv(infile_rd, outfile_rd)
+    print('Finished Writing file:'
+        , '\nFilename:', outfile_rd
+        , '\n=============================================================')
+
+if __name__ == '__main__':
+    main()
+#%% end of script    
+#=======================================================================