import pandas as pd from sparrow import Protein import numpy as np import matplotlib.pyplot as plt import re import os import subprocess import warnings import sys from scipy.stats import ttest_rel warnings.filterwarnings("ignore", "You are using `torch.load` with `weights_only=False`*.") warnings.filterwarnings("ignore", "FutureWarning: The behavior of DataFrame concatenation with empty or all-NA entries is deprecated*.") def makeFig(df,disorder,upID): ''' Make a figure of disorder for a given protein ''' print('making figure for '+upID) fig,ax = plt.subplots(figsize=[7,3]) ax.plot(np.arange(len(disorder)),disorder) ax.scatter(df['pos'],np.ones(len(df))*0.5,s=0.5,c='r') ax.text(0.95,0.95,upID,ha='right',transform = ax.transAxes) ax.set_ylabel('disorder') ax.set_xlabel('residue') fig.savefig('./output/'+upID+'/map.png') plt.close(fig) def assignGroup(aa): ''' Helper function to assign chemistries to amino acids ''' groups={'polar':['Q','N','S','T','G','H'], 'apolar':['I','L','V','A','M'], 'aro':['W','F','Y'], 'neg':['D','E'], 'pos':['K','R'], 'cys':['C'], 'pro':['P']} result_keys = [key for key, value_list in groups.items() if aa in value_list][0] return(result_keys) def assignChange(change): ''' Helped function to assign the mutation change type ''' if bool(re.search(r'(\w+)>(\1)',change))==True: return('none') else: return(change) def get_windows_containing_position(varID, wtSeq, pos, mut, length=50): """ Get all tiles of a given length that contain the mutation mut. Predict the Re of each tile with and without the mutant. Return a dataframe with all the predictions. """ if length > len(wtSeq): return(pd.DataFrame()) wtRe = [] mutRe = [] mutSeq = wtSeq[:pos]+mut+wtSeq[pos+1:] print("wtSeq: %s \nmutSeq: %s \n orig: %s \n pos: %i \n mut: %s"%(wtSeq,mutSeq,wtSeq[pos],pos,mut)) # Calculate the range of possible window starts # For each window, we need: start_pos <= position < start_pos + length earliest_start = max(0, pos - length+1 ) latest_start = min(pos, len(wtSeq) - length) print('earliest: %i latest: %i'%(earliest_start,latest_start)) # Get all possible windows for start in range(earliest_start, latest_start + 1): wtWindow = wtSeq[start:start + length] mutWindow = mutSeq[start:start + length] wtP=Protein(wtWindow) mutP=Protein(mutWindow) wtRe.append(wtP.predictor.end_to_end_distance(use_scaled=True)) mutRe.append(mutP.predictor.end_to_end_distance(use_scaled=True)) print('seq: %s / Re=%.3f'%(mutWindow,mutRe[-1])) wtRe_local = np.array(wtRe) mutRe_local = np.array(mutRe) return(pd.DataFrame({'mut_Re_avg_'+str(length):(mutRe_local).mean(), 'wt_Re_avg_'+str(length):(wtRe_local).mean(), 'mut_vs_wt_Re_min_'+str(length):(mutRe_local/wtRe_local).min(), 'mut_vs_wt_Re_max_'+str(length):(mutRe_local/wtRe_local).max(), 'mut_vs_wt_Re_avg_'+str(length):(mutRe_local/wtRe_local).mean(), 'mut_vs_wt_Re_std_'+str(length):(mutRe_local/wtRe_local).std(), 'mut_vs_wt_Re_delta_'+str(length):(mutRe_local-wtRe_local).mean(), 'mut_vs_wt_Re_logFC_'+str(length):np.log2((mutRe_local/wtRe_local).mean()), 'mut_vs_wt_Re_pval_'+str(length):-np.log10(ttest_rel(mutRe_local,wtRe_local).pvalue), 'mut_vs_wt_Re_N_'+str(length):len(mutRe_local)},index=[varID])) def getVariant(varID): ''' Pull information from a variant ID Input file is ./data/clinvar_20240917.vcf returns disease, significance, and review status ''' results = subprocess.run("awk '{if ($3=="+str(varID)+") print $NF}' ./data/clinvar_20240917.vcf", shell=True, capture_output=True, text=True) info = results.stdout disease = re.search(r'CLNDN=(.*?);',info).group(1).replace(",", ";") signi = re.search(r'CLNSIG=(.*?);',info).group(1).replace(",", ";") revstat = re.search(r'CLNREVSTAT=(.*?);',info).group(1).replace(",", ";") return(disease,signi,revstat)#[disease, signi]) def getProteinClinvar(upID,geneID): ''' Main function. Accepts uniprot ID and geneID for a single protein returns all missense ClinVar mutations and predictions of missense mutations for this protein ''' thisProtVars=pd.DataFrame() outdir = './output/'+upID+"/" os.makedirs(outdir,exist_ok=True) os.makedirs(outdir+'fasta/',exist_ok=True) hgvs_out = outdir+'hgvs.tsv' uniprot_out = outdir+'uniprot.tsv' with open(hgvs_out,'w') as outfile: subprocess.run("head -n1 ./data/hgvs_missense.txt",shell=True, stdout=outfile) subprocess.run("cat ./data/hgvs_missense.txt |awk '{if ($3 == "+str(geneID)+") print $0}'",shell=True, stdout=outfile) with open(uniprot_out,'w') as outfile: subprocess.run("head -n1 ./data/uniprotkb_proteome_UP000005640_2024_08_16.tsv",shell=True, stdout=outfile) result = subprocess.run("grep "+str(upID)+" ./data/uniprotkb_proteome_UP000005640_2024_08_16.tsv",shell=True, stdout=outfile) sl_uniprot = pd.read_table(uniprot_out,header=0) sl_hgvs = pd.read_table(hgvs_out,header=0,index_col=None) print(len(sl_hgvs)) sl_hgvs = sl_hgvs.drop_duplicates('VariationID') print(len(sl_hgvs)) # now start pulling data upName = sl_uniprot['Entry Name'].values[0].split('_')[0] seq = sl_uniprot['Sequence'].values[0] N_res = len(seq) P = Protein(seq) disorder = np.array(P.predictor.disorder()) N_disordered = disorder[disorder>0.7].size N_ordered = disorder[disorder<0.3].size wtRe = P.predictor.end_to_end_distance(use_scaled=True) for hgvIdx in sl_hgvs.index: varID = sl_hgvs.loc[hgvIdx]['VariationID'] print("%s: on varID %s, %i / %i" % (upName, varID, hgvIdx, np.max(sl_hgvs.index))) change = sl_hgvs.loc[hgvIdx]['ProteinChange'] orig = re.search(r'^.*\.(.*?)[0-9].*$',change).group(1) pos = int(re.search(r'\-?\d+',change).group(0)) mut = re.search(r'^.*.[0-9](.*)$',change).group(1) if (mut not in list(amino_acid_map.keys()))|(orig not in list(amino_acid_map.keys())): print(orig+" or "+mut+" not an amino acid") continue if (pos > len(seq)): print(str(pos)+" beyond original length "+str(len(seq))) continue if (seq[pos-1]!=amino_acid_map[orig]): print("original aa not the same") continue (disease,signi,revstat)=getVariant(varID) Re_local_df=pd.DataFrame() for tile_length in [20,25,30,35,40,45,50,60]: print('tile length %i'%tile_length) df = get_windows_containing_position(varID, seq, pos, amino_acid_map[mut],tile_length) if df.empty: continue else: # df.columns=[a+"_"+b for a,b in zip(df.columns,[str(tile_length)]*len(df.columns))] Re_local_df = pd.concat([Re_local_df,df],axis=1) thisProtVars = pd.concat([thisProtVars, pd.concat([pd.DataFrame({ 'upID':upID, 'geneID':geneID, 'N_res':N_res, 'N_disordered':N_disordered, 'N_ordered':N_ordered, 'disease':disease, 'signi':signi, 'revstat':revstat, 'res_disorder':disorder.tolist()[pos-1], 'change':change, 'changeType':assignChange(assignGroup(amino_acid_map[orig])+">"+ assignGroup(amino_acid_map[mut])), 'orig':amino_acid_map[orig], 'pos':pos, 'mut':amino_acid_map[mut]},index=[varID]),Re_local_df],axis=1)]) thisProtVars.to_csv(outdir+'vars.csv') # thisProtVars.to_csv('tmp.csv') # makeFig(thisProtVars,disorder,upID) print('%s finished. total passed %i variants' % (upID, len(thisProtVars))) amino_acid_map = { 'Ala': 'A', 'Arg': 'R', 'Asn': 'N', 'Asp': 'D', 'Cys': 'C', 'Gln': 'Q', 'Glu': 'E', 'Gly': 'G', 'His': 'H', 'Ile': 'I', 'Leu': 'L', 'Lys': 'K', 'Met': 'M', 'Phe': 'F', 'Pro': 'P', 'Ser': 'S', 'Thr': 'T', 'Trp': 'W', 'Tyr': 'Y', 'Val': 'V' } # ''' Main function first user input is Uniprot ID second user input is gene ID ''' upID = sys.argv[1] geneID = sys.argv[2] getProteinClinvar(upID,geneID)