combined and output all ors
This commit is contained in:
parent
d8b272b0ae
commit
a9498f8e08
2 changed files with 185 additions and 144 deletions
|
|
@ -14,7 +14,7 @@ Created on Wed Jun 10 11:13:49 2020
|
|||
#%% specify dirs
|
||||
import os, sys
|
||||
import pandas as pd
|
||||
#import numpy as np
|
||||
import numpy as np
|
||||
import re
|
||||
import argparse
|
||||
|
||||
|
|
@ -26,8 +26,8 @@ from find_missense import find_missense
|
|||
#%% 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) # case sensitive
|
||||
arg_parser.add_argument('-d', '--drug', help = 'drug name', default = 'pyrazinamide')
|
||||
arg_parser.add_argument('-g', '--gene', help = 'gene name (case sensitive)', default = 'pncA') # case sensitive
|
||||
|
||||
arg_parser.add_argument('-s', '--start_coord', help = 'start of coding region (cds) of gene', default = 2288681) # pnca cds
|
||||
arg_parser.add_argument('-e', '--end_coord', help = 'end of coding region (cds) of gene', default = 2289241) # pnca cds
|
||||
|
|
@ -37,6 +37,8 @@ args = arg_parser.parse_args()
|
|||
#%% variables
|
||||
#gene = 'pncA'
|
||||
#drug = 'pyrazinamide'
|
||||
#start_cds = 2288681
|
||||
#end_cds = 2289241
|
||||
|
||||
# cmd variables
|
||||
gene = args.gene
|
||||
|
|
@ -94,7 +96,7 @@ or_df.columns
|
|||
info_df2 = pd.read_csv(gene_info, sep = '\t', header = 0) #303, 10
|
||||
mis_mut_cover = (info_df2['chromosome_number'].nunique()/info_df2['chromosome_number'].count()) * 100
|
||||
print('*****RESULT*****'
|
||||
, '\nPercentage of missense mut in pncA:', mis_mut_cover
|
||||
, '\nPercentage of MISsense mut in pncA:', mis_mut_cover
|
||||
, '\n*****RESULT*****') #65.7%
|
||||
|
||||
# large file
|
||||
|
|
@ -117,8 +119,6 @@ genomic_pos_min = info_df['chromosome_number'].min()
|
|||
genomic_pos_max = info_df['chromosome_number'].max()
|
||||
|
||||
# genomic coord for pnca coding region
|
||||
#start_cds = 2288681
|
||||
#end_cds = 2289241
|
||||
cds_len = (end_cds-start_cds) + 1
|
||||
pred_prot_len = (cds_len/3) - 1
|
||||
|
||||
|
|
@ -134,6 +134,7 @@ if genomic_pos_min <= start_cds and genomic_pos_max >= end_cds:
|
|||
print ('PASS: coding region for gene included in snp_info.txt')
|
||||
else:
|
||||
print('FAIL: coding region for gene not included in info file snp_info.txt')
|
||||
sys.exit('ERROR: coding region of gene not included in the info file')
|
||||
|
||||
#%% Extracting ref allele and alt allele as single letters
|
||||
# info_df has some of these params as more than a single letter, which means that
|
||||
|
|
@ -162,7 +163,7 @@ del(orig_len, ncols_add)
|
|||
# check dtypes
|
||||
or_df.dtypes
|
||||
info_df.dtypes
|
||||
or_df.info()
|
||||
#or_df.info()
|
||||
|
||||
# pandas documentation where it mentions: "Pandas uses the object dtype for storing strings"
|
||||
# check how many unique chr_num in info_df are in or_df
|
||||
|
|
@ -175,23 +176,22 @@ or_df['chromosome_number'].isin(genomic_pos_df['chr_pos']).sum() #182
|
|||
if or_df['chromosome_number'].isin(genomic_pos_df['chr_pos']).sum() == len(or_df):
|
||||
print('PASS: all genomic locs in or_df have meta datain info.txt')
|
||||
else:
|
||||
print('FAIL: some genomic locs or_df chr number DO NOT have meta data in snp_info.txt')
|
||||
|
||||
sys.exit('FAIL: some genomic locs or_df chr number DO NOT have meta data in snp_info.txt')
|
||||
#%% Perform merge
|
||||
|
||||
#join_type = 'inner'
|
||||
#join_type = 'outer'
|
||||
join_type = 'left'
|
||||
#join_type = 'right'
|
||||
#my_join = 'inner'
|
||||
#my_join = 'outer'
|
||||
my_join = 'left'
|
||||
#my_join = 'right'
|
||||
|
||||
#dfm1 = pd.merge(or_df, info_df, on ='chromosome_number', how = join_type, indicator = True) # not unique!
|
||||
dfm1 = pd.merge(or_df, info_df, on = ['chromosome_number', 'ref_allele', 'alt_allele'], how = join_type, indicator = True)
|
||||
#dfm1 = pd.merge(or_df, info_df, on ='chromosome_number', how = my_join, indicator = True) # not unique!
|
||||
dfm1 = pd.merge(or_df, info_df, on = ['chromosome_number', 'ref_allele', 'alt_allele'], how = my_join, indicator = True)
|
||||
dfm1['_merge'].value_counts()
|
||||
|
||||
# count no. of missense mutations ONLY
|
||||
dfm1.snp_info.str.count(r'(missense.*)').sum()
|
||||
|
||||
dfm2 = pd.merge(or_df, info_df2, on = ['chromosome_number', 'ref_allele', 'alt_allele'], how = join_type, indicator = True)
|
||||
dfm2 = pd.merge(or_df, info_df2, on = ['chromosome_number', 'ref_allele', 'alt_allele'], how = my_join, indicator = True)
|
||||
dfm2['_merge'].value_counts()
|
||||
|
||||
# count no. of nan
|
||||
|
|
@ -241,12 +241,90 @@ else:
|
|||
, '\nGot no. of cols:', len(dfm2_mis.columns))
|
||||
sys.exit()
|
||||
|
||||
#%% formatting data for output
|
||||
print('no of cols preformatting data:', len(dfm2_mis.columns))
|
||||
|
||||
#1) Add column: OR for kinship calculated from beta coeff
|
||||
print('converting beta coeff to OR by exponent function\n:'
|
||||
, dfm2_mis['beta'].head())
|
||||
dfm2_mis['or_kin'] = np.exp(dfm2_mis['beta'])
|
||||
print(dfm2_mis['or_kin'].head())
|
||||
|
||||
#2) rename af column
|
||||
dfm2_mis.rename(columns = {'af': 'af_kin'
|
||||
, 'beta': 'beta_kin'
|
||||
, 'p_wald': 'pwald_kin'
|
||||
, 'se': 'se_kin', 'logl_H1': 'logl_H1_kin'
|
||||
, 'l_remle': 'l_remle_kin'}, inplace = True)
|
||||
|
||||
#3) drop some not required cols (including duplicate if you want)
|
||||
|
||||
#3a) drop duplicate columns
|
||||
dfm2_mis2 = dfm2_mis.T.drop_duplicates().T #changes dtypes in cols, so not used
|
||||
dup_cols = set(dfm2_mis.columns).difference(dfm2_mis2.columns)
|
||||
print('Duplicate columns identified:', dup_cols)
|
||||
dup_cols = {'alt_allele0', 'ps'} # didn't want to remove tot_diff
|
||||
|
||||
print('removing duplicate columns: kept one of the dup_cols i.e tot_diff')
|
||||
dfm2_mis.drop(list(dup_cols), axis = 1, inplace = True)
|
||||
print(dfm2_mis.columns)
|
||||
|
||||
#3b) other not useful columns
|
||||
dfm2_mis.drop(['chromosome_text', 'chr', 'symbol', '_merge', ], axis = 1, inplace = True)
|
||||
dfm2_mis.rename(columns = {'ref_allele1': 'reference_allele'}, inplace = True)
|
||||
|
||||
print(dfm2_mis.columns)
|
||||
|
||||
#4) reorder columns
|
||||
orkin_linked = dfm2_mis[['mutationinformation',
|
||||
'wild_type',
|
||||
'position',
|
||||
'mutant_type',
|
||||
'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',
|
||||
'or_kin',
|
||||
# 'ors_logistic',
|
||||
# 'ors_chi_cus',
|
||||
# 'ors_fisher',
|
||||
'pwald_kin',
|
||||
# 'pvals_logistic',
|
||||
# 'pvals_fisher',
|
||||
# 'ci_lb_fisher',
|
||||
# 'ci_ub_fisher' ,
|
||||
'beta_kin',
|
||||
'se_kin',
|
||||
'logl_H1_kin',
|
||||
'l_remle_kin',
|
||||
# 'stat_chi',
|
||||
# 'pvals_chi',
|
||||
'n_diff',
|
||||
'tot_diff',
|
||||
'n_miss']]
|
||||
|
||||
|
||||
# sanity check after reassigning columns
|
||||
if orkin_linked.shape == dfm2_mis.shape and set(orkin_linked.columns) == set(dfm2_mis.columns):
|
||||
print('PASS: Successfully formatted df with rearranged columns')
|
||||
else:
|
||||
sys.exit('FAIL: something went wrong when rearranging columns!')
|
||||
|
||||
#%% write file
|
||||
print('Writing output file:\n', outfile_or_kin
|
||||
print('\n====================================================================='
|
||||
, '\nWriting output file:\n', outfile_or_kin
|
||||
, '\nNo.of rows:', len(dfm2_mis)
|
||||
, '\nNo. of cols:', len(dfm2_mis.columns))
|
||||
dfm2_mis.to_csv(outfile_or_kin, index = False)
|
||||
orkin_linked.to_csv(outfile_or_kin, index = False)
|
||||
|
||||
#%% diff b/w allele0 and 1: or_df
|
||||
#https://stackoverflow.com/questions/40348541/pandas-diff-with-string
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue