390 lines
15 KiB
Python
390 lines
15 KiB
Python
import requests, sys, json
|
|
from prettify import warn, green, red, blue
|
|
from grantham_matrix import min_distance_LP
|
|
|
|
from argparse import ArgumentParser
|
|
|
|
AA_SINGLE = {
|
|
"A": ['ALANINE', 'Ala'],
|
|
"R": ['ARGININE', 'Arg'],
|
|
"N": ['ASPARAGINE', 'Asn'],
|
|
"D": ['ASPARTIC ACID', 'Asp'],
|
|
"C": ['CYSTEINE', 'Cys'],
|
|
"E": ['GLUTAMIC ACID', 'Glu'],
|
|
"Q": ['GLUTAMINE', 'Gln'],
|
|
"G": ['GLYCINE', 'Gly'],
|
|
"H": ['HISTIDINE', 'His'],
|
|
"I": ['ISOLEUCINE', 'Ile'],
|
|
"L": ['LEUCINE', 'Leu'],
|
|
"K": ['LYSINE', 'Lys'],
|
|
"M": ['METHIONINE', 'Met'],
|
|
"F": ['PHENYLALANINE', 'Phe'],
|
|
"P": ['PROLINE', 'Pro'],
|
|
"S": ['SERINE', 'Ser'],
|
|
"T": ['THREONINE', 'Thr'],
|
|
"W": ['TRYPTOPHAN', 'Trp'],
|
|
"Y": ['TYROSINE', 'Tys'],
|
|
"V": ['VALINE', 'Val'],
|
|
"U": ['SELENOCYSTEINE', 'Sec']
|
|
}
|
|
|
|
AA_TRIPLET = {t : s for (s, [_,t]) in list(AA_SINGLE.items())}
|
|
AA_TRIPLET['Tyr'] = "Y"
|
|
AA_TRIPLET['Sec'] = "U"
|
|
|
|
def parse_args():
|
|
parser = ArgumentParser(
|
|
description='Compare AAs in overlay graphs with UniProt')
|
|
parser.add_argument('-p', '--num_uniprot_entries', help="number of uniprot entries per rhea entry", type=int, default=1)
|
|
parser.add_argument('-t', '--num_candidates', help="number of candidate rhea entries", type=int, default=5)
|
|
parser.add_argument('-md', '--max_distance', type=int, default=500)
|
|
parser.add_argument('-v', '--verbose', action="store_true", default=False)
|
|
return parser.parse_args(sys.argv[1:])
|
|
|
|
def read_tsv_rhea(filepath, top=50):
|
|
with open(filepath) as f:
|
|
glob = f.read()
|
|
|
|
# rhea dict :: rhea_master_id -> attrs
|
|
rhea_dict = {}
|
|
# skip header line
|
|
for line in glob.split("\n")[1:top]:
|
|
values = line.split("\t")
|
|
rhea_master_id = int(values[1])
|
|
rhea_sub_ids = values[2].split("; ")
|
|
predicted_ec4s = values[5]
|
|
top_mcsa_ids = values[14].split("; ")
|
|
expected_amino_acids_overlay = values[19].split("; ")
|
|
expected_amino_acids_all_mcsa_entry = values[20].split("; ")
|
|
|
|
rhea_dict[rhea_master_id] = {
|
|
"sub_ids": rhea_sub_ids,
|
|
"predicted_ec4s": predicted_ec4s,
|
|
"top_mcsa_ids": top_mcsa_ids,
|
|
"expected_amino_acids_overlay" : expected_amino_acids_overlay,
|
|
"expected_amino_acids_all_mcsa_entry" : expected_amino_acids_all_mcsa_entry
|
|
}
|
|
return rhea_dict
|
|
|
|
def read_tsv_mcsa(filepath, top=50):
|
|
with open(filepath) as f:
|
|
glob = f.read()
|
|
|
|
# mcsa dict :: msca id -> attrs
|
|
mcsa_dict = {}
|
|
# skip header line
|
|
for line in glob.split("\n")[1:top]:
|
|
values = line.split("\t")
|
|
rhea_master_id = int(values[1])
|
|
rhea_sub_ids = values[2].split("; ")
|
|
predicted_ec4s = values[5]
|
|
top_mcsa_ids = values[14].split("; ")
|
|
expected_amino_acids_overlay = values[19].split("; ")
|
|
expected_amino_acids_all_mcsa_entry = values[20].split("; ")
|
|
|
|
for mcsa_id in top_mcsa_ids:
|
|
mcsa_dict[mcsa_id] = {
|
|
'rhea_master_id' : rhea_master_id,
|
|
'rhea_sub_ids' : rhea_sub_ids,
|
|
'predicted_ec4s' : predicted_ec4s,
|
|
'expected_amino_acids_overlay' : expected_amino_acids_overlay,
|
|
'expected_amino_acids_all_mcsa_entry' : expected_amino_acids_all_mcsa_entry
|
|
}
|
|
|
|
return mcsa_dict
|
|
|
|
def get_req_mcsa(mcsa_dict):
|
|
if not mcsa_dict:
|
|
print(red("No MCSA IDs to look up. (mcsa_dict is empty)"))
|
|
return mcsa_dict
|
|
url=f"https://www.ebi.ac.uk/thornton-srv/m-csa/api/entries/?format=json"
|
|
|
|
mcsa_ids = ",".join(sorted([x for x in list(mcsa_dict) if x]))
|
|
url_ids = url+f"&entries.mcsa_ids={mcsa_ids}"
|
|
print(url_ids)
|
|
|
|
response = requests.get(url_ids)
|
|
if not response.ok:
|
|
print(red(f"\tResponse not ok for MCSA."))
|
|
return mcsa_dict
|
|
|
|
data = json.loads(response.text)
|
|
print("number of results for get req mcsa:", len(data['results']))
|
|
return data['results']
|
|
|
|
def get_req_rhea(rheaID, rhea_dict):
|
|
url= f"https://www.rhea-db.org/rhea/?"
|
|
parameter = {
|
|
"query":f'{rheaID}',
|
|
"columns":"uniprot",
|
|
"format":'tsv',
|
|
"limit":1,
|
|
}
|
|
|
|
response = requests.get(url, params=parameter)
|
|
if not response.ok:
|
|
print(red(f"\tResponse not ok for {rheaID}. Skipping."))
|
|
return rhea_dict
|
|
|
|
data = response.text
|
|
if rheaID not in rhea_dict:
|
|
rhea_dict[rheaID] = {}
|
|
|
|
line = data.split("\n")[1]
|
|
if line:
|
|
num_uniprot = int(line)
|
|
else:
|
|
num_uniprot = 0
|
|
|
|
rhea_dict[rheaID]["num_uniprot"] = num_uniprot
|
|
return rhea_dict
|
|
|
|
def get_req_uniprot(rheaID):
|
|
print(f"Sending GET request to UniProt for rhea-id {blue(rheaID)}...")
|
|
base_url= "https://rest.uniprot.org/uniprotkb/stream?"
|
|
|
|
parameter = {
|
|
"query":f'((cc_catalytic_activity:"rhea:{rheaID}") AND (fragment:false) AND (reviewed:true))',
|
|
"format":'json',
|
|
"fields": [
|
|
"ft_binding",
|
|
"ft_act_site",
|
|
"sequence"
|
|
]
|
|
}
|
|
headers = {
|
|
"accept": "application/json"
|
|
}
|
|
response = requests.get(base_url, headers=headers, params=parameter)
|
|
|
|
if not response.ok:
|
|
response.raise_for_status()
|
|
|
|
data = response.json()
|
|
#print(data)
|
|
return data
|
|
|
|
def get_req_uniprotid(uniprotID):
|
|
print(f"Sending GET request to UniProt for UniProt ID {blue(uniprotID)}...")
|
|
params = {
|
|
"fields": [
|
|
"ft_binding",
|
|
"ft_act_site",
|
|
"ec",
|
|
"sequence"
|
|
]
|
|
}
|
|
headers = {
|
|
"accept": "application/json"
|
|
}
|
|
base_url = f"https://rest.uniprot.org/uniprotkb/{uniprotID}"
|
|
|
|
response = requests.get(base_url, headers=headers, params=params)
|
|
if not response.ok:
|
|
response.raise_for_status()
|
|
sys.exit()
|
|
|
|
data = response.json()
|
|
return data
|
|
|
|
def identify_aas_in_actsite(active_sites, sequence):
|
|
aas = []
|
|
for site in active_sites:
|
|
start = int(site['location']['start']['value'])-1
|
|
start_exact = site['location']['start']['modifier'] == "EXACT"
|
|
end = int(site['location']['end']['value'])-1
|
|
end_exact = site['location']['end']['modifier'] == "EXACT"
|
|
if not start_exact or not end_exact:
|
|
print(red("\tActive site locations not exact. Returning False."))
|
|
return False
|
|
|
|
for i in range(start,end+1):
|
|
aa_multi_code = AA_SINGLE[sequence[i]][1]
|
|
aas.append(f"{aa_multi_code}{i}")
|
|
|
|
return sorted(aas)
|
|
|
|
def parse_entries_by_rhea(num_rhea_candidates, num_uniprot_entries, verb=False):
|
|
rhea_dict = read_tsv_rhea("./rhea_unclassified_unique_ec4_candidates.tsv", top=num_rhea_candidates)
|
|
for rhea_master_id in rhea_dict:
|
|
get_req_rhea(rhea_master_id, rhea_dict)
|
|
|
|
num_candidates = sum([1 for r in rhea_dict if rhea_dict[r]["num_uniprot"]<=num_uniprot_entries])
|
|
print(f"\nFound {blue(num_candidates)} rhea entries with <={num_uniprot_entries} protein match(es) in UniProt.\n")
|
|
|
|
num_hits = 0
|
|
for rhea_master_id in rhea_dict:
|
|
if rhea_dict[rhea_master_id]["num_uniprot"]<=num_uniprot_entries:
|
|
data = get_req_uniprot(rhea_master_id)
|
|
result_match = False
|
|
for result in data['results']:
|
|
if result['features']:
|
|
# has active site info
|
|
sites = result['features']
|
|
sequence = result['sequence']['value']
|
|
aas = identify_aas_in_actsite(sites, sequence)
|
|
mcsa_aas = rhea_dict[rhea_master_id]['expected_amino_acids_all_mcsa_entry']
|
|
mcsa_aas_short = list(set([aa[:3] for aa in mcsa_aas]))
|
|
sog_aas = rhea_dict[rhea_master_id]['expected_amino_acids_overlay']
|
|
sog_aas_short = list(set([aa[:3] for aa in sog_aas]))
|
|
if aas == mcsa_aas:
|
|
print(green(f"\tMCSA hit."))
|
|
result_match = True
|
|
num_hits+=1
|
|
elif aas == sog_aas:
|
|
print(green(f"\tsOG hit."))
|
|
result_match = True
|
|
num_hits+=1
|
|
if verb:
|
|
print(f"\tMCSA ID: {rhea_dict[rhea_master_id]['top_mcsa_ids']}, UniProt ID: {result['primaryAccession']} ({result['entryType']})")
|
|
print_AA_differences(mcsa_aas, sog_aas, aas)
|
|
if not result_match:
|
|
#print(red(f"\tNo match found."))
|
|
print("")
|
|
|
|
print(f"\nFound a total of {warn(num_hits)} matching active sites.")
|
|
|
|
def parse_entries_by_mcsa(num_candidates, max_distance, verb=False):
|
|
#mcsa_dict = read_tsv_mcsa("./rhea_unclassified_unique_ec4_candidates.tsv", num_candidates)
|
|
mcsa_dict = [str(x) for x in list(range(1,num_candidates+1))]
|
|
print(f"Number of mcsa IDs: {len(mcsa_dict)}")
|
|
if len(mcsa_dict) < 100:
|
|
mcsaid_groups = [list(mcsa_dict)]
|
|
else: mcsaid_groups = list(zip(*(iter(mcsa_dict),) * 100))
|
|
distances = []
|
|
num_hits = 0
|
|
hits = []
|
|
for mcsa_dict in mcsaid_groups:
|
|
results = get_req_mcsa(mcsa_dict)
|
|
|
|
for record in results:
|
|
mcsa_id = record['mcsa_id']
|
|
residues = record['residues']
|
|
all_mcsa_ecs = record['all_ecs']
|
|
mcsa_aas = []
|
|
for residue in residues:
|
|
sequences = residue['residue_sequences']
|
|
for sequence in sequences:
|
|
aa_code = sequence['code']
|
|
aa_index = sequence['resid']
|
|
mcsa_aas.append(f"{aa_code}{aa_index}")
|
|
|
|
uniprot_ids = record['reference_uniprot_id']
|
|
for uniprot_id in uniprot_ids.split(", "):
|
|
data = get_req_uniprotid(uniprot_id)
|
|
result_match = False
|
|
if 'features' in data:
|
|
# has active site info
|
|
sites = data['features']
|
|
sequence = data['sequence']['value']
|
|
all_uniprot_ecs = [x['value'] for x in data['proteinDescription']['recommendedName']['ecNumbers']]
|
|
aas = identify_aas_in_actsite(sites, sequence)
|
|
#mcsa_aas = mcsa_dict[str(mcsa_id)]['expected_amino_acids_all_mcsa_entry']
|
|
mcsa_aas_short = list(set([aa[:3] for aa in mcsa_aas]))
|
|
if verb:
|
|
print(f"\tMCSA ID: {warn(mcsa_id)}, UniProt ID: {uniprot_id} ({data['entryType']})")
|
|
(distance, pairs, unpaired) = get_AA_distances(mcsa_aas, [], aas)
|
|
if distance <=max_distance:
|
|
hits.append((f"{mcsa_id}\t{uniprot_id}\t{mcsa_aas}\t{aas}\t{pairs}\t{unpaired}\t{distance}",distance))
|
|
if distance==0:
|
|
num_hits+=1
|
|
print_AA_differences(mcsa_aas, [], aas)
|
|
print_AA_distances(mcsa_aas, [], aas)
|
|
distances.append((mcsa_id, uniprot_id, distance))
|
|
print(red(f""))
|
|
|
|
print(f"\nFound a total of {warn(num_hits)} matching active sites.")
|
|
print(f"Num distances under {max_distance}:", blue(len([1 for (_,_,x) in distances if x<=max_distance])))
|
|
print(f"Smallest distance:", green(min(distances, key=lambda x: x[2])))
|
|
print(f"Largest distance:", red(max(distances, key=lambda x: x[2])))
|
|
print(f"Average distance:", warn(sum([x for (_,_,x) in distances])/len(distances)))
|
|
|
|
# Print MCSA ID and UniProt ID for distances <= max_distance
|
|
hits = sorted(hits, key=lambda x: x[1])
|
|
hits = [x for (x,_) in hits]
|
|
hits.insert(0,"MCSA ID\tUniProt ID\tMCSA AAs\tUniProt AAs\tPairs\tUnpaired\tDistance")
|
|
blob = "\n".join(hits)
|
|
with open(f"./distance_under_{max_distance}_from_{num_candidates}_MCSA_entries.tsv", 'w') as f:
|
|
f.write(blob)
|
|
|
|
def main():
|
|
args = parse_args()
|
|
#parse_entries_by_rhea(args.num_candidates, args.num_uniprot_entries, args.verbose)
|
|
parse_entries_by_mcsa(args.num_candidates, args.max_distance, verb=True)
|
|
|
|
|
|
def print_AA_differences(mcsa_aas, sog_aas, uniprot_aas):
|
|
# set to remove duplications, eg. His93 His93
|
|
mcsa_aas2 = sorted(list(set(mcsa_aas)))
|
|
sog_aas2 = sorted(list(set(sog_aas)))
|
|
uni_aas2 = sorted(list(set(uniprot_aas)))
|
|
|
|
mcsa_short = list(set([aa[:3] for aa in mcsa_aas2]))
|
|
sog_short = list(set([aa[:3] for aa in sog_aas2]))
|
|
uni_short = list(set([aa[:3] for aa in uni_aas2]))
|
|
|
|
mcsa_counts = {b : sum([1 for a in mcsa_aas2 if a[:3]==b]) for b in mcsa_short}
|
|
sog_counts = {b : sum([1 for a in sog_aas2 if a[:3]==b]) for b in sog_short}
|
|
uni_counts = {b : sum([1 for a in uni_aas2 if a[:3]==b]) for b in uni_short}
|
|
|
|
for shortcode in mcsa_short+sog_short+uni_short:
|
|
if shortcode not in mcsa_short:
|
|
mcsa_counts[shortcode] = 0
|
|
if shortcode not in sog_short:
|
|
sog_counts[shortcode] = 0
|
|
if shortcode not in uni_short:
|
|
uni_counts[shortcode] = 0
|
|
|
|
mcsa_f = {b : uni_counts[b] for b in mcsa_counts} # the amount that should be green
|
|
sog_f = {b : uni_counts[b] for b in mcsa_counts}
|
|
uni_f = {b : mcsa_counts[b] for b in mcsa_counts}
|
|
|
|
pretty_mcsa = []
|
|
pretty_sog = []
|
|
pretty_uni = []
|
|
|
|
for a in mcsa_aas2:
|
|
sc = a[:3]
|
|
pretty_mcsa.append(green(a) if mcsa_f[sc]>0 else red(a))
|
|
mcsa_f[sc]-=1
|
|
for a in sog_aas2:
|
|
sc = a[:3]
|
|
pretty_sog.append(green(a) if sog_f[sc]>0 else red(a))
|
|
sog_f[sc]-=1
|
|
for a in uni_aas2:
|
|
sc = a[:3]
|
|
pretty_uni.append(green(a) if uni_f[sc]>0 else red(a))
|
|
uni_f[sc]-=1
|
|
|
|
|
|
print("MCSA: ", end="")
|
|
[print(f"\t{x}", end="") for x in pretty_mcsa]
|
|
print("\nsOG: ", end="")
|
|
[print(f"\t{x}", end="") for x in pretty_sog]
|
|
print("\nUniP: ", end="")
|
|
[print(f"\t{x}", end="") for x in pretty_uni]
|
|
print("")
|
|
|
|
def get_AA_distances(mcsa_aas, sog_aas, uniprot_aas):
|
|
mcsa_short = list([aa[:3] for aa in list(set(mcsa_aas))])
|
|
uni_short = list([aa[:3] for aa in list(set(uniprot_aas))])
|
|
aa_set1 = [AA_TRIPLET[t] for t in mcsa_short]
|
|
aa_set2 = [AA_TRIPLET[t] for t in uni_short]
|
|
(objVal, pairs, unpaired) = min_distance_LP(aa_set1, aa_set2)
|
|
pairs = [(AA_SINGLE[a[0]][1], AA_SINGLE[b[0]][1]) for a,b in pairs]
|
|
unpaired = [AA_SINGLE[a[0]][1] for a in unpaired]
|
|
return (objVal, pairs, unpaired)
|
|
|
|
def print_AA_distances(mcsa_aas, sog_aas, uniprot_aas):
|
|
# print distance measure between MCSA and UniProt
|
|
(objVal, pairs, unpaired) = get_AA_distances(mcsa_aas, sog_aas, uniprot_aas)
|
|
pairs = [(AA_SINGLE[a[0]][1], AA_SINGLE[b[0]][1]) for a,b in pairs]
|
|
unpaired = [AA_SINGLE[a[0]][1] for a in unpaired]
|
|
print(f"Distance measure between MCSA and UniProt: {blue(objVal)}")
|
|
print(f"Pairs")
|
|
print(f"\t{'\t'.join([green(x) for (x,_) in pairs])}")
|
|
print(f"\t{'\t'.join([green(x) for (_,x) in pairs])}")
|
|
print(f"Unpaired: {blue(unpaired)}")
|
|
|
|
if __name__=="__main__":
|
|
main() |