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

script to combine ors and afs

Browse source
This commit is contained in:
Tanushree Tunstall 2020-06-22 13:07:26 +01:00
parent c98ca7c8ae
commit 28e52d4194
2 changed files with 440 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

341
scripts/combine_afs_ors.py Executable file
View file

@ -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

99
scripts/find_missense.py Executable file
View file

@ -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()