diff --git a/scripts/combine_afs_ors.py b/scripts/combine_afs_ors.py new file mode 100755 index 0000000..253260a --- /dev/null +++ b/scripts/combine_afs_ors.py @@ -0,0 +1,341 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +''' +Created on Tue Aug 6 12:56:03 2019 + +@author: tanu +''' +# FIXME: change filename 4 (mcsm normalised data) +# to be consistent like (pnca_complex_mcsm_norm.csv) : changed manually, but ensure this is done in the mcsm pipeline +#======================================================================= +# Task: combine 2 dfs with aa position as linking column +# This is done in 2 steps: +# merge 1: + +# useful link +# https://stackoverflow.com/questions/23668427/pandas-three-way-joining-multiple-dataframes-on-columns +#======================================================================= +#%% load packages +import sys, os +import pandas as pd +import numpy as np +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() + +# local import +from reference_dict import my_aa_dict # CHECK DIR STRUC THERE! +#======================================================================= +#%% 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) # case sensitive +args = arg_parser.parse_args() +#======================================================================= +#%% variable assignment: input and output +drug = 'pyrazinamide' +gene = 'pncA' +gene_match = gene + '_p.' + +# cmd variables +#drug = args.drug +#gene = args.gene + +#========== +# dir +#========== +datadir = homedir + '/' + 'git/Data' +outdir = datadir + '/' + drug + '/' + 'output' + +#======= +# input +#======= +in_filename_afor = gene.lower() + '_af_or.csv' +# FIXME +in_filename_afor_kin = gene.lower() + '_af_or_kinship.csv' +# needs to contain OR. it only has beta! + +infile1 = outdir + '/' + in_filename_afor +infile2 = outdir + '/' + in_filename_afor_kin + + +print('Input file1:', infile1 + , '\nInput file2:', infile2 + , '\n===================================================================') + +#======= +# output +#======= +out_filename = gene.lower() + '_metadata_afs_ors.csv' +outfile = outdir + '/' + out_filename +print('Output file:', outfile + , '\n===================================================================') + + +del(in_filename_afor, in_filename_afor_kin, outfile, datadir, outdir) +#%% end of variable assignment for input and output files +#======================================================================= +#%% format mutations +# mut_format: gene.abc1cde | 1A>1B + +#======================== +# read input csv files to combine +#======================== +afor_df = pd.read_csv(infile1, sep = ',') +afor_df_ncols = len(afor_df.columns) +afor_df_nrows = len(afor_df) +print('No. of rows in', infile1, ':', afor_df_nrows + , '\nNo. of cols in', infile1, ':', afor_df_ncols) + +afor_kin_df = pd.read_csv(infile2, sep = ',') +afor_kin_df_nrows = len(afor_kin_df) +afor_kin_df_ncols = len(afor_kin_df.columns) +print('No. of rows in', infile2, ':', afor_kin_df_nrows + , '\nNo. of cols in', infile2, ':', afor_kin_df_ncols) + +#======= +# Iterate through the dict, create a lookup dict i.e +# lookup_dict = {three_letter_code: one_letter_code}. +# lookup dict should be the key and the value (you want to create a column for) +# Then use this to perform the mapping separetly for wild type and mutant cols. +# The three letter code is extracted using a string match match from the dataframe and then converted +# to 'pandas series'since map only works in pandas series +#======= +gene_regex = gene_match.lower()+'(\w{3})' +print('gene regex being used:', gene_regex) + +# initialise a sub dict that is lookup dict for three letter code to 1-letter code +# adding three more cols +lookup_dict = dict() +for k, v in my_aa_dict.items(): + lookup_dict[k] = v['one_letter_code'] +# wt = gene_LF1['mutation'].str.extract('gene_p.(\w{3})').squeeze() # converts to a series that map works on + wt = afor_df['mutation'].str.extract(gene_regex).squeeze() + afor_df['wild_type'] = wt.map(lookup_dict) + mut = afor_df['mutation'].str.extract('\d+(\w{3})$').squeeze() + afor_df['mutant_type'] = mut.map(lookup_dict) + +# extract position info from mutation column separetly using string match +afor_df['position'] = afor_df['mutation'].str.extract(r'(\d+)') + +# combine the wild_type+poistion+mutant_type columns to generate +# mutationinformation (matches mCSM output field) +# Remember to use .map(str) for int col types to allow string concatenation + +afor_df['mutationinformation'] = afor_df['wild_type'] + afor_df['position'].map(str) + afor_df['mutant_type'] +print('Created column: mutationinformation' + , '\n=====================================================================' + , afor_df['mutationinformation'].head(10)) + +# sanity check +ncols_add = 4 # beware of hardcoding (3 cols for mcsm style mut + 1 for concatenating them all) +if len(afor_df.columns) == afor_df_ncols + ncols_add: + afor_df_ncols = len(afor_df.columns) # update afor_df_ncols after adding cols + print('PASS: successfully added', ncols_add, 'cols' + , '\nold length:', afor_df_ncols + , '\nnew length:', len(afor_df.columns)) +else: + print('FAIL: failed to add cols:' + , '\nExpected cols:', afor_df_ncols + ncols_add + , '\nGot:', len(afor_df.columns)) + sys.exit() +#%% Detect mutation format to see if you apply this func +# FIXME +#afor_df.iloc[[0]].str.match('pnca_') +#afor_df.dtypes + +#foo = afor_df.loc[:, afor_df.dtypes == object] + +genomic_mut_regex = gene_match.lower()+'\w{3}\d+\w{3}' +print('gene regex being used:', genomic_mut_regex) +afor_df[(afor_df == genomic_mut_regex).any(axis = 1)] + +#%% Finding common col to merge on +# Define merging column: multiple cols have been used for merge else the common cols +# get suffixes '_x' and '_y' attached +# also, couldn't include 'position' in merging_cols since data types don't match +merging_cols = ['wild_type', 'mutant_type', 'mutationinformation'] +ncommon_cols= len(merging_cols) + +# checking cross-over of mutations in the two dfs to merge +ndiff1 = afor_kin_df_nrows - afor_df['mutationinformation'].isin(afor_kin_df['mutationinformation']).sum() +print(ndiff1) +ndiff2 = afor_kin_df_nrows - afor_kin_df['mutationinformation'].isin(afor_df['mutationinformation']).sum() +print(ndiff2) + +# Define join type +#my_join = 'inner' +#my_join = 'right' +##my_join = 'left' +my_join = 'outer' + +# sanity check: how many muts from afor_kin_df are in afor_df. should be a complete subset +if ndiff2 == 0: + print('PASS: all muts in afor_kin_df are present in afor_df' + , '\nProceeding with combining the dfs...') + + combined_df = pd.merge(afor_df, afor_kin_df, on = merging_cols, how = my_join) + + if my_join == 'outer': + expected_rows = afor_df_nrows + ndiff1 + expected_cols = (afor_df_ncols + afor_kin_df_ncols) - ncommon_cols + if len(combined_df) == expected_rows and len(combined_df.columns) == expected_cols: + print('PASS: successfully combined dfs with:', my_join, 'join') + else: + print('FAIL: ', my_join, 'join') + + print('\nExpected no. of rows:', expected_rows + , '\nGot:', len(combined_df) + , '\nExpected no. of cols:', expected_cols + , '\nGot:', len(combined_df.columns)) + + elif my_join == 'inner': + expected_rows = afor_kin_df_nrows + expected_cols = (afor_df_ncols + afor_kin_df_ncols) - ncommon_cols + if len(combined_df) == expected_rows and len(combined_df.columns) == expected_cols: + print('PASS: successfully combined dfs with:', my_join, 'join') + else: + print('FAIL: ', my_join, 'join') + + print('\nExpected no. of rows:', expected_rows + , '\nGot:', len(combined_df) + , '\nExpected no. of cols:', expected_cols + , '\nGot:', len(combined_df.columns)) + + elif my_join == 'left': + expected_rows = afor_df_nrows + ndiff1 + expected_cols = (afor_df_ncols + afor_kin_df_ncols) - ncommon_cols + if len(combined_df) == expected_rows and len(combined_df.columns) == expected_cols: + print('PASS: successfully combined dfs with:', my_join, 'join') + else: + print('FAIL: ', my_join, 'join') + + print('\nExpected no. of rows:', expected_rows + , '\nGot:', len(combined_df) + , '\nExpected no. of cols:', expected_cols + , '\nGot:', len(combined_df.columns)) + + + elif my_join == 'right': + expected_rows = afor_kin_df_nrows + expected_cols = (afor_df_ncols + afor_kin_df_ncols) - ncommon_cols + if len(combined_df) == expected_rows and len(combined_df.columns) == expected_cols: + print('PASS: successfully combined dfs with:', my_join, 'join') + else: + print('FAIL: ', my_join, 'join') + + print('\nExpected no. of rows:', expected_rows + , '\nGot:', len(combined_df) + , '\nExpected no. of cols:', expected_cols + , '\nGot:', len(combined_df.columns)) + + else: + print('FAIL: failed to combine dfs, expected rows and cols not matched') +else: + print('FAIL: numbers mismatch, mutations present in afor_kin_df but not in afor_df') + + +#%% check duplicate cols: ones containing suffix '_x' or '_y' +# should only be position +foo = combined_df.filter(regex = r'.*_x|_y', axis = 1) +print(foo.columns) # should only be position + +# drop position col containing suffix '_y' and then rename col without suffix +combined_or_df = combined_df.drop(combined_df.filter(regex = r'.*_y').columns, axis = 1) +combined_or_df['position_x'].head() + +# renaming columns +combined_or_df.rename(columns = {'position_x': 'position'}, inplace = True) +combined_or_df['position'].head() + +# recheck +foo = combined_or_df.filter(regex = r'.*_x|_y', axis = 1) +print(foo.columns) # should only be position + +combined_or_df['af'].head() +combined_or_df.rename(columns = {'af': 'af_kin'}, inplace = True) +combined_or_df['af_kin'] +#%% calculate OR for kinship +combined_or_df['or_kin'] = np.exp(combined_or_df['beta']) + +# drop duplicate columns +#if combined_or_df['alternate_allele'].equals(combined_or_df['alt_allele0']): +# combined_or_df.drop('alternate_allele', axis = 1, inplace = True) + +combined_or_df2 = combined_or_df.T.drop_duplicates().T# changes dtypes in cols +dup_cols = set(combined_or_df.columns).difference(combined_or_df2.columns) +#tot_diff is equal to n_diff + +# drop some not required cols +combined_or_df.drop(list(dup_cols), axis = 1, inplace = True) +print(combined_or_df.columns) +combined_or_df.drop(['chromosome_text', 'chr', 'symbol', '_merge', ], axis = 1, inplace = True) + +combined_or_df.rename(columns = {'ref_allele1': 'reference_allele'}, inplace = True) + +print(combined_or_df.columns) +#%% reorder columns +#https://stackoverflow.com/questions/13148429/how-to-change-the-order-of-dataframe-columns + +# check af: curiosity + +# setting column's order +output_df = combined_or_df[['mutation', 'wild_type', 'position', 'mutant_type', 'mutationinformation' + , 'chr_num_allele', 'ref_allele', 'alt_allele' + , 'mut_info', 'mut_type', 'gene_id', 'gene_number', 'mut_region' + , 'reference_allele', 'alternate_allele', 'chromosome_number' + , 'afs', 'af_kin', 'ors_logistic', 'ors_chi_cus', 'or_kin', 'ors_fisher' + , 'pvals_logistic', 'pvals_fisher', 'p_wald', 'ci_lb_fisher', 'ci_ub_fisher' + , 'beta', 'se', 'logl_H1', 'l_remle','stat_chi', 'pvals_chi', 'n_diff' , 'n_miss']] + +#%% output combined or df +#=============== +# writing file +#=============== +print('Writing file...') + +#combined_or_df.to_csv(outfile, header = True, index = False) +output_df.to_csv(outfile, header = True, index = False) + +print('Finished writing file:', outfile + , '\nNo. of rows:', len(combined_or_df) + , '\nNo. of cols:', len(combined_or_df.columns) + , '\n=========================================================') + + + + + + + +#%% +df = pd.DataFrame() +column_names = ['x','y','z','mean'] +for col in column_names: + df[col] = np.random.randint(0,100, size=10000) +df.head() + + +# drop duplicate col with dup values not necessarily colnames +df['xdup'] = df['x'] +df + +df.T.drop_duplicates().T + +import math +math.exp(0) + +df['expX'] = np.exp(df['x']) # math doesn't understand series dtype +df + + + + + + diff --git a/scripts/find_missense.py b/scripts/find_missense.py new file mode 100755 index 0000000..51339cb --- /dev/null +++ b/scripts/find_missense.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python3 + +import pandas as pd + +DEBUG = False +#%% +#def find_missense(test_df, ref_allele1, alt_allele0): +def find_missense(test_df, ref_allele_column, alt_allele_column, n_diff_colname = 'n_diff', tot_diff_colname = 'tot_diff', ref_a_colname = 'ref_allele', alt_a_colname = 'alt_allele'): + + + """Find mismatches in pairwise comparison of strings b/w col_a and col_b + + Case insensitive, converts strings to uppercase before comparison + + @test_df: df containing columns to compare + @type: pandas df + + @ref_allele_column: column containing ref allele str + @type: str (converts to uppercase) + + @alt_allele_column: column containing alt_allele str + @type: str (converts to uppercase) + + @n_diff_colname: user defined colname for no. of char diff b/w ref_allele_str and alt_allele_str + @type: str + + @tot_diff_colname: user defined colname abs diff to indicate if strings are of equal length + @type: str + + @ref_a_colname: user defined colname containing extracted referece allele + @type: str + + @alt_a_colname: user defined colname containing extracted alt allele + @type: str + + returns df: with 4 columns. If column names clash, the function column + name will override original column + @rtype: pandas df + """ + for ind, val in test_df.iterrows(): + if DEBUG: + print('index:', ind, 'value:', val + , '\n============================================================') + ref_a = val[ref_allele_column].upper() + alt_a = val[alt_allele_column].upper() + if DEBUG: + print('ref_allele_string:', ref_a, 'alt_allele_string:', alt_a) + difference = sum(1 for e in zip(ref_a, alt_a) if e[0] != e[1]) + test_df.at[ind, n_diff_colname] = difference # adding column + tot_difference = difference + abs(len(ref_a) - len(alt_a)) + test_df.at[ind, tot_diff_colname] = tot_difference # adding column + if difference != tot_difference: + print('WARNING: lengths of ref_allele and alt_allele differ at index:', ind + , '\nNon-missense muts detected') + + # Now finding the mismatched char + ref_aln = '' + alt_aln = '' + if ref_a == alt_a: + ##test_df.at[ind, 'ref_allele'] = 'no_change' # adding column + ##test_df.at[ind, 'alt_allele'] = 'no_change' # adding column + test_df.at[ind, ref_a_colname] = 'no_change' # adding column + test_df.at[ind, alt_a_colname] = 'no_change' # adding column + elif len(ref_a) == len(alt_a) and len(ref_a) > 0: + print('ref:', ref_a, 'alt:', alt_a) + for n in range(len(ref_a)): + if ref_a[n] != alt_a[n]: + ref_aln += ref_a[n] + alt_aln += alt_a[n] + ##test_df.at[ind, 'ref_allele'] = ref_aln + ##test_df.at[ind, 'alt_allele'] = alt_aln + test_df.at[ind, ref_a_colname] = ref_aln + test_df.at[ind, alt_a_colname] = alt_aln + print('ref:', ref_aln) + print('alt:', alt_aln) + else: + ##test_df.at[ind, 'ref_allele'] = 'ERROR_Not_nsSNP' + ##test_df.at[ind, 'alt_allele'] = 'ERROR_Not_nsSNP' + test_df.at[ind, ref_a_colname] = 'ERROR_Not_nsSNP' + test_df.at[ind, alt_a_colname] = 'ERROR_Not_nsSNP' + + return test_df +#======================================== +# a representative example +eg_df = {'chromosome_number': [2288719, 2288766, 2288775, 2288779, 2288827, 1111111, 2222222], + 'ref_allele1': ['Tc', 'AG', 'AGCACCCTG', 'CCCTGGTGGCC', 'CACA', 'AA', 'CAT'], + 'alt_allele0': ['CC', 'CA', 'GGCACCCTGZ','TCCTGGTGGCCAAD', 'TACA', 'AA', 'TCZ']} + +# snippet of actual data +#eg_df = pd.read_csv('pnca_assoc.txt', sep = '\t', nrows = 10, header = 0, index_col = False) +eg_df = pd.DataFrame(eg_df) + +def main(): + #find_missense(eg_df, ref_allele1 = 'ref_allele', alt_allele0 = 'alt_allele') + find_missense(test_df = eg_df, ref_allele_column = 'ref_allele1', alt_allele_column = 'alt_allele0') + print(eg_df) + +if __name__ == '__main__': + main()