tidy script for linking or_kinship with missense variant info

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(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_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
+            ##df.at[ind, 'ref_allele'] = 'no_change' # adding column
-            ##test_df.at[ind, 'alt_allele'] = 'no_change' # adding column
+            ##df.at[ind, 'alt_allele'] = 'no_change' # adding column
-            test_df.at[ind, ref_a_colname] = 'no_change' # adding column
+            df.at[ind, ref_a_colname] = 'no_change' # adding column
-            test_df.at[ind, alt_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
+            ##df.at[ind, 'ref_allele'] = ref_aln
-            ##test_df.at[ind, 'alt_allele'] = alt_aln
+            ##df.at[ind, 'alt_allele'] = alt_aln
-            test_df.at[ind, ref_a_colname] = ref_aln
+            df.at[ind, ref_a_colname] = ref_aln
-            test_df.at[ind, alt_a_colname] = alt_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'
+            ##df.at[ind, 'ref_allele'] = 'ERROR_Not_nsSNP'
-            ##test_df.at[ind, 'alt_allele'] = 'ERROR_Not_nsSNP'
+            ##df.at[ind, 'alt_allele'] = 'ERROR_Not_nsSNP'
-            test_df.at[ind, ref_a_colname] = 'ERROR_Not_nsSNP'
+            df.at[ind, ref_a_colname] = 'ERROR_Not_nsSNP'
-            test_df.at[ind, alt_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__':

scripts/or_kinship_link.py

@@ -74,13 +74,8 @@ if not outdir:
 #=======
 # input
 #=======
-info_filename  = 'snp_info.txt'
+gene_info_filename  = 'ns'+ gene.lower()+ '_snp_info.txt'
-snp_info = datadir + '/' + info_filename
+#gene_info_filename  = 'ns'+ gene.lower()+ '_snp_info.csv'
-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 = 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 = '\t', header = 0) #447, 10 
-info_df2  = pd.read_csv(gene_info, sep = ',', header = 0) #303, 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:    
+# v6: 61.07
-    print ('PASS: coding region for gene included in snp_info.txt')
+# v4: 65.7%   
-else:
-    sys.exit('FAIL: coding region for gene not included in info file snp_info.txt')
-    
 #%%  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')
+      , '\nExtracting missense muts as single letters from: find_missense()')
-
-#find_missense(or_df, 'ref_allele1', 'alt_allele0')
-# adds columns to df passed
 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 
+ncols_added_func = 4 
-if or_df.shape[1] == df_ncols + ncols_from_func:
+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 before merging
-# check dtypes
+#or_df.dtypes
-or_df.dtypes
+or_df.info()
-info_df.dtypes
-#or_df.info()
 
-# pandas documentation where it mentions: "Pandas uses the object dtype for storing strings" 
+#info_df2.dtypes
-# check how many unique chr_num in info_df are in or_df
+info_df2.info()
-genomic_pos_df['chr_pos'].isin(or_df['chromosome_number']).sum() #144
 
-# check how many chr_num in or_df are in info_df: should be ALL of them
+#%% perform merge: or_df and snp_info
-or_df['chromosome_number'].isin(genomic_pos_df['chr_pos']).sum() #182
+print('Preparing dfs for merging... Finding common cols to merge')
 
-# sanity check 2
+# find common columns
-if or_df['chromosome_number'].isin(genomic_pos_df['chr_pos']).sum() == len(or_df):
+#merging_cols = ['chromosome_number', 'ref_allele', 'alt_allele']
-    print('PASS: all genomic locs in or_df have meta datain info.txt')
+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')
+    print('FAIL: unable to make dtypes consistent required for merging! Check dtypes')
-#%% perform merge
+    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
+# count no. of nan in 'mut_type'. ideally should be 0
-print('No. of NA in dfm2:', dfm2['mut_type'].isna().sum()) 
+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', 
+ 'gene_name', 
- 'mut_region', 
- 'reference_allele', 
- 'alternate_allele', 
  'chromosome_number',
  #'afs
  'af_kin',
@@ -439,11 +398,13 @@ column_order = ['mutation',
 # 'n_diff',
 # 'tot_diff',
  'n_miss',
- 'wt_3let',
+ #'wt_3let',
- 'mt_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'
+          , '\nColumn names match to perform reordering'
           , '\nApplying column order to df...' )
     orkin_linked = dfm2_mis[column_order]
 else:
@@ -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