From b4dbad7e54b6768f31fd9b5960017f991659b4a9 Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Mon, 23 Mar 2020 17:40:19 +0000 Subject: [PATCH] delete old file --- meta_data_analysis/pnca_data_extraction.py | 983 --------------------- 1 file changed, 983 deletions(-) delete mode 100755 meta_data_analysis/pnca_data_extraction.py diff --git a/meta_data_analysis/pnca_data_extraction.py b/meta_data_analysis/pnca_data_extraction.py deleted file mode 100755 index 72c7860..0000000 --- a/meta_data_analysis/pnca_data_extraction.py +++ /dev/null @@ -1,983 +0,0 @@ -#!/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? - -#%% load libraries -################### -# load libraries -import os, sys -import pandas as pd -import numpy as np - -#from pandas.api.types import is_string_dtype -#from pandas.api.types import is_numeric_dtype - -#======================================================== -# TASK: extract ALL pncA_p. mutations from GWAS data -# Input data file has the following format: each row = unique sample id -# id,country,lineage,sublineage,drtype,pyrazinamide,dr_mutations_pyrazinamide,other_mutations_pyrazinamide... -# 0,sampleID,USA,lineage2,lineage2.2.1,Drug-resistant,0.0,WT,pncA_p.POS; pncA_c.POS... -# where multiple mutations and multiple mutation types are separated by ';'. We are interested in the -# protein coding region i.e mutation with the 'p.' format. - -# the 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 - -# output files: -# 0) pnca_common_ids.csv -# 1) pnca_ambiguous_muts.csv -# 2) pnca_mcsm_snps.csv -# 3) pnca_metadata.csv -# 4) pnca_comp_snps.csv <---deleted> - -# 4) pnca_all_muts_msa.csv -# 5) pnca_mutational_positons.csv -#======================================================== -#%% specify homedir as python doesn't recognise tilde -homedir = os.path.expanduser('~') - -# my working dir -os.getcwd() -os.chdir(homedir + '/git/LSHTM_analysis/meta_data_analysis') -os.getcwd() - -# import aa dict -from reference_dict import my_aa_dict #CHECK DIR STRUC THERE! -#======================================================== - -#%% variable assignment: input and output paths & filenames - -drug = 'pyrazinamide' -gene = 'pncA' -gene_match = gene + '_p.' - -#======= -# input dir -#======= -#indir = 'git/Data/pyrazinamide/input/original' -indir = 'git/Data' + '/' + drug + '/' + 'input/original' -#========= -# output dir -#========= -# several output files -# output filenames in respective sections at the time of outputting files -#outdir = 'git/Data/pyrazinamide/output' -outdir = 'git/Data' + '/' + drug + '/' + 'output' - -#%%end of variable assignment for input and output files -#============================================================================== -#%% Read files - -in_filename = 'original_tanushree_data_v2.csv' -infile = homedir + '/' + indir + '/' + in_filename -print('Reading input master file:', infile) - -master_data = pd.read_csv(infile, sep = ',') - -# column names -#list(master_data.columns) - -# extract elevant columns to extract from meta data related to the pyrazinamide -meta_data = master_data[['id' - ,'country' - ,'lineage' - ,'sublineage' - ,'drtype' - , 'pyrazinamide' - , 'dr_mutations_pyrazinamide' - , 'other_mutations_pyrazinamide' - ]] - -del(master_data) - -# checks and results -total_samples = meta_data['id'].nunique() -print('RESULT: Total samples:', total_samples) -print('======================================================================') - -# counts NA per column -meta_data.isna().sum() -print('No. of NAs/column:' + '\n', meta_data.isna().sum()) -print('======================================================================') -# glance -meta_data.head() - -# equivalent of table in R -# pyrazinamide counts -meta_data.pyrazinamide.value_counts() -print('RESULT: Sus and Res samples:\n', meta_data.pyrazinamide.value_counts()) -print('======================================================================') - -# clear variables -del(indir, in_filename,infile) -#del(outdir) -#%% -# !!!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 pncA_p.mutations - -#======== -# First: counting pncA_p. mutations in dr_mutations_pyrazinamide column -#======== -print('Now counting WT & pncA_p. muts within the column: dr_mutations_pyrazinamide') - -# drop na and extract a clean df -clean_df = meta_data.dropna(subset=['dr_mutations_pyrazinamide']) - -# sanity check: count na -na_count = meta_data['dr_mutations_pyrazinamide'].isna().sum() - -if len(clean_df) == (total_samples - na_count): - print('PASS: clean_df extracted: length is', len(clean_df), - '\nNo.of NA s=', na_count, '/', total_samples) -else: - print('FAIL: dropping NA failed') -print('======================================================================') - -dr_pnca_count = 0 -wt = 0 -id_dr = [] -id2_dr = [] - -for i, id in enumerate(clean_df.id): -# print (i, id) -# id_dr.append(id) -# count_pnca_dr = clean_df.dr_mutations_pyrazinamide.iloc[i].count('pncA_p.') #works 2201 - count_pnca_dr = clean_df.dr_mutations_pyrazinamide.iloc[i].count(gene_match) #works 2201 - if count_pnca_dr > 0: - id_dr.append(id) - if count_pnca_dr > 1: - id2_dr.append(id) -# print(id, count_pnca_dr) - dr_pnca_count = dr_pnca_count + count_pnca_dr - count_wt = clean_df.dr_mutations_pyrazinamide.iloc[i].count('WT') - wt = wt + count_wt - -print('RESULTS:') -print('Total WT in dr_mutations_pyrazinamide:', wt) -print('Total matches of', gene_match, 'in dr_mutations_pyrazinamide:', dr_pnca_count) -print('Total samples with > 1', gene_match, 'muts in dr_mutations_pyrazinamide:', len(id2_dr) ) -print('======================================================================') - -del(i, id, wt, id2_dr, clean_df, na_count, count_pnca_dr, count_wt) - -#======== -# Second: counting pncA_p. mutations in dr_mutations_pyrazinamide column -#======== -print('Now counting WT & pncA_p. muts within the column: other_mutations_pyrazinamide') - -# drop na and extract a clean df -clean_df = meta_data.dropna(subset=['other_mutations_pyrazinamide']) - -# sanity check: count na -na_count = meta_data['other_mutations_pyrazinamide'].isna().sum() - -if len(clean_df) == (total_samples - na_count): - print('PASS: clean_df extracted: length is', len(clean_df), - '\nNo.of NA s=', na_count, '/', total_samples) -else: - print('FAIL: dropping NA failed') -print('======================================================================') - -other_pnca_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_pnca_other = clean_df.other_mutations_pyrazinamide.iloc[i].count('pncA_p.') - count_pnca_other = clean_df.other_mutations_pyrazinamide.iloc[i].count(gene_match) - if count_pnca_other > 0: - id_other.append(id) - if count_pnca_other > 1: - id2_other.append(id) -# print(id, count_pnca_other) - other_pnca_count = other_pnca_count + count_pnca_other - count_wt = clean_df.other_mutations_pyrazinamide.iloc[i].count('WT') - wt_other = wt_other + count_wt -print('RESULTS:') -print('Total WT in other_mutations_pyrazinamide:', wt_other) -print('Total matches of', gene_match, 'in other_mutations_pyrazinamide:', other_pnca_count) -print('Total samples with > 1', gene_match, 'muts in other_mutations_pyrazinamide:', len(id2_other) ) -print('======================================================================') - -print('Predicting total no. of rows in your curated df:', dr_pnca_count + other_pnca_count ) -expected_rows = dr_pnca_count + other_pnca_count - -del(i, id, wt_other, clean_df, na_count, id2_other, count_pnca_other, count_wt) - -#%% -############ -# extracting dr and other muts separately along with the common cols -############# -print('======================================================================') -print('Extracting dr_muts in a dr_mutations_pyrazinamide with other meta_data') -print('gene to extract:', gene_match ) - -#=============== -# dr mutations: extract pncA_p. entries with meta data and ONLY dr_muts col -#=============== -# FIXME: replace pyrazinamide with variable containing the drug name -# !!! important !!! -meta_data_dr = meta_data[['id' - ,'country' - ,'lineage' - ,'sublineage' - ,'drtype' - , 'pyrazinamide' - , 'dr_mutations_pyrazinamide' - ]] -print("expected dim should be:", len(meta_data), (len(meta_data.columns)-1) ) -print("actual dim:", meta_data_dr.shape ) -print('======================================================================') - -# Extract within this the gene of interest using string match -#meta_pnca_dr = meta_data.loc[meta_data.dr_mutations_pyrazinamide.str.contains('pncA_p.*', na = False)] -meta_pnca_dr = meta_data_dr.loc[meta_data_dr.dr_mutations_pyrazinamide.str.contains(gene_match, na = False)] - -dr_id = meta_pnca_dr['id'].unique() -print('RESULT: No. of samples with dr muts in pncA:', len(dr_id)) -print('checking RESULT:', - '\nexpected len =', len(id_dr), - '\nactual len =', len(meta_pnca_dr) ) - -if len(id_dr) == len(meta_pnca_dr): - print('PASS: lengths match') -else: - print('FAIL: length mismatch') -print('======================================================================') - -dr_id = pd.Series(dr_id) - -#================= -# other mutations: extract pncA_p. entries -#================== -print('======================================================================') -print('Extracting dr_muts in a other_mutations_pyrazinamide with other meta_data') -# FIXME: replace pyrazinamide with variable containing the drug name -# !!! important !!! -meta_data_other = meta_data[['id' - ,'country' - ,'lineage' - ,'sublineage' - ,'drtype' - , 'pyrazinamide' - , 'other_mutations_pyrazinamide' - ]] - -print("expected dim should be:", len(meta_data), (len(meta_data.columns)-1) ) -print("actual dim:", meta_data_other.shape ) -print('======================================================================') - -# Extract within this the gene of interest using string match -meta_pnca_other = meta_data_other.loc[meta_data_other.other_mutations_pyrazinamide.str.contains(gene_match, na = False)] - -other_id = meta_pnca_other['id'].unique() -print('RESULT: No. of samples with other muts:', len(other_id)) -print('checking RESULT:', - '\nexpected len =', len(id_other), - '\nactual len =', len(meta_pnca_other) ) - -if len(id_other) == len(meta_pnca_other): - print('PASS: lengths match') -else: - print('FAIL: length mismatch') -print('======================================================================') - -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 -dr_id.isin(other_id).sum() == other_id.isin(dr_id).sum() - -# check if the 24 Ids that are 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 -print(common_ids['id'].equals(common_ids2['id2'])) - -# good sanity check: use it later to check pnca_sample_counts -expected_pnca_samples = ( len(meta_pnca_dr) + len(meta_pnca_other) - common_mut_ids ) -print("expected no. of pnca samples:", expected_pnca_samples) -print('======================================================================') -#%% write file -#print(outdir) -out_filename0 = gene.lower() + '_' + 'common_ids.csv' -outfile0 = homedir + '/' + outdir + '/' + out_filename0 - -#FIXME: CHECK line len(common_ids) -print('Writing file: common ids:\n', - '\nFilename:', out_filename0, - '\nPath:', homedir +'/'+ outdir, - '\nExpected no. of rows:', len(common_ids) ) - -common_ids.to_csv(outfile0) -print('======================================================================') -del(out_filename0) - - -# clear variables -del(dr_id, other_id, meta_data_dr, meta_data_other, common_ids, common_mut_ids, common_ids2) - -#%% Now extract "all" pncA mutations: i.e 'pncA_p.*' -print("extracting all pncA mutations from dr_ and other_ cols using string match:", gene_match) -#meta_pnca_all = meta_data.loc[meta_data.dr_mutations_pyrazinamide.str.contains(gene_match) | meta_data.other_mutations_pyrazinamide.str.contains(gene_match) ] -meta_pnca_all = meta_data.loc[meta_data.dr_mutations_pyrazinamide.str.contains(gene_match, na = False) | meta_data.other_mutations_pyrazinamide.str.contains(gene_match, na = False) ] -print('======================================================================') - -extracted_pnca_samples = meta_pnca_all['id'].nunique() -print("RESULT: actual no. of pnca samples extracted:", extracted_pnca_samples) -print('======================================================================') - -# sanity check: length of pnca samples -print('Performing sanity check:') -if extracted_pnca_samples == expected_pnca_samples: - print('No. of pnca samples:', len(meta_pnca_all), - '\nPASS: expected & actual no. of pnca samples match') -else: - print("FAIL: Debug please!") -print('======================================================================') - -# count NA in pyrazinamide column -pnca_na = meta_pnca_all['pyrazinamide'].isna().sum() -print("No. of pnca samples without pza testing i.e NA in pza column:",pnca_na) - -# use it later to check number of complete samples from LF data -comp_pnca_samples = len(meta_pnca_all) - pnca_na -print("comp pnca samples tested for pza:", comp_pnca_samples) -print('======================================================================') - -# Comment: This is still dirty data since these -# are samples that have pncA_p. muts, but can have others as well -# since the format for mutations is mut1; mut2, etc. -print('This is still dirty data: samples have pncA_p. muts, but may have others as well', - '\nsince the format for mutations is mut1; mut2, etc.') -print('======================================================================') - -#%% tidy_split():Function to split mutations on specified delimiter: ';' -#https://stackoverflow.com/questions/41476150/removing-space-from-dataframe-columns-in-pandas - -print('Performing tidy_spllit(): to separate the mutations into indivdual rows') -# define the split function -def tidy_split(df, column, sep='|', keep=False): - """ - Split the values of a column and expand so the new DataFrame has one split - value per row. Filters rows where the column is missing. - - Params - ------ - df : pandas.DataFrame - dataframe with the column to split and expand - column : str - the column to split and expand - sep : str - the string used to split the column's values - keep : bool - whether to retain the presplit value as it's own row - - Returns - ------- - pandas.DataFrame - Returns a dataframe with the same columns as `df`. - """ - indexes = list() - new_values = list() - #df = df.dropna(subset=[column])#!!!!-----see this incase you need to uncomment based on use case - for i, presplit in enumerate(df[column].astype(str)): - values = presplit.split(sep) - if keep and len(values) > 1: - indexes.append(i) - new_values.append(presplit) - for value in values: - indexes.append(i) - new_values.append(value) - new_df = df.iloc[indexes, :].copy() - new_df[column] = new_values - return new_df - -#%% end of tidy_split() -#========= -# DF1: dr_mutations_pyrazinamide -#========= -######## -# tidy_split(): on 'dr_mutations_pyrazinamide' column and remove leading white spaces -######## -col_to_split1 = 'dr_mutations_pyrazinamide' -print ('Firstly, applying tidy split on dr df:', meta_pnca_dr.shape, - '\ncolumn name:', col_to_split1) -print('======================================================================') -# apply tidy_split() -dr_WF0 = tidy_split(meta_pnca_dr, col_to_split1, sep = ';') -# remove leading white space else these are counted as distinct mutations as well -dr_WF0['dr_mutations_pyrazinamide'] = dr_WF0['dr_mutations_pyrazinamide'].str.lstrip() - -# extract only the samples/rows with pncA_p. -dr_pnca_WF0 = dr_WF0.loc[dr_WF0.dr_mutations_pyrazinamide.str.contains(gene_match)] - -print('lengths after tidy split and extracting', gene_match, 'muts:', - '\nold length:' , len(meta_pnca_dr), - '\nlen after split:', len(dr_WF0), - '\ndr_pnca_WF0 length:', len(dr_pnca_WF0), - '\nexpected len:', dr_pnca_count) - -if len(dr_pnca_WF0) == dr_pnca_count: - print('PASS: length of dr_pnca_WF0 match with expected length') -else: - print('FAIL: lengths mismatch') - -print('======================================================================') - -# count the freq of "dr_muts" samples -dr_muts_df = dr_pnca_WF0 [['id', 'dr_mutations_pyrazinamide']] -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) -print('======================================================================') - -# sanity check: length of pnca samples -if len(dr_pnca_WF0) == c1.sum(): - print('PASS: WF data has expected length', - '\nlength of dr_pnca WFO:', c1.sum() ) -else: - print("FAIL: Debug please!") - -print('======================================================================') - -#!!! 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_pnca_WF0.assign(mutation_info = 'dr_mutations_pyrazinamide') -print("dim of dr_df:", dr_df.shape) -print('======================================================================') -print('End of tidy split() on dr_muts, and added an extra column relecting mut_category') -print('======================================================================') - -#%% -#========= -# DF2: other_mutations_pyrazinamdie -#========= -######## -# tidy_split(): on 'other_mutations_pyrazinamide' column and remove leading white spaces -######## -col_to_split2 = 'other_mutations_pyrazinamide' -print ("applying second tidy split separately on df:", meta_pnca_other.shape, '\n' - , "column name:", col_to_split2) -print('======================================================================') - -# apply tidy_split() -other_WF1 = tidy_split(meta_pnca_other, col_to_split2, sep = ';') -# remove the leading white spaces in the column -other_WF1['other_mutations_pyrazinamide'] = other_WF1['other_mutations_pyrazinamide'].str.strip() - -# extract only the samples/rows with pncA_p. -other_pnca_WF1 = other_WF1.loc[other_WF1.other_mutations_pyrazinamide.str.contains(gene_match)] - -print('lengths after tidy split and extracting', gene_match, 'muts:', - '\nold length:' , len(meta_pnca_other), - '\nlen after split:', len(other_WF1), - '\nother_pnca_WF1 length:', len(other_pnca_WF1), - '\nexpected len:', other_pnca_count) - -if len(other_pnca_WF1) == other_pnca_count: - print('PASS: length of dr_pnca_WF0 match with expected length') -else: - print('FAIL: lengths mismatch') - -print('======================================================================') -# count the freq of "other muts" samples -other_muts_df = other_pnca_WF1 [['id', 'other_mutations_pyrazinamide']] -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 pnca samples -if len(other_pnca_WF1) == c2.sum(): - print('PASS: WF data has expected length', - '\nlength of other_pnca WFO:', c2.sum() ) -else: - print("FAIL: Debug please!") -print('======================================================================') - -#!!! 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_pnca_WF1.assign(mutation_info = 'other_mutations_pyrazinamide') -print("dim of other_df:", other_df.shape) -print('======================================================================') -print('End of tidy split() on other_muts, and added an extra column relecting mut_category') -print('======================================================================') -#%% -#========== -# 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') - -dr_df.columns -dr_df.rename(columns = {'dr_mutations_pyrazinamide': 'mutation'}, inplace = True) -dr_df.columns - -other_df.columns -other_df.rename(columns = {'other_mutations_pyrazinamide': 'mutation'}, inplace = True) -other_df.columns - -print('======================================================================') -print('Now appending the two dfs:', - '\ndr_df dim:', dr_df.shape, - '\nother_df dim:', other_df.shape, - '\ndr_df length:', len(dr_df), - '\nother_df length:', len(other_df), - '\nexpected length:', len(dr_df) + len(other_df) ) -print('======================================================================') - -# 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: - print("FAIL: Debug please!") - -# concatenate (axis = 0): Rbind -pnca_LF0 = pd.concat([dr_df, other_df], ignore_index = True, axis = 0) - -# checking colnames and length after concat -print("checking colnames AFTER concatenating the two dfs...") -if (set(dr_df.columns) == set(pnca_LF0.columns)): - print('PASS: column names match') -else: - print("FAIL: Debug please!") - -print("checking length AFTER concatenating the two dfs...") - -if len(pnca_LF0) == len(dr_df) + len(other_df): - print("PASS:length of df after concat match") -else: - print("FAIL: length mismatch") -print('======================================================================') -#%% -########### -# This is hopefully clean data, but just double check -# Filter LF data so that you only have -# mutations corresponding to pncA_p.* (string match pattern) -# this will be your list you run OR calcs -########### -print('length of pnca_LF0:', len(pnca_LF0), - '\nThis should be what you need. But just double check and extract', gene_match, - '\nfrom LF0 (concatenated data)') - -print('using string match:', gene_match) - -print('Double checking and creating df: pnca_LF1') -pnca_LF1 = pnca_LF0[pnca_LF0['mutation'].str.contains(gene_match)] - -if len(pnca_LF0) == len(pnca_LF1): - print('PASS: length of pnca_LF0 and pnca_LF1 match', - '\nconfirming extraction and concatenating worked correctly') -else: - print('FAIL: BUT NOT FATAL!', - '\npnca_LF0 and pnca_LF1 lengths differ', - '\nsuggesting error in extraction process' - ' use pnca_LF1 for downstreama analysis') -print('======================================================================') -print('length of dfs pre and post processing...', - '\npnca WF data (unique samples in each row):', extracted_pnca_samples, - '\npnca LF data (unique mutation in each row):', len(pnca_LF1)) -print('======================================================================') -#%% -# final sanity check -print('Verifying whether extraction process worked correctly...') -if len(pnca_LF1) == expected_rows: - print('PASS: extraction process performed correctly', - '\nexpected length:', expected_rows, - '\ngot:', len(pnca_LF1), - '\nRESULT: Total no. of mutant sequences for logo plot:', len(pnca_LF1) ) - -else: - print('FAIL: extraction process has bugs', - '\nexpected length:', expected_rows, - '\ngot:', len(pnca_LF1), - '\Debug please') -#%% -print('Perfmorning some more sanity checks...') -# From LF1: -# no. of unique muts -distinct_muts = pnca_LF1.mutation.value_counts() -print("Distinct mutations:", len(distinct_muts)) - -# no. of samples contributing the unique muta -pnca_LF1.id.nunique() -print("No.of samples contributing to distinct muts:", pnca_LF1.id.nunique() ) - -# no. of dr and other muts -mut_grouped = pnca_LF1.groupby('mutation_info').mutation.nunique() -print("No.of unique dr and other muts:", pnca_LF1.groupby('mutation_info').mutation.nunique() ) - -# sanity check -if len(distinct_muts) == mut_grouped.sum() : - print("PASS:length of LF1 is as expected ") -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 ...') - muts_split = list(pnca_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 !!!! -# sanity check: There should not be any common muts -# i.e 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') -else: - print('PASS: NO ambiguous muts detected', - '\nMuts are unique to dr_ and other_ mutation class') - -# 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:\n', 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: - print('Error: ambiguous muts present, but extraction failed. Debug!') - -print('======================================================================') -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('No. of ambigiuous muts detected:'+ str(len(common_muts)), - 'list of ambiguous mutations (see below):', *common_muts, sep = '\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())) - -#%% clear variables -del(id_dr, id_other, meta_data, meta_pnca_dr, meta_pnca_other, mut_grouped, muts_split, other_WF1, other_df, other_muts_df, other_pnca_count, pnca_LF0, pnca_na) - -del(c1, c2, col_to_split1, col_to_split2, comp_pnca_samples, dr_WF0, dr_df, dr_muts_df, dr_pnca_WF0, dr_pnca_count, expected_pnca_samples, other_pnca_WF1) - -#%% end of data extraction and some files writing. Below are some more files writing. - -#%%: 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_filename1 = gene.lower() + '_' + 'ambiguous_muts.csv' -outfile1 = homedir + '/' + outdir + '/' + out_filename1 -print('Writing file: ambiguous muts...', - '\nFilename:', out_filename1, - '\nPath:', homedir +'/'+ outdir) - -#common_muts = ['pncA_p.Val180Phe','pncA_p.Gln10Pro'] # test -inspect = pnca_LF1[pnca_LF1['mutation'].isin(common_muts)] -inspect.to_csv(outfile1) - -print('Finished writing:', out_filename1, '\nExpected no. of rows (no. of samples with the ambiguous muts present):', dr_muts.isin(other_muts).sum() + other_muts.isin(dr_muts).sum()) -print('======================================================================') -del(out_filename1) - - -#%% -#=========== -# Split 'mutation' column into three: wild_type, position and -# mutant_type separately. Then map three letter code to one using -# reference_dict. -# First: Import reference dict -# Second: convert to mutation to lowercase for compatibility with dict -#=========== -from reference_dict import my_aa_dict # CHECK DIR STRUC THERE! -pnca_LF1['mutation'] = pnca_LF1.loc[:, 'mutation'].str.lower() -#======= -# 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 -#======= -lookup_dict = dict() -for k, v in my_aa_dict.items(): - lookup_dict[k] = v['one_letter_code'] - wt = pnca_LF1['mutation'].str.extract('pnca_p.(\w{3})').squeeze() # converts to a series that map works on - pnca_LF1['wild_type'] = wt.map(lookup_dict) - mut = pnca_LF1['mutation'].str.extract('\d+(\w{3})$').squeeze() - pnca_LF1['mutant_type'] = mut.map(lookup_dict) - -# extract position info from mutation column separetly using string match -pnca_LF1['position'] = pnca_LF1['mutation'].str.extract(r'(\d+)') - -# clear variables -del(k, v, wt, mut, lookup_dict) -print('======================================================================') -#========= -# 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. -#========= -# initialise a sub dict that is lookup dict for three letter code to aa prop -lookup_dict = dict() - -for k, v in my_aa_dict.items(): - lookup_dict[k] = v['aa_prop_water'] - #print(lookup_dict) - wt = pnca_LF1['mutation'].str.extract('pnca_p.(\w{3})').squeeze() # converts to a series that map works on - pnca_LF1['wt_prop_water'] = wt.map(lookup_dict) - mut = pnca_LF1['mutation'].str.extract('\d+(\w{3})$').squeeze() - pnca_LF1['mut_prop_water'] = mut.map(lookup_dict) - -# added two more cols - -# clear variables -del(k, v, wt, mut, lookup_dict) -print('======================================================================') -#======== -# 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. -#========= -# initialise a sub dict that is lookup dict for three letter code to aa prop -lookup_dict = dict() - -for k, v in my_aa_dict.items(): - lookup_dict[k] = v['aa_prop_polarity'] - #print(lookup_dict) - wt = pnca_LF1['mutation'].str.extract('pnca_p.(\w{3})').squeeze() # converts to a series that map works on - pnca_LF1['wt_prop_polarity'] = wt.map(lookup_dict) - mut = pnca_LF1['mutation'].str.extract(r'\d+(\w{3})$').squeeze() - pnca_LF1['mut_prop_polarity'] = mut.map(lookup_dict) - -# added two more cols - -# clear variables -del(k, v, wt, mut, lookup_dict) -print('======================================================================') - -#======== -# iterate through the dict, create a lookup dict that i.e -# lookup_dict = {three_letter_code: aa_taylor} -# Do this for both wild_type and mutant as above. -# caution: taylor mapping will create a list within a column -#========= -#lookup_dict = dict() - -#for k, v in my_aa_dict.items(): -# lookup_dict[k] = v['aa_taylor'] -# #print(lookup_dict) -# wt = pnca_LF1['mutation'].str.extract('pnca_p.(\w{3})').squeeze() # converts to a series that map works on -# pnca_LF1['wt_taylor'] = wt.map(lookup_dict) -# mut = pnca_LF1['mutation'].str.extract(r'\d+(\w{3})$').squeeze() -# pnca_LF1['mut_taylor'] = mut.map(lookup_dict) - -# added two more cols -# clear variables -#del(k, v, wt, mut, lookup_dict) -#print('======================================================================') - -#======== -# 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. -#========= -lookup_dict = dict() - -for k, v in my_aa_dict.items(): - lookup_dict[k] = v['aa_calcprop'] - #print(lookup_dict) - wt = pnca_LF1['mutation'].str.extract('pnca_p.(\w{3})').squeeze() # converts to a series that map works on - pnca_LF1['wt_calcprop'] = wt.map(lookup_dict) - mut = pnca_LF1['mutation'].str.extract(r'\d+(\w{3})$').squeeze() - pnca_LF1['mut_calcprop'] = mut.map(lookup_dict) - -# added two more cols -# clear variables -del(k, v, wt, mut, lookup_dict) -print('======================================================================') -######## -# 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 -######### -pnca_LF1['Mutationinformation'] = pnca_LF1['wild_type'] + pnca_LF1.position.map(str) + pnca_LF1['mutant_type'] -print('Created column: Mutationinformation') -print('======================================================================') -#%% Write file: mCSM muts -snps_only = pd.DataFrame(pnca_LF1['Mutationinformation'].unique()) -snps_only.head() -# assign column name -snps_only.columns = ['Mutationinformation'] - -# count how many positions this corresponds to -pos_only = pd.DataFrame(pnca_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') -else: - print('FAIL: SNP has NA, Possible mapping issues from dict?', - '\nDebug please!') -print('======================================================================') - -out_filename2 = gene.lower() + '_' + 'mcsm_snps.csv' -outfile2 = homedir + '/' + outdir + '/' + out_filename2 - -print('Writing file: mCSM style muts', - '\nmutation format (SNP): {Wt}{Mut}', - '\nNo. of distinct muts:', len(snps_only), - '\nNo. of distinct positions:', len(pos_only), - '\nFilename:', out_filename2, - '\nPath:', homedir +'/'+ outdir) - -snps_only.to_csv(outfile2, header = False, index = False) - -print('Finished writing:', out_filename2, - '\nNo. of rows:', len(snps_only) ) -print('======================================================================') -del(out_filename2) - - -#%% Write file: pnca_metadata (i.e pnca_LF1) -# where each row has UNIQUE mutations NOT unique sample ids -out_filename3 = gene.lower() + '_' + 'metadata.csv' -outfile3 = homedir + '/' + outdir + '/' + out_filename3 -print('Writing file: LF formatted data', - '\nFilename:', out_filename3, - '\nPath:', homedir +'/'+ outdir) - -pnca_LF1.to_csv(outfile3, header = True, index = False) -print('Finished writing:', out_filename3, - '\nNo. of rows:', len(pnca_LF1), - '\nNo. of cols:', len(pnca_LF1.columns) ) -print('======================================================================') -del(out_filename3) - - -#%% write file: mCSM style but with repitions for MSA and logo plots -all_muts_msa = pd.DataFrame(pnca_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') -else: - print('FAIL: SNP has NA, Possible mapping issues from dict?', - '\nDebug please!') -print('======================================================================') - -out_filename4 = gene.lower() + '_' + 'all_muts_msa.csv' -outfile4 = homedir + '/' + outdir + '/' + out_filename4 - -print('Writing file: mCSM style muts for msa', - '\nmutation format (SNP): {Wt}{Mut}', - '\nNo.of lines of msa:', len(all_muts_msa), - '\nFilename:', out_filename4, - '\nPath:', homedir +'/'+ outdir) - -all_muts_msa_sorted.to_csv(outfile4, header = False, index = False) - -print('Finished writing:', out_filename4, - '\nNo. of rows:', len(all_muts_msa) ) -print('======================================================================') -del(out_filename4) - - -#%% write file for mutational positions -# count how many positions this corresponds to -pos_only = pd.DataFrame(pnca_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_filename5 = gene.lower() + '_' + 'mutational_positons.csv' -outfile5 = homedir + '/' + outdir + '/' + out_filename5 - -print('Writing file: mutational positions', - '\nNo. of distinct positions:', len(pos_only_sorted), - '\nFilename:', out_filename5, - '\nPath:', homedir +'/'+ outdir) - -pos_only_sorted.to_csv(outfile5, header = True, index = False) - -print('Finished writing:', out_filename5, - '\nNo. of rows:', len(pos_only_sorted) ) -print('======================================================================') -del(out_filename5) - - -#%% end of script -print('======================================================================') -print(u'\u2698' * 50, - '\nEnd of script: Data extraction and writing files' - '\n' + u'\u2698' * 50 ) -#%%