From 48773a19efdd6eb1372d9493c6dc4e9be4e45d85 Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Fri, 14 Aug 2020 16:41:11 +0100 Subject: [PATCH] tidy script for linking or_kinship with missense variant info --- scripts/find_missense.py | 40 +++--- scripts/or_kinship_link.py | 256 ++++++++++++++++++------------------- 2 files changed, 143 insertions(+), 153 deletions(-) diff --git a/scripts/find_missense.py b/scripts/find_missense.py index 71f2c5d..eb4789d 100755 --- a/scripts/find_missense.py +++ b/scripts/find_missense.py @@ -4,15 +4,15 @@ 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'): +#def find_missense(df, ref_allele1, alt_allele0): +def find_missense(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 + @df: df containing columns to compare @type: pandas df @ref_allele_column: column containing ref allele str @@ -37,7 +37,7 @@ def find_missense(test_df, ref_allele_column, alt_allele_column, n_diff_colname name will override original column @rtype: pandas df """ - for ind, val in test_df.iterrows(): + for ind, val in df.iterrows(): if DEBUG: print('index:', ind, 'value:', val , '\n============================================================') @@ -46,9 +46,9 @@ def find_missense(test_df, ref_allele_column, alt_allele_column, n_diff_colname 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 + 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 + 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') @@ -57,29 +57,29 @@ def find_missense(test_df, ref_allele_column, alt_allele_column, n_diff_colname 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 + ##df.at[ind, 'ref_allele'] = 'no_change' # adding column + ##df.at[ind, 'alt_allele'] = 'no_change' # adding column + df.at[ind, ref_a_colname] = 'no_change' # adding column + 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 + ##df.at[ind, 'ref_allele'] = ref_aln + ##df.at[ind, 'alt_allele'] = alt_aln + df.at[ind, ref_a_colname] = ref_aln + 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' + ##df.at[ind, 'ref_allele'] = 'ERROR_Not_nsSNP' + ##df.at[ind, 'alt_allele'] = 'ERROR_Not_nsSNP' + df.at[ind, ref_a_colname] = 'ERROR_Not_nsSNP' + df.at[ind, alt_a_colname] = 'ERROR_Not_nsSNP' - return test_df + return df #======================================== # a representative example #eg_df = pd.read_csv('pnca_assoc.txt', sep = '\t', nrows = 10, header = 0, index_col = False) @@ -90,7 +90,7 @@ def find_missense(test_df, ref_allele_column, alt_allele_column, n_diff_colname #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') +# find_missense(df = eg_df, ref_allele_column = 'ref_allele1', alt_allele_column = 'alt_allele0') # print(eg_df) #if __name__ == '__main__': diff --git a/scripts/or_kinship_link.py b/scripts/or_kinship_link.py index a0fe7f3..8c6b569 100755 --- a/scripts/or_kinship_link.py +++ b/scripts/or_kinship_link.py @@ -74,13 +74,8 @@ if not outdir: #======= # input #======= -info_filename = 'snp_info.txt' -snp_info = datadir + '/' + info_filename -print('Info file: ', snp_info - , '\n============================================================') - -#gene_info_filename = 'ns'+ gene.lower()+ '_snp_info.txt' # without headers -gene_info_filename = 'ns'+ gene.lower()+ '_snp_info.csv' +gene_info_filename = 'ns'+ gene.lower()+ '_snp_info.txt' +#gene_info_filename = 'ns'+ gene.lower()+ '_snp_info.csv' gene_info = indir + '/' + gene_info_filename print('gene info file: ', gene_info , '\n============================================================') @@ -102,57 +97,22 @@ print('Output file: ', outfile_or_kin #%% read files: preformatted using bash # or file: '...assoc.txt' # FIXME: call bash script from here - or_df = pd.read_csv(gene_or, sep = '\t', header = 0, index_col = False) # 182, 12 (without filtering for missense muts, it was 212 i.e we 30 muts weren't missense) or_df.head() or_df.columns -#%% snp_info file: master and gene specific ones + +#%% snp_info file: master and gene specific ones # gene info -#info_df2 = pd.read_csv('nssnp_info_pnca.txt', sep = '\t', header = 0) #303, 10 -info_df2 = pd.read_csv(gene_info, sep = ',', header = 0) #303, 10 +info_df2 = pd.read_csv(gene_info, sep = '\t', header = 0) #447, 10 +#info_df2 = pd.read_csv(gene_info, sep = ',', header = 0) #447 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 - , '\n*****RESULT*****') #65.7% - -# large file -#info_df = pd.read_csv('snp_info.txt', sep = '\t', header = None) #12010 -info_df = pd.read_csv(snp_info, sep = '\t') #12010 -info_df.columns -#info_df.columns = ['chromosome_number', 'ref_allele', 'alt_allele', 'snp_info'] #12009, 4 - -info_df['chromosome_number'].nunique() #10257 -mut_cover = (info_df['chromosome_number'].nunique()/info_df['chromosome_number'].count()) * 100 -print('*****RESULT*****' - ,'\nPercentage of mutations in pncA:', mut_cover - , '\n*****RESULT*****') #85.4% - -# extract unique chr position numbers -genomic_pos = info_df['chromosome_number'].unique() -genomic_pos_df = pd.DataFrame(genomic_pos, columns = ['chr_pos']) -genomic_pos_df.dtypes - -genomic_pos_min = info_df['chromosome_number'].min() -genomic_pos_max = info_df['chromosome_number'].max() - -# genomic coord for pnca coding region -cds_len = (end_cds-start_cds) + 1 -pred_prot_len = (cds_len/3) - 1 - -# mindblowing: difference b/w bitwise (&) and 'and' -# DO NOT want &: is this bit set to '1' in both variables? Is this what you want? -#if (genomic_pos_min <= start_cds) & (genomic_pos_max >= end_cds): -print('*****RESULT*****' - , '\nlength of coding region:', cds_len, 'bp' - , '\npredicted protein length:', pred_prot_len, 'aa' , '\n*****RESULT*****') -if genomic_pos_min <= start_cds and genomic_pos_max >= end_cds: - print ('PASS: coding region for gene included in snp_info.txt') -else: - sys.exit('FAIL: coding region for gene not included in info file snp_info.txt') - +# v6: 61.07 +# v4: 65.7% #%% 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 # when you try to merge ONLY using chromosome_number, then it messes up... and is WRONG. @@ -161,104 +121,101 @@ else: df_ncols = len(or_df.columns) print('Dim of df:',or_df.shape - , '\nExtracting missense muts as single letters') - -#find_missense(or_df, 'ref_allele1', 'alt_allele0') -# adds columns to df passed + , '\nExtracting missense muts as single letters from: find_missense()') find_missense(or_df, ref_allele_column = 'ref_allele1', alt_allele_column = 'alt_allele0') - print('Dim of revised df:', or_df.shape, ' after extraction of missense muts') # FIXME: import this from function -ncols_from_func = 4 -if or_df.shape[1] == df_ncols + ncols_from_func: +ncols_added_func = 4 +if or_df.shape[1] == df_ncols + ncols_added_func: print('PASS: Succesfuly extracted ref and alt alleles for missense muts') else: print('FAIL: No. of cols mismatch' ,'\nOriginal length:', df_ncols - , '\nExpected no. of cols:', df_ncols + ncols_from_func + , '\nExpected no. of cols:', df_ncols + ncols_added_func , '\nGot:', or_df.shape[1] , '\nCheck hardcoded value of ncols_add?') +if (or_df['tot_diff'] == 1).sum() == len(or_df) and (or_df['n_diff'] == 1).sum() == len(or_df) and or_df['n_diff'].equals(or_df['tot_diff']): + print('PASS: missene muts correctly extracted from source') +else: + print('FAIL: n_diff and tot_diff differ, check source data') sys.exit() -del(df_ncols, ncols_from_func) +del(df_ncols, ncols_added_func) -#%% TRY MERGE -# check dtypes -or_df.dtypes -info_df.dtypes -#or_df.info() +#%% check dtypes before merging +#or_df.dtypes +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 -genomic_pos_df['chr_pos'].isin(or_df['chromosome_number']).sum() #144 +#info_df2.dtypes +info_df2.info() -# check how many chr_num in or_df are in info_df: should be ALL of them -or_df['chromosome_number'].isin(genomic_pos_df['chr_pos']).sum() #182 +#%% perform merge: or_df and snp_info +print('Preparing dfs for merging... Finding common cols to merge') -# sanity check 2 -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') +# find common columns +#merging_cols = ['chromosome_number', 'ref_allele', 'alt_allele'] +merging_cols = or_df.columns[or_df.columns.isin(info_df2.columns)].to_list() + +print('No. of common cols identified:', len(merging_cols) + , '\nColumns to merge on:', merging_cols + , '\nChecking dtypes in merging_cols...' + , '\n=================================================') + +# make sure chromosome_number dtypes are consisent +or_df[merging_cols].dtypes == info_df2[merging_cols].dtypes + +# info_df2 contains multiple chromosome number in the column, so it is not +# possible to convert this to int. Therefore, converting to string in or_df column +if not (or_df[merging_cols].dtypes == info_df2[merging_cols].dtypes).all(): + print('Data types not same, converting chromsome_number to str in or_df') + or_df['chromosome_number'] = or_df['chromosome_number'].astype(str) + print('Checking after converting dtype in or_df') +if (or_df[merging_cols].dtypes == info_df2[merging_cols].dtypes).all(): + print('PASS: dfs ready to merge..') else: - sys.exit('FAIL: some genomic locs or_df chr number DO NOT have meta data in snp_info.txt') -#%% perform merge + print('FAIL: unable to make dtypes consistent required for merging! Check dtypes') + sys.exit() +# %% sanity check and perform merge #my_join = 'inner' #my_join = 'outer' my_join = 'left' #my_join = 'right' -merging_cols = ['chromosome_number', 'ref_allele', 'alt_allele'] +expected_cols = or_df.shape[1] + info_df2.shape[1] - len(merging_cols) -print('Merging 2 dfs: or_df and info_df using join type:', my_join +print('Merging 2 dfs: or_df and info_df' + , '\nJoin type:', my_join , '\nColumns to merge on:', merging_cols + , '\nExpected cols after merging:', expected_cols , '\n=================================================') -#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 = merging_cols, how = my_join, indicator = True) -dfm1['_merge'].value_counts() -# count no. of missense mutations ONLY -print('Expected no. of missense mutations:', dfm1.snp_info.str.count(r'(missense.*)').sum()) - -# Merge with info_df2 has this has extra columns due to bash preformatting -# These extra columns are just 'snp_info' column split on '|' -print('Merging with info_df2 as it has,', len(set(info_df2.columns).difference(info_df.columns)) - , 'extra columns relevant for downstream analyses:\n\n' - , set(info_df2.columns).difference(info_df.columns)) dfm2 = pd.merge(or_df, info_df2, on = merging_cols, how = my_join, indicator = True) dfm2['_merge'].value_counts() +expected_cols = expected_cols + 1 # due to indicator = T -# count no. of nan -print('No. of NA in dfm2:', dfm2['mut_type'].isna().sum()) +# count no. of nan in 'mut_type'. ideally should be 0 +if not dfm2['mut_type'].isna().sum() > 0: + print('Good merging, no NA detected') +else: + print('OR detected without metadata' + , '\nNo. of NA in dfm2:', dfm2['mut_type'].isna().sum() + , '\nWriting these to output file to check later with jody!') + dfm2_missing_info = dfm2[dfm2['mut_type'].isna()] + missing_or_metadata = "or_kinship_missing_metadata.csv" + outfile_missing_or_metadata = outdir + '/' + str(dfm2['mut_type'].isna().sum()) + '_' + missing_or_metadata + print('\noutput file:', outfile_missing_or_metadata) + dfm2_missing_info.to_csv(outfile_missing_or_metadata, index = False) + print('Finsihed writing file' + , '\nDim:', dfm2_missing_info.shape) -# drop nan from dfm2_mis +#PENDING Jody's reply +# !!!!!!!! +# drop nan from dfm2_mis as these are not useful +print('Dropping NAs before further processing...') dfm2_mis = dfm2[dfm2['mut_type'].notnull()] - -#%% sanity check -# count no. of missense muts - -#if len(dfm1) - dfm1.snp_info.str.count(r'(missense.*)').sum() == dfm2['mut_type'].isna().sum(): -if dfm2_mis.shape[0] == dfm1.snp_info.str.count(r'(missense.*)').sum(): - print('PASSED: numbers cross checked' - , '\nTotal no. of missense mutations:', dfm1.snp_info.str.count(r'(missense.*)').sum() - , '\nNo. of mutations falsely assumed to be missense:', len(dfm1) - dfm1.snp_info.str.count(r'(missense.*)').sum()) -else: - print('FAIL: numbers mismatch' - , '\Expected no. of rows:',dfm1.snp_info.str.count(r'(missense.*)').sum() - , '\nGot:', dfm2_mis.shape[0] - , '\nExpected no. of cols:', dfm1.shape[1] + len(set(info_df2.columns).difference(info_df.columns))-1) - -# two ways to filter to get only missense muts -test = dfm1[dfm1['snp_info'].str.count('missense.*')>0] -dfm1_mis = dfm1[dfm1['snp_info'].str.match('(missense.*)') == True] -test.equals(dfm1_mis) - -if dfm1_mis[['chromosome_number', 'ref_allele', 'alt_allele']].equals(dfm2_mis[['chromosome_number', 'ref_allele', 'alt_allele']]): - print('PASS: Further cross checks successful') -else: - sys.exit('FAIL: Second cross check unsuccessful!') - -del(test, dfm1_mis) +# !!!!!!!! #%% extract mut info into three cols df_ncols = len(dfm2_mis.columns) @@ -269,19 +226,20 @@ ncols_add = 0 if not 'wild_type' in dfm2_mis.columns: print('Extracting and adding column: wild_type' , '\n===============================================================') - dfm2_mis['wild_type'] = dfm2_mis['mut_info'].str.extract('(\w{1})>') + dfm2_mis['wild_type'] = dfm2_mis['mut_info_f1'].str.extract('(\w{1})>') ncols_add+=1 if not 'position' in dfm2_mis.columns: print('Extracting and adding column: position' , '\n===============================================================') - dfm2_mis['position'] = dfm2_mis['mut_info'].str.extract('(\d+)') + dfm2_mis['position'] = dfm2_mis['mut_info_f1'].str.extract('(\d+)') + #dfm2_mis['position'] = dfm2_mis[:,'mut_info_f1'].str.extract('(\d+)') ncols_add+=1 if not 'mutant_type' in dfm2_mis.columns: print('Extracting and adding column: mutant_type' , '\n================================================================') - dfm2_mis['mutant_type'] = dfm2_mis['mut_info'].str.extract('>\d+(\w{1})') + dfm2_mis['mutant_type'] = dfm2_mis['mut_info_f1'].str.extract('>\d+(\w{1})') ncols_add+=1 if not 'mutationinformation' in dfm2_mis.columns: @@ -293,7 +251,7 @@ if not 'mutationinformation' in dfm2_mis.columns: print('No. of cols added:', ncols_add) if len(dfm2_mis.columns) == df_ncols + ncols_add: - print('PASS: mcsm style muts present in df' + print('PASS: mcsm style muts added to df' , '\n===============================================================') else: print('FAIL: No. of cols mismatch' @@ -349,7 +307,7 @@ del(df_ncols, ncols_add) #%%==============================!!!!!!!======================================= #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 +dfm2_mis2 = dfm2_mis.T.drop_duplicates().T #changes dtypes in cols, only used for sanity check dup_cols = set(dfm2_mis.columns).difference(dfm2_mis2.columns) print('Total no of duplicate columns:', len(dup_cols) , '\nDuplicate columns identified:', dup_cols @@ -359,7 +317,7 @@ print('Total no of duplicate columns:', len(dup_cols) #print('removing duplicate columns: kept one of the dup_cols i.e tot_diff') df_ncols = dfm2_mis.shape[1] -print('Removing duplicate columns' +print('Removing', len(dup_cols), 'duplicate columns:', dup_cols , '\nOriginal dim:', dfm2_mis.shape) dfm2_mis.drop(list(dup_cols), axis = 1, inplace = True) @@ -379,7 +337,8 @@ else: del(df_ncols) #3b) other not useful columns -cols_to_drop = ['chromosome_text', 'n_diff', 'chr', 'symbol', '_merge' ] +print('Dropping other redundant or unnecessary columns...') +cols_to_drop = ['chromosome_text', 'n_diff', 'chr', '_merge' , 'mut_region' , 'reference_allele', 'alternate_allele'] df_ncols = dfm2_mis.shape[1] dfm2_mis.drop(cols_to_drop, axis = 1, inplace = True) #dfm2_mis.rename(columns = {'ref_allele1': 'reference_allele'}, inplace = True) @@ -403,21 +362,21 @@ del(df_ncols) print('Reordering', dfm2_mis.shape[1], 'columns' , '\n===============================================') -column_order = ['mutation', +#dfm2_mis.columns + +column_order = [#'mutation', 'mutationinformation', 'wild_type', 'position', 'mutant_type', - 'chr_num_allele', + #'chr_num_allele', 'ref_allele', 'alt_allele', - 'mut_info', + 'mut_info_f1', + 'mut_info_f2', 'mut_type', 'gene_id', - 'gene_number', - 'mut_region', - 'reference_allele', - 'alternate_allele', + 'gene_name', 'chromosome_number', #'afs 'af_kin', @@ -439,12 +398,14 @@ column_order = ['mutation', # 'n_diff', # 'tot_diff', 'n_miss', - 'wt_3let', - 'mt_3let'] + #'wt_3let', + #'mt_3let' + ] -if len(column_order) == dfm2_mis.shape[1]: +if (len(column_order) == dfm2_mis.shape[1] and pd.DataFrame(column_order).isin(dfm2_mis.columns).all()).all(): print('PASS: Column order generated for', len(column_order), 'columns' - , '\nApplying column order to df...' ) + , '\nColumn names match to perform reordering' + , '\nApplying column order to df...' ) orkin_linked = dfm2_mis[column_order] else: print('FAIL: Mismatch in no. of cols to reorder' @@ -459,12 +420,41 @@ if orkin_linked.shape == dfm2_mis.shape and set(orkin_linked.columns) == set(dfm else: sys.exit('FAIL: something went wrong when rearranging columns!') + +# converting position and chromosome number to numeric +orkin_linked.dtypes +#orkin_linked['chromosome_number'] = pd.to_numeric(orkin_linked['chromosome_number']) + +orkin_linked[['chromosome_number', 'position']] = orkin_linked[['chromosome_number', 'position']].apply(pd.to_numeric) +orkin_linked.dtypes + +# write frequency of position counts +orkin_pos = pd.DataFrame(orkin_linked['position']) +z = orkin_pos['position'].value_counts() +z1 = z.to_dict() +orkin_pos['or_pos_count'] = orkin_pos['position'].map(z1) +orkin_pos['or_pos_count'].value_counts() + +orkin_linked['position'] +foo = orkin_linked['position'].value_counts() + +# order df by position +orkin_linked_o = orkin_linked.sort_values(by = ['position']) +bar = orkin_linked_o['position'].value_counts() + +if (foo == bar).all(): + print('PASS: df reorderd by position for output' + , '\nWriting output file') +else: + print('FAIL: could not reorder by position') + sys.exit() + #%% write file print('\n=====================================================================' , '\nWriting output file:\n', outfile_or_kin , '\nNo. of rows:', len(dfm2_mis) , '\nNo. of cols:', len(dfm2_mis.columns)) -orkin_linked.to_csv(outfile_or_kin, index = False) +orkin_linked_o.to_csv(outfile_or_kin, index = False) #%% diff b/w allele0 and 1: or_df #https://stackoverflow.com/questions/40348541/pandas-diff-with-string