223 lines
7.9 KiB
Python
223 lines
7.9 KiB
Python
import requests, sys, json
|
|
from prettify import warn, green, red, blue
|
|
|
|
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'],
|
|
}
|
|
|
|
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_rhea_candidates', help="number of candidate rhea entries", type=int, default=1)
|
|
parser.add_argument('-v', '--verbose', action="store_true", default=False)
|
|
return parser.parse_args(sys.argv[1:])
|
|
|
|
def read_tsv(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 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 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 main():
|
|
args = parse_args()
|
|
rhea_dict = read_tsv("./rhea_unclassified_unique_ec4_candidates.tsv", top=args.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"]<=args.num_uniprot_entries])
|
|
print(f"\nFound {blue(num_candidates)} rhea entries with <={args.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"]<=args.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 args.verbose:
|
|
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(f"\nFound a total of {warn(num_hits)} matching active sites.")
|
|
|
|
def print_AA_differences(mcsa_aas, sog_aas, uniprot_aas):
|
|
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("")
|
|
|
|
if __name__=="__main__":
|
|
main() |