1205 lines
47 KiB
Python
Executable file
1205 lines
47 KiB
Python
Executable file
#!/usr/bin/env python3
|
|
# -*- coding: utf-8 -*-
|
|
'''
|
|
Created on Tue Aug 6 12:56:03 2019
|
|
|
|
@author: tanu
|
|
'''
|
|
|
|
# FIXME: include error checking to enure you only
|
|
# concentrate on positions that have structural info?
|
|
|
|
# FIXME: import dirs.py to get the basic dir paths available
|
|
#=======================================================================
|
|
# TASK: extract ALL <gene> matched mutations from GWAS data
|
|
# Input data file has the following format: each row = unique sample id
|
|
# id,country,lineage,sublineage,drtype,drug,dr_muts_col,other_muts_col...
|
|
# 0,sampleID,USA,lineage2,lineage2.2.1,Drug-resistant,0.0,WT,gene_match<wt>POS<mut>; pncA_c.<wt>POS<mut>...
|
|
# where multiple mutations and multiple mutation types are separated by ';'.
|
|
# We are interested in the protein coding region i.e mutation with the<gene>_'p.' format.
|
|
# This script splits the mutations on the ';' and extracts protein coding muts only
|
|
# where each row is a separate mutation
|
|
# sample ids AND mutations are NOT unique, but the COMBINATION (sample id + mutation) = unique
|
|
|
|
# NOTE
|
|
#drtype is renamed to 'resistance' in the 35k dataset
|
|
|
|
# output files: all lower case
|
|
# 0) <gene>_common_ids.csv
|
|
# 1) <gene>_ambiguous_muts.csv
|
|
# 2) <gene>_mcsm_snps.csv
|
|
# 3) <gene>_metadata.csv
|
|
# 4) <gene>_all_muts_msa.csv
|
|
# 5) <gene>_mutational_positons.csv
|
|
|
|
# FIXME
|
|
## Make all cols lowercase
|
|
## change WildPos: wild_pos
|
|
## Add an extra col: wild_chain_pos
|
|
## output df: <gene>_linking_df.csv
|
|
#containing the following cols
|
|
#1. Mutationinformation
|
|
#2. wild_type
|
|
#3. position
|
|
#4. mutant_type
|
|
#5. chain
|
|
#6. wild_pos
|
|
#7. wild_chain_pos
|
|
|
|
|
|
#=======================================================================
|
|
#%% load libraries
|
|
import os, sys
|
|
import re
|
|
import pandas as pd
|
|
#import numpy as np
|
|
import argparse
|
|
#=======================================================================
|
|
#%% homdir and curr dir and local imports
|
|
homedir = os.path.expanduser('~')
|
|
# set working dir
|
|
os.getcwd()
|
|
os.chdir(homedir + '/git/LSHTM_analysis/scripts')
|
|
os.getcwd()
|
|
|
|
# import aa dict
|
|
from reference_dict import my_aa_dict # CHECK DIR STRUC THERE!
|
|
from tidy_split import tidy_split
|
|
#=======================================================================
|
|
#%% 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
|
|
args = arg_parser.parse_args()
|
|
#=======================================================================
|
|
#%% variable assignment: input and output paths & filenames
|
|
drug = args.drug
|
|
gene = args.gene
|
|
|
|
#drug = 'pyrazinamide'
|
|
#gene = 'pncA'
|
|
|
|
gene_match = gene + '_p.'
|
|
print('mut pattern for gene', gene, ':', gene_match)
|
|
|
|
nssnp_match = gene_match+'[A-Z]{3}[0-9]+[A-Z]{3}'
|
|
print('nsSNP for gene', gene, ':', nssnp_match)
|
|
|
|
wt_regex = gene_match.lower()+'(\w{3})'
|
|
print('wt regex:', wt_regex)
|
|
|
|
mut_regex = r'\d+(\w{3})$'
|
|
print('mt regex:', mut_regex)
|
|
|
|
pos_regex = r'(\d+)'
|
|
print('position regex:', pos_regex)
|
|
|
|
# building cols to extract
|
|
dr_muts_col = 'dr_mutations_' + drug
|
|
other_muts_col = 'other_mutations_' + drug
|
|
|
|
print('Extracting columns based on variables:\n'
|
|
, drug
|
|
, '\n'
|
|
, dr_muts_col
|
|
, '\n'
|
|
, other_muts_col
|
|
, '\n===============================================================')
|
|
#=======================================================================
|
|
#%% input and output dirs and files
|
|
#=======
|
|
# dirs
|
|
#=======
|
|
datadir = homedir + '/' + 'git/Data'
|
|
indir = datadir + '/' + drug + '/' + 'input'
|
|
outdir = datadir + '/' + drug + '/' + 'output'
|
|
|
|
#=======
|
|
# input
|
|
#=======
|
|
#in_filename_master_master = 'original_tanushree_data_v2.csv' #19k
|
|
in_filename_master = 'mtb_gwas_meta_v3.csv' #33k
|
|
infile_master = datadir + '/' + in_filename_master
|
|
print('Input file: ', infile_master
|
|
, '\n============================================================')
|
|
|
|
#=======
|
|
# output
|
|
#=======
|
|
# several output files: in respective sections at the time of outputting files
|
|
print('Output filename: in the respective sections'
|
|
, '\nOutput path: ', outdir
|
|
, '\n=============================================================')
|
|
|
|
#%%end of variable assignment for input and output files
|
|
#=======================================================================
|
|
#%% Read input file
|
|
master_data = pd.read_csv(infile_master, sep = ',')
|
|
|
|
# column names
|
|
#list(master_data.columns)
|
|
|
|
# extract elevant columns to extract from meta data related to the drug
|
|
if in_filename_master == 'original_tanushree_data_v2.csv':
|
|
meta_data = master_data[['id'
|
|
, 'country'
|
|
, 'lineage'
|
|
, 'sublineage'
|
|
, 'drtype' #19k only
|
|
, drug
|
|
, dr_muts_col
|
|
, other_muts_col]]
|
|
|
|
if in_filename_master == 'mtb_gwas_meta_v3.csv':
|
|
core_cols = ['id'
|
|
, 'country'
|
|
, 'country2'
|
|
, 'geographic_source'
|
|
, 'region'
|
|
, 'date'
|
|
, 'strain'
|
|
, 'lineage'
|
|
, 'sublineage' #drtype renamed to resistance
|
|
, 'resistance'
|
|
, 'location'
|
|
, 'host_body_site'
|
|
, 'environment_material'
|
|
, 'host_status'
|
|
, 'hiv_status'
|
|
, 'HIV_status'
|
|
, 'isolation_source']
|
|
|
|
variable_based_cols = [drug
|
|
, dr_muts_col
|
|
, other_muts_col]
|
|
|
|
cols_to_extract = core_cols + variable_based_cols
|
|
print('Extracting', len(cols_to_extract), 'columns from master data')
|
|
|
|
meta_data = master_data[cols_to_extract]
|
|
del(master_data, variable_based_cols, cols_to_extract)
|
|
|
|
print('Extracted meta data from filename:', in_filename_master
|
|
, '\nDim:', meta_data.shape)
|
|
|
|
# checks and results
|
|
total_samples = meta_data['id'].nunique()
|
|
print('RESULT: Total samples:', total_samples
|
|
, '\n===========================================================')
|
|
|
|
# counts NA per column
|
|
meta_data.isna().sum()
|
|
print('No. of NAs/column:' + '\n', meta_data.isna().sum()
|
|
, '\n===========================================================')
|
|
#
|
|
# glance
|
|
#meta_data.head()
|
|
#total_samples - NA pyrazinamide = ?
|
|
# 19K: 19265-6754 = 12511
|
|
# 33K: 33681 - 23823 = 9858
|
|
|
|
# equivalent of table in R
|
|
# drug counts: complete samples for OR calcs
|
|
meta_data[drug].value_counts()
|
|
print('RESULT: Sus and Res samples:\n', meta_data[drug].value_counts()
|
|
, '\n===========================================================')
|
|
|
|
#%%
|
|
# IMPORTANT sanity check:
|
|
# This is to find out how many samples have 1 and more than 1 mutation,so you
|
|
# can use it to check if your data extraction process for dr_muts
|
|
# and other_muts has worked correctly AND also to check the dim of the
|
|
# final formatted data.
|
|
# This will have: unique COMBINATION of sample id and <gene_match> mutations
|
|
|
|
#========
|
|
# First: counting <gene_match> mutations in dr_muts_col column
|
|
#========
|
|
print('Now counting WT &', nssnp_match, 'muts within the column:', dr_muts_col)
|
|
# drop na and extract a clean df
|
|
clean_df = meta_data.dropna(subset=[dr_muts_col])
|
|
|
|
# sanity check: count na
|
|
na_count = meta_data[dr_muts_col].isna().sum()
|
|
|
|
if len(clean_df) == (total_samples - na_count):
|
|
print('PASS: clean_df extracted: length is', len(clean_df)
|
|
, '\nNo.of NAs in', dr_muts_col, '=', na_count, '/', total_samples
|
|
, '\n==========================================================')
|
|
else:
|
|
sys.exit('FAIL: Could not drop NAs')
|
|
|
|
dr_gene_count = 0
|
|
wt = 0
|
|
id_dr = []
|
|
id2_dr = []
|
|
#nssnp_match_regex = re.compile(nssnp_match)
|
|
|
|
for i, id in enumerate(clean_df.id):
|
|
#print (i, id)
|
|
#id_dr.append(id)
|
|
#count_gene_dr = clean_df[dr_muts_col].iloc[i].count(gene_match) # can include stop muts
|
|
count_gene_dr = len(re.findall(nssnp_match, clean_df[dr_muts_col].iloc[i], re.IGNORECASE))
|
|
#print(count_gene_dr)
|
|
if count_gene_dr > 0:
|
|
id_dr.append(id)
|
|
if count_gene_dr > 1:
|
|
id2_dr.append(id)
|
|
#print(id, count_gene_dr)
|
|
dr_gene_count = dr_gene_count + count_gene_dr
|
|
count_wt = clean_df[dr_muts_col].iloc[i].count('WT')
|
|
wt = wt + count_wt
|
|
|
|
print('RESULTS:')
|
|
print('Total WT in dr_muts_col:', wt)
|
|
print('Total matches of', gene, 'SNP matches in', dr_muts_col, ':', dr_gene_count)
|
|
print('Total samples with > 1', gene, 'muts in dr_muts_col:', len(id2_dr) )
|
|
print('=================================================================')
|
|
|
|
del(clean_df, na_count, i, id, wt, id2_dr, count_gene_dr, count_wt)
|
|
|
|
#========
|
|
# Second: counting <gene_match> mutations in dr_muts_col column
|
|
#========
|
|
print('Now counting WT &', nssnp_match, 'muts within the column:', other_muts_col)
|
|
# drop na and extract a clean df
|
|
clean_df = meta_data.dropna(subset=[other_muts_col])
|
|
|
|
# sanity check: count na
|
|
na_count = meta_data[other_muts_col].isna().sum()
|
|
|
|
if len(clean_df) == (total_samples - na_count):
|
|
print('PASS: clean_df extracted: length is', len(clean_df)
|
|
, '\nNo.of NAs =', na_count, '/', total_samples
|
|
, '\n=========================================================')
|
|
else:
|
|
sys.exit('FAIL: Could not drop NAs')
|
|
|
|
other_gene_count = 0
|
|
wt_other = 0
|
|
id_other = []
|
|
id2_other = []
|
|
|
|
for i, id in enumerate(clean_df.id):
|
|
#print (i, id)
|
|
#id_other.append(id)
|
|
#count_gene_other = clean_df[other_muts_col].iloc[i].count(gene_match) # can include stop muts
|
|
count_gene_other = len(re.findall(nssnp_match, clean_df[other_muts_col].iloc[i], re.IGNORECASE))
|
|
if count_gene_other > 0:
|
|
id_other.append(id)
|
|
if count_gene_other > 1:
|
|
id2_other.append(id)
|
|
#print(id, count_gene_other)
|
|
other_gene_count = other_gene_count + count_gene_other
|
|
count_wt = clean_df[other_muts_col].iloc[i].count('WT')
|
|
wt_other = wt_other + count_wt
|
|
|
|
print('RESULTS:')
|
|
print('Total WT in other_muts_col:', wt_other)
|
|
print('Total matches of', gene, 'SNP matches in', other_muts_col, ':', other_gene_count)
|
|
print('Total samples with > 1', gene, 'muts in other_muts_col:', len(id2_other) )
|
|
print('=================================================================')
|
|
|
|
print('Predicting total no. of rows in the curated df:', dr_gene_count + other_gene_count
|
|
, '\n===================================================================')
|
|
expected_rows = dr_gene_count + other_gene_count
|
|
|
|
del(i, id, wt_other, clean_df, na_count, id2_other, count_gene_other, count_wt)
|
|
|
|
#%%
|
|
############
|
|
# extracting dr and other muts separately along with the common cols
|
|
#############
|
|
print('Extracting dr_muts from col:', dr_muts_col, 'with other meta_data')
|
|
print('muts to extract:', nssnp_match )
|
|
|
|
#===============
|
|
# dr mutations: extract gene_match entries with meta data and ONLY dr_muts col
|
|
#===============
|
|
if in_filename_master == 'original_tanushree_data_v2.csv':
|
|
meta_data_dr = meta_data[['id'
|
|
,'country'
|
|
,'lineage'
|
|
,'sublineage'
|
|
,'drtype'
|
|
, drug
|
|
, dr_muts_col]]
|
|
if in_filename_master == 'mtb_gwas_meta_v3.csv':
|
|
dr_based_cols = [drug, dr_muts_col]
|
|
cols_to_extract = core_cols + dr_based_cols
|
|
print('Extracting', len(cols_to_extract), 'columns from meta data')
|
|
meta_data_dr = meta_data[cols_to_extract]
|
|
del(dr_based_cols, cols_to_extract)
|
|
|
|
if meta_data_dr.shape[0] == len(meta_data) and meta_data_dr.shape[1] == (len(meta_data.columns)-1):
|
|
print('PASS: Dimensions match'
|
|
, '\n===============================================================')
|
|
else:
|
|
print('FAIL: Dimensions mismatch:'
|
|
, 'Expected dim:', len(meta_data), (len(meta_data.columns)-1)
|
|
, '\nGot:', meta_data_dr.shape
|
|
, '\n===============================================================')
|
|
sys.exit()
|
|
|
|
# Extract within this the nsSNPs for gene of interest using string match
|
|
#meta_gene_dr = meta_data_dr.loc[meta_data_dr[dr_muts_col].str.contains(gene_match, na = False)]
|
|
meta_gene_dr = meta_data_dr.loc[meta_data_dr[dr_muts_col].str.contains(nssnp_match, na = False, regex = True, case = False)]
|
|
print('gene_snp_match in dr:', len(meta_gene_dr))
|
|
|
|
dr_id = meta_gene_dr['id'].unique()
|
|
print('RESULT: No. of samples with dr muts in pncA:', len(dr_id))
|
|
|
|
if len(id_dr) == len(meta_gene_dr):
|
|
print('PASS: lengths match'
|
|
, '\n===============================================================')
|
|
else:
|
|
print('FAIL: length mismatch'
|
|
, '\nExpected len:', len(id_dr)
|
|
, '\nGot:', len(meta_gene_dr))
|
|
sys.exit()
|
|
|
|
dr_id = pd.Series(dr_id)
|
|
|
|
#=================
|
|
# other mutations: extract nssnp_match entries from other_muts_col
|
|
#==================
|
|
print('Extracting other_muts from:', other_muts_col,'with other meta_data')
|
|
print('muts to extract:', nssnp_match)
|
|
|
|
if in_filename_master == 'original_tanushree_data_v2.csv':
|
|
meta_data_other = meta_data[['id'
|
|
, 'country'
|
|
, 'lineage'
|
|
, 'sublineage'
|
|
, 'drtype'
|
|
, drug
|
|
, other_muts_col]]
|
|
|
|
if in_filename_master == 'mtb_gwas_meta_v3.csv':
|
|
other_based_cols = [drug, other_muts_col]
|
|
cols_to_extract = core_cols + other_based_cols
|
|
print('Extracting', len(cols_to_extract), 'columns from meta data')
|
|
meta_data_other = meta_data[cols_to_extract]
|
|
del(other_based_cols, cols_to_extract)
|
|
|
|
if meta_data_other.shape[0] == len(meta_data) and meta_data_other.shape[1] == (len(meta_data.columns)-1):
|
|
print('PASS: Dimensions match'
|
|
, '\n===============================================================')
|
|
else:
|
|
print('FAIL: Dimensions mismatch:'
|
|
, 'Expected dim:', len(meta_data), (len(meta_data.columns)-1)
|
|
, '\nGot:', meta_data_other.shape
|
|
, '\n===============================================================')
|
|
sys.exit()
|
|
|
|
# Extract within this nSSNP for gene of interest using string match
|
|
#meta_gene_other = meta_data_other.loc[meta_data_other[other_muts_col].str.contains(gene_match, na = False)]
|
|
meta_gene_other = meta_data_other.loc[meta_data_other[other_muts_col].str.contains(nssnp_match, na = False, regex = True, case = False)]
|
|
print('gene_snp_match in other:', len(meta_gene_other))
|
|
|
|
other_id = meta_gene_other['id'].unique()
|
|
print('RESULT: No. of samples with other muts:', len(other_id))
|
|
|
|
if len(id_other) == len(meta_gene_other):
|
|
print('PASS: lengths match'
|
|
, '\n==============================================================')
|
|
else:
|
|
print('FAIL: length mismatch'
|
|
, '\nExpected len:', len(id_other)
|
|
, '\nGot:', len(meta_gene_other))
|
|
sys.exit()
|
|
|
|
other_id = pd.Series(other_id)
|
|
#%% Find common IDs
|
|
print('Now extracting common_ids...')
|
|
common_mut_ids = dr_id.isin(other_id).sum()
|
|
print('RESULT: No. of common ids:', common_mut_ids)
|
|
|
|
# sanity checks
|
|
# check if True: should be since these are common
|
|
if dr_id.isin(other_id).sum() == other_id.isin(dr_id).sum():
|
|
print('PASS: Cross check on no. of common ids')
|
|
else:
|
|
sys.exit('FAIL: Cross check on no. of common ids failed')
|
|
|
|
# check if the common are indeed the same!
|
|
# bit of a tautology, but better safe than sorry!
|
|
common_ids = dr_id[dr_id.isin(other_id)]
|
|
common_ids = common_ids.reset_index()
|
|
common_ids.columns = ['index', 'id']
|
|
|
|
common_ids2 = other_id[other_id.isin(dr_id)]
|
|
common_ids2 = common_ids2.reset_index()
|
|
common_ids2.columns = ['index', 'id2']
|
|
|
|
# should be True
|
|
if common_ids['id'].equals(common_ids2['id2']):
|
|
print('PASS: Further cross checks on common ids')
|
|
else:
|
|
sys.exit('FAIL: Further cross checks on common ids')
|
|
|
|
# good sanity check: use it later to check gene_sample_counts
|
|
expected_gene_samples = (len(meta_gene_dr) + len(meta_gene_other) - common_mut_ids)
|
|
print('Expected no. of gene samples:', expected_gene_samples
|
|
, '\n=================================================================')
|
|
#%% write file
|
|
#print(outdir)
|
|
out_filename_cid = gene.lower() + '_common_ids.csv'
|
|
outfile_cid = outdir + '/' + out_filename_cid
|
|
|
|
print('Writing file:'
|
|
, '\nFile:', outfile_cid
|
|
, '\nNo. of rows:', len(common_ids)
|
|
, '\n=============================================================')
|
|
|
|
common_ids.to_csv(outfile_cid, index = False)
|
|
|
|
del(out_filename_cid)
|
|
|
|
# clear variables
|
|
del(dr_id, other_id, meta_data_dr, meta_data_other, common_ids, common_mut_ids, common_ids2)
|
|
|
|
#%% Now extract 'all' gene specific nsSNP mutations: i.e 'nssnp_match'
|
|
print('Extracting nsSNP match:', gene, 'mutations from cols:\n'
|
|
, dr_muts_col, 'and', other_muts_col, 'using string match:'
|
|
, '\n===================================================================')
|
|
#meta_gene_all = meta_data.loc[meta_data[dr_muts_col].str.contains(gene_match, na = False) | meta_data[other_muts_col].str.contains(gene_match, na = False) ]
|
|
meta_gene_all = meta_data.loc[meta_data[dr_muts_col].str.contains(nssnp_match, na = False, regex = True, case = False) | meta_data[other_muts_col].str.contains(nssnp_match, na = False, regex = True, case = False) ]
|
|
|
|
extracted_gene_samples = meta_gene_all['id'].nunique()
|
|
print('RESULT: actual no. of gene samples extracted:', extracted_gene_samples
|
|
, '\n===================================================================')
|
|
|
|
# sanity check: length of gene samples
|
|
print('Performing sanity check:')
|
|
if extracted_gene_samples == expected_gene_samples:
|
|
print('PASS: expected & actual no. of nssnp gene samples match'
|
|
, '\nNo. of gene samples:', len(meta_gene_all)
|
|
, '\n=========================================================')
|
|
else:
|
|
sys.exit('FAIL: Length mismatch in gene samples!')
|
|
|
|
# count NA in drug column
|
|
gene_na = meta_gene_all[drug].isna().sum()
|
|
print('No. of gene samples without', drug, 'testing:', gene_na)
|
|
|
|
# use it later to check number of complete samples from LF data
|
|
comp_gene_samples = len(meta_gene_all) - gene_na
|
|
print('Complete gene samples tested for', drug, ':', comp_gene_samples)
|
|
print('=================================================================')
|
|
|
|
# Comment: This is still dirty data since these
|
|
# are samples that have nsSNP muts, but can have others as well
|
|
# since the format for mutations is mut1; mut2, etc.
|
|
print('This is still dirty data: samples have ', nssnp_match, 'muts but may have others as well'
|
|
, '\nsince the format for mutations is mut1; mut2, etc.'
|
|
, '\n=============================================================')
|
|
|
|
print('Performing tidy_split(): to separate the mutations into indivdual rows')
|
|
|
|
#=========
|
|
# DF1: dr_muts_col
|
|
#=========
|
|
########
|
|
# tidy_split(): on dr_muts_col column and remove leading white spaces
|
|
########
|
|
col_to_split1 = dr_muts_col
|
|
print ('Firstly, applying tidy split on dr muts df', meta_gene_dr.shape
|
|
, '\ncolumn name to apply tidy_split():', col_to_split1
|
|
, '\n============================================================')
|
|
# apply tidy_split()
|
|
dr_WF0 = tidy_split(meta_gene_dr, col_to_split1, sep = ';')
|
|
# remove leading white space else these are counted as distinct mutations as well
|
|
dr_WF0[dr_muts_col] = dr_WF0[dr_muts_col].str.lstrip()
|
|
|
|
# extract only the samples/rows with nssnp_match
|
|
#dr_gene_WF0 = dr_WF0.loc[dr_WF0[dr_muts_col].str.contains(gene_match)]
|
|
dr_gene_WF0 = dr_WF0.loc[dr_WF0[dr_muts_col].str.contains(nssnp_match, regex = True, case = False)]
|
|
|
|
print('Lengths after tidy split and extracting', nssnp_match, 'muts:'
|
|
, '\nOld length:' , len(meta_gene_dr)
|
|
, '\nLength after split:', len(dr_WF0)
|
|
, '\nLength of nssnp df:', len(dr_gene_WF0)
|
|
, '\nExpected len:', dr_gene_count
|
|
, '\n=============================================================')
|
|
if len(dr_gene_WF0) == dr_gene_count:
|
|
print('PASS: length matches expected length'
|
|
, '\n===============================================================')
|
|
else:
|
|
sys.exit('FAIL: lengths mismatch')
|
|
|
|
# count the freq of 'dr_muts' samples
|
|
dr_muts_df = dr_gene_WF0 [['id', dr_muts_col]]
|
|
print('Dim of dr_muts_df:', dr_muts_df.shape)
|
|
|
|
# add freq column
|
|
dr_muts_df['dr_sample_freq'] = dr_muts_df.groupby('id')['id'].transform('count')
|
|
#dr_muts_df['dr_sample_freq'] = dr_muts_df.loc[dr_muts_df.groupby('id')].transform('count')
|
|
print('Revised dim of dr_muts_df:', dr_muts_df.shape)
|
|
|
|
c1 = dr_muts_df.dr_sample_freq.value_counts()
|
|
print('Counting no. of sample frequency:\n', c1
|
|
, '\n===================================================================')
|
|
|
|
# sanity check: length of gene samples
|
|
if len(dr_gene_WF0) == c1.sum():
|
|
print('PASS: WF data has expected length'
|
|
, '\nLength of dr_gene WFO:', c1.sum()
|
|
, '\n===============================================================')
|
|
else:
|
|
sys.exit('FAIL: length mismatch!')
|
|
|
|
# Important: Assign 'column name' on which split was performed as an extra column
|
|
# This is so you can identify if mutations are dr_type or other in the final df
|
|
dr_df = dr_gene_WF0.assign(mutation_info = dr_muts_col)
|
|
print('Dim of dr_df:', dr_df.shape
|
|
, '\n=============================================================='
|
|
, '\nEnd of tidy split() on dr_muts, and added an extra column relecting mut_category'
|
|
, '\n===============================================================')
|
|
|
|
#%%
|
|
#=========
|
|
# DF2: other_mutations_pyrazinamdie
|
|
#=========
|
|
########
|
|
# tidy_split(): on other_muts_col column and remove leading white spaces
|
|
########
|
|
col_to_split2 = other_muts_col
|
|
print ('applying second tidy split() separately on other muts df', meta_gene_other.shape
|
|
, '\ncolumn name to apply tidy_split():', col_to_split2
|
|
, '\n============================================================')
|
|
|
|
# apply tidy_split()
|
|
other_WF1 = tidy_split(meta_gene_other, col_to_split2, sep = ';')
|
|
# remove the leading white spaces in the column
|
|
other_WF1[other_muts_col] = other_WF1[other_muts_col].str.strip()
|
|
|
|
# extract only the samples/rows with nssnp_match
|
|
#other_gene_WF1 = other_WF1.loc[other_WF1[other_muts_col].str.contains(gene_match)]
|
|
other_gene_WF1 = other_WF1.loc[other_WF1[other_muts_col].str.contains(nssnp_match, regex = True, case = False)]
|
|
|
|
print('Lengths after tidy split and extracting', gene_match, 'muts:',
|
|
'\nOld length:' , len(meta_gene_other),
|
|
'\nLength after split:', len(other_WF1),
|
|
'\nLength of nssnp df:', len(other_gene_WF1),
|
|
'\nExpected len:', other_gene_count
|
|
, '\n=============================================================')
|
|
if len(other_gene_WF1) == other_gene_count:
|
|
print('PASS: length matches expected length'
|
|
, '\n===============================================================')
|
|
else:
|
|
sys.exit('FAIL: lengths mismatch')
|
|
|
|
# count the freq of 'other muts' samples
|
|
other_muts_df = other_gene_WF1 [['id', other_muts_col]]
|
|
print('Dim of other_muts_df:', other_muts_df.shape)
|
|
|
|
# add freq column
|
|
other_muts_df['other_sample_freq'] = other_muts_df.groupby('id')['id'].transform('count')
|
|
print('Revised dim of other_muts_df:', other_muts_df.shape)
|
|
|
|
c2 = other_muts_df.other_sample_freq.value_counts()
|
|
print('Counting no. of sample frequency:\n', c2)
|
|
print('=================================================================')
|
|
# sanity check: length of gene samples
|
|
if len(other_gene_WF1) == c2.sum():
|
|
print('PASS: WF data has expected length'
|
|
, '\nLength of other_gene WFO:', c2.sum()
|
|
, '\n===============================================================')
|
|
else:
|
|
sys.exit('FAIL: Length mismatch')
|
|
|
|
# Important: Assign 'column name' on which split was performed as an extra column
|
|
# This is so you can identify if mutations are dr_type or other in the final df
|
|
other_df = other_gene_WF1.assign(mutation_info = other_muts_col)
|
|
print('dim of other_df:', other_df.shape
|
|
, '\n==============================================================='
|
|
, '\nEnd of tidy split() on other_muts, and added an extra column relecting mut_category'
|
|
, '\n===============================================================')
|
|
|
|
#%%
|
|
#==========
|
|
# Concatentating the two dfs: equivalent of rbind in R
|
|
#==========
|
|
# Important: Change column names to allow concat:
|
|
# dr_muts.. & other_muts : 'mutation'
|
|
print('Now concatenating the two dfs by row'
|
|
, '\nFirst assigning a common colname: "mutation" to the col containing muts'
|
|
, '\nThis is done for both dfs'
|
|
, '\n===================================================================')
|
|
|
|
dr_df.columns
|
|
dr_df.rename(columns = {dr_muts_col: 'mutation'}, inplace = True)
|
|
dr_df.columns
|
|
|
|
other_df.columns
|
|
other_df.rename(columns = {other_muts_col: 'mutation'}, inplace = True)
|
|
other_df.columns
|
|
|
|
if len(dr_df.columns) == len(other_df.columns):
|
|
print('Checking dfs for concatening by rows:'
|
|
, '\nDim of dr_df:', dr_df.shape
|
|
, '\nDim of other_df:', other_df.shape
|
|
, '\nExpected nrows:', len(dr_df) + len(other_df)
|
|
, '\n=============================================================')
|
|
else:
|
|
sys.exit('FAIL: No. of cols mismatch for concatenating')
|
|
|
|
# checking colnames before concat
|
|
print('Checking colnames BEFORE concatenating the two dfs...')
|
|
if (set(dr_df.columns) == set(other_df.columns)):
|
|
print('PASS: column names match necessary for merging two dfs')
|
|
else:
|
|
sys.exit('FAIL: Colnames mismatch for concatenating!')
|
|
|
|
# concatenate (axis = 0): Rbind
|
|
print('Now appending the two dfs: Rbind')
|
|
gene_LF_comb = pd.concat([dr_df, other_df], ignore_index = True, axis = 0)
|
|
|
|
print('Finding stop mutations in concatenated df')
|
|
stop_muts = gene_LF_comb['mutation'].str.contains('\*').sum()
|
|
if stop_muts == 0:
|
|
print('PASS: No stop mutations detected')
|
|
else:
|
|
print('stop mutations detected'
|
|
, '\nNo. of stop muts:', stop_muts, '\n'
|
|
, gene_LF_comb.groupby(['mutation_info'])['mutation'].apply(lambda x: x[x.str.contains('\*')].count())
|
|
, '\nNow removing them')
|
|
|
|
gene_LF0_nssnp = gene_LF_comb[~gene_LF_comb['mutation'].str.contains('\*')]
|
|
print('snps only: by subtracting stop muts:', len(gene_LF0_nssnp))
|
|
|
|
gene_LF0 = gene_LF_comb[gene_LF_comb['mutation'].str.contains(nssnp_match, case = False)]
|
|
print('snps only by direct regex:', len(gene_LF0))
|
|
|
|
if gene_LF0_nssnp.equals(gene_LF0):
|
|
print('PASS: regex for extracting nssnp worked correctly & stop mutations successfully removed'
|
|
, '\nUsing the regex extracted df')
|
|
else:
|
|
sys.exit('FAIL: posssibly regex or no of stop mutations'
|
|
, 'Regex being used:', nssnp_match)
|
|
#sys.exit()
|
|
|
|
# checking colnames and length after concat
|
|
print('Checking colnames AFTER concatenating the two dfs...')
|
|
if (set(dr_df.columns) == set(gene_LF0.columns)):
|
|
print('PASS: column names match'
|
|
, '\n=============================================================')
|
|
else:
|
|
sys.exit('FAIL: Colnames mismatch AFTER concatenating')
|
|
|
|
print('Checking concatenated df')
|
|
if len(gene_LF0) == (len(dr_df) + len(other_df))- stop_muts :
|
|
print('PASS:length of df after concat match'
|
|
, '\n===============================================================')
|
|
else:
|
|
sys.exit('FAIL: length mismatch')
|
|
|
|
#%%
|
|
###########
|
|
# This is hopefully clean data, but just double check
|
|
# Filter LF data so that you only have
|
|
# mutations corresponding to nssnp_match (case insensitive)
|
|
# this will be your list you run OR calcs
|
|
###########
|
|
print('Length of gene_LF0:', len(gene_LF0)
|
|
, '\nThis should be what we need. But just double checking and extracting nsSNP for', gene
|
|
, '\nfrom LF0 (concatenated data) using case insensitive regex match:', nssnp_match)
|
|
|
|
gene_LF1 = gene_LF0[gene_LF0['mutation'].str.contains(nssnp_match, regex = True, case = False)]
|
|
|
|
if len(gene_LF0) == len(gene_LF1):
|
|
print('PASS: length of gene_LF0 and gene_LF1 match',
|
|
'\nConfirming extraction and concatenating worked correctly'
|
|
, '\n==========================================================')
|
|
else:
|
|
print('FAIL: BUT NOT FATAL!'
|
|
, '\ngene_LF0 and gene_LF1 lengths differ'
|
|
, '\nsuggesting error in extraction process'
|
|
, ' use gene_LF1 for downstreama analysis'
|
|
, '\n=========================================================')
|
|
|
|
print('Length of dfs pre and post processing...'
|
|
, '\ngene WF data (unique samples in each row):', extracted_gene_samples
|
|
, '\ngene LF data (unique mutation in each row):', len(gene_LF1)
|
|
, '\n=============================================================')
|
|
|
|
|
|
#%% sanity check for extraction
|
|
# This ought to pass if nsnsp_match happens right at the beginning of creating 'expected_rows'
|
|
print('Verifying whether extraction process worked correctly...')
|
|
if len(gene_LF1) == expected_rows:
|
|
print('PASS: extraction process performed correctly'
|
|
, '\nExpected length:', expected_rows
|
|
, '\nGot:', len(gene_LF1)
|
|
, '\nRESULT: Total no. of mutant sequences for logo plot:', len(gene_LF1)
|
|
, '\n=========================================================')
|
|
else:
|
|
print('FAIL: extraction process has bugs'
|
|
, '\nExpected length:', expected_rows
|
|
, '\nGot:', len(gene_LF1)
|
|
, '\nDebug please'
|
|
, '\n=========================================================')
|
|
#%%
|
|
print('Performing some more sanity checks...')
|
|
|
|
# From LF1: useful for OR counts
|
|
# no. of unique muts
|
|
distinct_muts = gene_LF1.mutation.value_counts()
|
|
print('Distinct genomic mutations:', len(distinct_muts))
|
|
|
|
# no. of samples contributing the unique muts
|
|
gene_LF1.id.nunique()
|
|
print('No.of samples contributing to distinct genomic muts:', gene_LF1.id.nunique())
|
|
|
|
# no. of dr and other muts
|
|
mut_grouped = gene_LF1.groupby('mutation_info').mutation.nunique()
|
|
print('No.of unique dr and other muts:\n', gene_LF1.groupby('mutation_info').mutation.nunique())
|
|
|
|
# sanity check
|
|
if len(distinct_muts) == mut_grouped.sum() :
|
|
print('PASS:length of LF1 is as expected'
|
|
, '\n===============================================================')
|
|
else:
|
|
print('FAIL: Mistmatch in count of muts'
|
|
, '\nExpected count:', len(distinct_muts)
|
|
, '\nActual count:', mut_grouped.sum()
|
|
, '\nMuts should be distinct within dr* and other* type'
|
|
, '\nInspecting...possibly ambiguous muts'
|
|
, '\nNo. of possible ambiguous muts:', mut_grouped.sum() - len(distinct_muts)
|
|
, '\n=========================================================')
|
|
|
|
muts_split = list(gene_LF1.groupby('mutation_info'))
|
|
dr_muts = muts_split[0][1].mutation
|
|
other_muts = muts_split[1][1].mutation
|
|
print('splitting muts by mut_info:', muts_split)
|
|
print('no.of dr_muts samples:', len(dr_muts))
|
|
print('no. of other_muts samples', len(other_muts))
|
|
|
|
#%%
|
|
# IMPORTANT: The same mutation cannot be classed as a drug AND 'others'
|
|
if dr_muts.isin(other_muts).sum() & other_muts.isin(dr_muts).sum() > 0:
|
|
print('WARNING: Ambiguous muts detected in dr_ and other_ mutation category'
|
|
, '\n===============================================================')
|
|
else:
|
|
print('PASS: NO ambiguous muts detected'
|
|
, '\nMuts are unique to dr_ and other_ mutation class'
|
|
, '\n=========================================================')
|
|
|
|
# inspect dr_muts and other muts
|
|
if dr_muts.isin(other_muts).sum() & other_muts.isin(dr_muts).sum() > 0:
|
|
print('Finding ambiguous muts...'
|
|
, '\n========================================================='
|
|
, '\nTotal no. of samples in dr_muts present in other_muts:', dr_muts.isin(other_muts).sum()
|
|
, '\nThese are:', dr_muts[dr_muts.isin(other_muts)]
|
|
, '\n========================================================='
|
|
, '\nTotal no. of samples in other_muts present in dr_muts:', other_muts.isin(dr_muts).sum()
|
|
, '\nThese are:\n', other_muts[other_muts.isin(dr_muts)]
|
|
, '\n=========================================================')
|
|
else:
|
|
sys.exit('Error: ambiguous muts present, but extraction failed. Debug!')
|
|
|
|
print('Counting no. of ambiguous muts...')
|
|
|
|
if dr_muts[dr_muts.isin(other_muts)].nunique() == other_muts[other_muts.isin(dr_muts)].nunique():
|
|
common_muts = dr_muts[dr_muts.isin(other_muts)].value_counts().keys().tolist()
|
|
print('Distinct no. of ambigiuous muts detected:'+ str(len(common_muts))
|
|
, '\nlist of ambiguous mutations (see below):', *common_muts, sep = '\n')
|
|
print('\n===========================================================')
|
|
else:
|
|
print('Error: ambiguous muts detected, but extraction failed. Debug!'
|
|
, '\nNo. of ambiguous muts in dr:'
|
|
, len(dr_muts[dr_muts.isin(other_muts)].value_counts().keys().tolist())
|
|
, '\nNo. of ambiguous muts in other:'
|
|
, len(other_muts[other_muts.isin(dr_muts)].value_counts().keys().tolist())
|
|
, '\n=========================================================')
|
|
|
|
#%% clear variables
|
|
del(id_dr, id_other, meta_data, meta_gene_dr, meta_gene_other, mut_grouped, muts_split, other_WF1, other_df, other_muts_df, other_gene_count, gene_LF0, gene_na)
|
|
|
|
del(c1, c2, col_to_split1, col_to_split2, comp_gene_samples, dr_WF0, dr_df, dr_muts_df, dr_gene_WF0, dr_gene_count, expected_gene_samples, other_gene_WF1)
|
|
|
|
#%%: write file: ambiguous muts
|
|
# uncomment as necessary
|
|
#print(outdir)
|
|
#dr_muts.to_csv('dr_muts.csv', header = True)
|
|
#other_muts.to_csv('other_muts.csv', header = True)
|
|
|
|
out_filename_ambig_muts = gene.lower() + '_ambiguous_muts.csv'
|
|
outfile_ambig_muts = outdir + '/' + out_filename_ambig_muts
|
|
print('Writing file: ambiguous muts'
|
|
, '\nFilename:', outfile_ambig_muts)
|
|
|
|
inspect = gene_LF1[gene_LF1['mutation'].isin(common_muts)]
|
|
inspect.to_csv(outfile_ambig_muts, index = False)
|
|
|
|
print('Finished writing:', out_filename_ambig_muts
|
|
, '\nNo. of rows:', len(inspect)
|
|
, '\nNo. of cols:', len(inspect.columns)
|
|
, '\nNo. of rows = no. of samples with the ambiguous muts present:'
|
|
, dr_muts.isin(other_muts).sum() + other_muts.isin(dr_muts).sum()
|
|
, '\n=============================================================')
|
|
|
|
del(out_filename_ambig_muts)
|
|
#%% end of data extraction and some files writing. Below are some more files writing.
|
|
#=============================================================================
|
|
#%% Formatting df: read aa dict and pull relevant info
|
|
print('Now some more formatting:'
|
|
, '\nread aa dict and pull relevant info'
|
|
, '\nformat mutations:'
|
|
, '\nsplit mutation into mCSM style muts: '
|
|
, '\nFormatting mutation in mCSM style format: {WT}<POS>{MUT}'
|
|
, '\nassign aa properties: adding 2 cols at a time for each prop'
|
|
, '\n===================================================================')
|
|
|
|
# BEWARE hardcoding : only works as we are adding aa prop once for wt and once for mut
|
|
# in each lookup cycle
|
|
ncol_mutf_add = 3 # mut split into 3 cols
|
|
ncol_aa_add = 2 # 2 aa prop add (wt & mut) in each mapping
|
|
|
|
#===========
|
|
# Split 'mutation' column into three: wild_type, position and
|
|
# mutant_type separately. Then map three letter code to one using
|
|
# reference_dict imported at the beginning.
|
|
# After importing, convert to mutation to lowercase for compatibility with dict
|
|
#===========
|
|
gene_LF1['mutation'] = gene_LF1.loc[:, 'mutation'].str.lower()
|
|
|
|
print('wt regex being used:', wt_regex
|
|
, '\nmut regex being used:', mut_regex
|
|
, '\nposition regex being used:', pos_regex)
|
|
|
|
mylen0 = len(gene_LF1.columns)
|
|
#=======
|
|
# 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
|
|
#=======
|
|
print('Adding', ncol_mutf_add, 'more cols:\n')
|
|
|
|
# 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_regex).squeeze()converts to a series that map works on
|
|
wt = gene_LF1['mutation'].str.extract(wt_regex).squeeze()
|
|
gene_LF1['wild_type'] = wt.map(lookup_dict)
|
|
#mut = gene_LF1['mutation'].str.extract('\d+(\w{3})$').squeeze()
|
|
mut = gene_LF1['mutation'].str.extract(mut_regex).squeeze()
|
|
gene_LF1['mutant_type'] = mut.map(lookup_dict)
|
|
|
|
# extract position info from mutation column separetly using string match
|
|
#gene_LF1['position'] = gene_LF1['mutation'].str.extract(r'(\d+)')
|
|
gene_LF1['position'] = gene_LF1['mutation'].str.extract(pos_regex)
|
|
|
|
mylen1 = len(gene_LF1.columns)
|
|
|
|
# sanity checks
|
|
print('checking if 3-letter wt&mut residue extraction worked correctly')
|
|
if wt.isna().sum() & mut.isna().sum() == 0:
|
|
print('PASS: 3-letter wt&mut residue extraction worked correctly:'
|
|
, '\nNo NAs detected:'
|
|
, '\nwild-type\n', wt
|
|
, '\nmutant-type\n', mut
|
|
, '\ndim of df:', gene_LF1.shape)
|
|
else:
|
|
print('FAIL: 3-letter wt&mut residue extraction failed'
|
|
, '\nNo NAs detected:'
|
|
, '\nwild-type\n', wt
|
|
, '\nmutant-type\n', mut
|
|
, '\ndim of df:', gene_LF1.shape)
|
|
|
|
if mylen1 == mylen0 + ncol_mutf_add:
|
|
print('PASS: successfully added', ncol_mutf_add, 'cols'
|
|
, '\nold length:', mylen0
|
|
, '\nnew len:', mylen1)
|
|
else:
|
|
print('FAIL: failed to add cols:'
|
|
, '\nold length:', mylen0
|
|
, '\nnew len:', mylen1)
|
|
|
|
# clear variables
|
|
del(k, v, wt, mut, lookup_dict)
|
|
|
|
#=========
|
|
# iterate through the dict, create a lookup dict that i.e
|
|
# lookup_dict = {three_letter_code: aa_prop_water}
|
|
# Do this for both wild_type and mutant as above.
|
|
#=========
|
|
print('Adding', ncol_aa_add, 'more cols:\n')
|
|
|
|
# initialise a sub dict that is lookup dict for three letter code to aa prop
|
|
# adding two more cols
|
|
lookup_dict = dict()
|
|
for k, v in my_aa_dict.items():
|
|
lookup_dict[k] = v['aa_prop_water']
|
|
#print(lookup_dict)
|
|
wt = gene_LF1['mutation'].str.extract(wt_regex).squeeze()
|
|
gene_LF1['wt_prop_water'] = wt.map(lookup_dict)
|
|
#mut = gene_LF1['mutation'].str.extract('\d+(\w{3})$').squeeze()
|
|
mut = gene_LF1['mutation'].str.extract(mut_regex).squeeze()
|
|
gene_LF1['mut_prop_water'] = mut.map(lookup_dict)
|
|
|
|
mylen2 = len(gene_LF1.columns)
|
|
|
|
# sanity checks
|
|
print('checking if 3-letter wt&mut residue extraction worked correctly')
|
|
if wt.isna().sum() & mut.isna().sum() == 0:
|
|
print('PASS: 3-letter wt&mut residue extraction worked correctly:'
|
|
, '\nNo NAs detected:'
|
|
, '\nwild-type\n', wt
|
|
, '\nmutant-type\n', mut
|
|
, '\ndim of df:', gene_LF1.shape)
|
|
else:
|
|
print('FAIL: 3-letter wt&mut residue extraction failed'
|
|
, '\nNo NAs detected:'
|
|
, '\nwild-type\n', wt
|
|
, '\nmutant-type\n', mut
|
|
, '\ndim of df:', gene_LF1.shape)
|
|
|
|
if mylen2 == mylen1 + ncol_aa_add:
|
|
print('PASS: successfully added', ncol_aa_add, 'cols'
|
|
, '\nold length:', mylen1
|
|
, '\nnew len:', mylen2)
|
|
else:
|
|
print('FAIL: failed to add cols:'
|
|
, '\nold length:', mylen1
|
|
, '\nnew len:', mylen2)
|
|
|
|
# clear variables
|
|
del(k, v, wt, mut, lookup_dict)
|
|
|
|
#========
|
|
# iterate through the dict, create a lookup dict that i.e
|
|
# lookup_dict = {three_letter_code: aa_prop_polarity}
|
|
# Do this for both wild_type and mutant as above.
|
|
#=========
|
|
print('Adding', ncol_aa_add, 'more cols:\n')
|
|
|
|
# initialise a sub dict that is lookup dict for three letter code to aa prop
|
|
# adding two more cols
|
|
lookup_dict = dict()
|
|
for k, v in my_aa_dict.items():
|
|
lookup_dict[k] = v['aa_prop_polarity']
|
|
#print(lookup_dict)
|
|
wt = gene_LF1['mutation'].str.extract(wt_regex).squeeze()
|
|
gene_LF1['wt_prop_polarity'] = wt.map(lookup_dict)
|
|
#mut = gene_LF1['mutation'].str.extract(r'\d+(\w{3})$').squeeze()
|
|
mut = gene_LF1['mutation'].str.extract(mut_regex).squeeze()
|
|
gene_LF1['mut_prop_polarity'] = mut.map(lookup_dict)
|
|
|
|
mylen3 = len(gene_LF1.columns)
|
|
|
|
# sanity checks
|
|
print('checking if 3-letter wt&mut residue extraction worked correctly')
|
|
if wt.isna().sum() & mut.isna().sum() == 0:
|
|
print('PASS: 3-letter wt&mut residue extraction worked correctly:'
|
|
, '\nNo NAs detected:'
|
|
, '\nwild-type\n', wt
|
|
, '\nmutant-type\n', mut
|
|
, '\ndim of df:', gene_LF1.shape)
|
|
else:
|
|
print('FAIL: 3-letter wt&mut residue extraction failed'
|
|
, '\nNo NAs detected:'
|
|
, '\nwild-type\n', wt
|
|
, '\nmutant-type\n', mut
|
|
, '\ndim of df:', gene_LF1.shape)
|
|
|
|
if mylen3 == mylen2 + ncol_aa_add:
|
|
print('PASS: successfully added', ncol_aa_add, 'cols'
|
|
, '\nold length:', mylen1
|
|
, '\nnew len:', mylen2)
|
|
else:
|
|
print('FAIL: failed to add cols:'
|
|
, '\nold length:', mylen1
|
|
, '\nnew len:', mylen2)
|
|
|
|
# clear variables
|
|
del(k, v, wt, mut, lookup_dict)
|
|
|
|
#========
|
|
# iterate through the dict, create a lookup dict that i.e
|
|
# lookup_dict = {three_letter_code: aa_calcprop}
|
|
# Do this for both wild_type and mutant as above.
|
|
#=========
|
|
print('Adding', ncol_aa_add, 'more cols:\n')
|
|
|
|
lookup_dict = dict()
|
|
for k, v in my_aa_dict.items():
|
|
lookup_dict[k] = v['aa_calcprop']
|
|
#print(lookup_dict)
|
|
wt = gene_LF1['mutation'].str.extract(wt_regex).squeeze()
|
|
gene_LF1['wt_calcprop'] = wt.map(lookup_dict)
|
|
mut = gene_LF1['mutation'].str.extract(mut_regex).squeeze()
|
|
gene_LF1['mut_calcprop'] = mut.map(lookup_dict)
|
|
|
|
mylen4 = len(gene_LF1.columns)
|
|
|
|
# sanity checks
|
|
print('checking if 3-letter wt&mut residue extraction worked correctly')
|
|
if wt.isna().sum() & mut.isna().sum() == 0:
|
|
print('PASS: 3-letter wt&mut residue extraction worked correctly:'
|
|
, '\nNo NAs detected:'
|
|
, '\nwild-type\n', wt
|
|
, '\nmutant-type\n', mut
|
|
, '\ndim of df:', gene_LF1.shape)
|
|
else:
|
|
print('FAIL: 3-letter wt&mut residue extraction failed'
|
|
, '\nNo NAs detected:'
|
|
, '\nwild-type\n', wt
|
|
, '\nmutant-type\n', mut
|
|
, '\ndim of df:', gene_LF1.shape)
|
|
|
|
if mylen4 == mylen3 + ncol_aa_add:
|
|
print('PASS: successfully added', ncol_aa_add, 'cols'
|
|
, '\nold length:', mylen3
|
|
, '\nnew len:', mylen4)
|
|
else:
|
|
print('FAIL: failed to add cols:'
|
|
, '\nold length:', mylen3
|
|
, '\nnew len:', mylen4)
|
|
|
|
# clear variables
|
|
del(k, v, wt, mut, lookup_dict)
|
|
|
|
########
|
|
# 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
|
|
#########
|
|
gene_LF1['mutationinformation'] = gene_LF1['wild_type'] + gene_LF1.position.map(str) + gene_LF1['mutant_type']
|
|
print('Created column: mutationinformation'
|
|
, '\n=====================================================================\n'
|
|
, gene_LF1.mutationinformation.head(10))
|
|
|
|
#%% Write file: mCSM muts
|
|
snps_only = pd.DataFrame(gene_LF1['mutationinformation'].unique())
|
|
snps_only.head()
|
|
# assign column name
|
|
snps_only.columns = ['mutationinformation']
|
|
|
|
# count how many positions this corresponds to
|
|
pos_only = pd.DataFrame(gene_LF1['position'].unique())
|
|
|
|
print('Checking NA in snps...')# should be 0
|
|
if snps_only.mutationinformation.isna().sum() == 0:
|
|
print ('PASS: NO NAs/missing entries for SNPs'
|
|
, '\n===============================================================')
|
|
else:
|
|
sys.exit('FAIL: SNP has NA, Possible mapping issues from dict?')
|
|
|
|
out_filename_mcsmsnps = gene.lower() + '_mcsm_snps.csv'
|
|
outfile_mcsmsnps = outdir + '/' + out_filename_mcsmsnps
|
|
|
|
print('Writing file: mCSM style muts'
|
|
, '\nFile:', outfile_mcsmsnps
|
|
, '\nmutation format (SNP): {WT}<POS>{MUT}'
|
|
, '\nNo. of distinct muts:', len(snps_only)
|
|
, '\nNo. of distinct positions:', len(pos_only)
|
|
, '\n=============================================================')
|
|
|
|
snps_only.to_csv(outfile_mcsmsnps, header = False, index = False)
|
|
|
|
print('Finished writing:', outfile_mcsmsnps
|
|
, '\nNo. of rows:', len(snps_only)
|
|
, '\nNo. of cols:', len(snps_only.columns)
|
|
, '\n=============================================================')
|
|
del(out_filename_mcsmsnps)
|
|
|
|
#%% Write file: gene_metadata (i.e gene_LF1)
|
|
# where each row has UNIQUE mutations NOT unique sample ids
|
|
out_filename_metadata = gene.lower() + '_metadata.csv'
|
|
outfile_metadata = outdir + '/' + out_filename_metadata
|
|
print('Writing file: LF formatted data'
|
|
, '\nFile:', outfile_metadata
|
|
, '\n============================================================')
|
|
|
|
gene_LF1.to_csv(outfile_metadata, header = True, index = False)
|
|
print('Finished writing:', outfile_metadata
|
|
, '\nNo. of rows:', len(gene_LF1)
|
|
, '\nNo. of cols:', len(gene_LF1.columns)
|
|
, '\n=============================================================')
|
|
del(out_filename_metadata)
|
|
|
|
#%% write file: mCSM style but with repitions for MSA and logo plots
|
|
all_muts_msa = pd.DataFrame(gene_LF1['mutationinformation'])
|
|
all_muts_msa.head()
|
|
# assign column name
|
|
all_muts_msa.columns = ['mutationinformation']
|
|
|
|
# make sure it is string
|
|
all_muts_msa.columns.dtype
|
|
|
|
# sort the column
|
|
all_muts_msa_sorted = all_muts_msa.sort_values(by = 'mutationinformation')
|
|
|
|
# create an extra column with protein name
|
|
#all_muts_msa_sorted = all_muts_msa_sorted.assign(fasta_name = '3PL1')
|
|
#all_muts_msa_sorted.head()
|
|
|
|
# rearrange columns so the fasta name is the first column (required for mutate.script)
|
|
#all_muts_msa_sorted = all_muts_msa_sorted[['fasta_name', 'mutationinformation']]
|
|
all_muts_msa_sorted.head()
|
|
|
|
print('Checking NA in snps...')# should be 0
|
|
if all_muts_msa.mutationinformation.isna().sum() == 0:
|
|
print ('PASS: NO NAs/missing entries for SNPs'
|
|
, '\n===============================================================')
|
|
else:
|
|
sys.exit('FAIL: SNP has NA, Possible mapping issues from dict?')
|
|
|
|
out_filename_msa = gene.lower() +'_all_muts_msa.csv'
|
|
outfile_msa = outdir + '/' + out_filename_msa
|
|
|
|
print('Writing file: mCSM style muts for msa',
|
|
'\nFile:', outfile_msa,
|
|
'\nmutation format (SNP): {WT}<POS>{MUT}',
|
|
'\nNo.of lines of msa:', len(all_muts_msa))
|
|
|
|
all_muts_msa_sorted.to_csv(outfile_msa, header = False, index = False)
|
|
|
|
print('Finished writing:', outfile_msa
|
|
, '\nNo. of rows:', len(all_muts_msa)
|
|
, '\nNo. of cols:', len(all_muts_msa.columns)
|
|
, '\n=============================================================')
|
|
|
|
del(out_filename_msa)
|
|
|
|
#%% write file for mutational positions
|
|
# count how many positions this corresponds to
|
|
pos_only = pd.DataFrame(gene_LF1['position'].unique())
|
|
# assign column name
|
|
pos_only.columns = ['position']
|
|
# make sure dtype of column position is int or numeric and not string
|
|
pos_only.position.dtype
|
|
pos_only['position'] = pos_only['position'].astype(int)
|
|
pos_only.position.dtype
|
|
|
|
# sort by position value
|
|
pos_only_sorted = pos_only.sort_values(by = 'position', ascending = True)
|
|
|
|
out_filename_pos = gene.lower() + '_mutational_positons.csv'
|
|
outfile_pos = outdir + '/' + out_filename_pos
|
|
|
|
print('Writing file: mutational positions'
|
|
, '\nFile:', outfile_pos
|
|
, '\nNo. of distinct positions:', len(pos_only_sorted)
|
|
, '\n=============================================================')
|
|
|
|
pos_only_sorted.to_csv(outfile_pos, header = True, index = False)
|
|
|
|
print('Finished writing:', outfile_pos
|
|
, '\nNo. of rows:', len(pos_only_sorted)
|
|
, '\nNo. of cols:', len(pos_only_sorted.columns)
|
|
, '\n=============================================================')
|
|
|
|
del(out_filename_pos)
|
|
#=======================================================================
|
|
print(u'\u2698' * 50,
|
|
'\nEnd of script: Data extraction and writing files'
|
|
'\n' + u'\u2698' * 50 )
|
|
#%% end of script
|