<jupyter_start><jupyter_text>Mapping to Brenda to UniprotFirst I grabbed the "brenda_download.txt" file from the Brenda website (need to accept their user license)First, I used awk to extact the Uniprot accessions from this file and map them to Brenda's e.c. numbers. The regular expressions was based on that descibed here: https://www.uniprot.org/help/accession_numbers<jupyter_code>%%bash
###extract uniprot accessions out of brenda_download.txt
awk -vOFS="\t" '$0=="///"{current_ec="unknown"};$1=="ID"{current_ec=$2};patsplit($0,uniprot_accessions,/[OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2}/){for(i in uniprot_accessions)print(uniprot_accessions[i],current_ec)};' brenda_download.txt > brenda_uniprot_to_ec.tsv<jupyter_output><empty_output><jupyter_text>I took the Uniprot accessions from the first column of "brenda_uniprot_to_ec.tsv", and used the "retrieve ID" function on the Uniprot website, made sure the Pfam column was enabled, and exported it into "brenda_uniprot.tsv".Then we merge the data in that file with "brenda_uniprot_to_ec.tsv" as shown below to generate "uniprot_to_brenda_merged.tsv":<jupyter_code>import pandas
brenda_from_uniprot_df=pandas.read_table("brenda_uniprot.tsv",dtype="str")
brenda_from_uniprot_df.drop(columns=brenda_from_uniprot_df.columns[0],inplace=True)
brenda_from_uniprot_df
brenda_to_ec_df=pandas.read_table("brenda_uniprot_to_ec.tsv",dtype="str",names=["uniprot_accession","brenda_ec_number"])
brenda_to_ec_deduplicated=brenda_to_ec_df.groupby("uniprot_accession")['brenda_ec_number'].apply(lambda x:";".join(frozenset(x)))
brenda_uniprot_merged=pandas.merge(brenda_from_uniprot_df,brenda_to_ec_deduplicated,left_on=["Entry"],right_on=["uniprot_accession"])
brenda_uniprot_merged.set_index("Entry",inplace=True)
brenda_uniprot_merged.to_csv("uniprot_to_brenda_merged.tsv",sep="\t")<jupyter_output><empty_output><jupyter_text>Mapping MIBiG to NCBI protein accessions, and those in turn to UniprotDownloaded the MIBiG json files from the mibig website and extracted them into a folder. These json files do not give uniprot accessions, but only contain NCBI nucleotide accessions, and sometimes start and stop locations. So, we need to map these NCBI nucleotide accessions to NCBI protein accessions (Using NCBI's efetch tool), and then map NCBI protein accessions to Uniprot (Using Uniprot's ID mapping service)You ought to fill in your own e-mail address and api-key for NCBI:<jupyter_code>email_for_ncbi="<EMAIL>"
ncbi_api_key="get_this_from_your_NCBI_account"
#Run this cell only once

#The NCBI efetch service will sometimes time out (maybe when it's too busy?) which stops the loop in the next cell.
#Using this cache, we do not need to request data we have already requested before.
cache={}#map urls to response text
#If this cell fails, you may want to try running it again and see if it gets further along.

###this cell queries NCBI for all the protein IDs in Mibig

import requests
import pandas
import glob
import json
import re
import sys

session = requests.Session()
session.headers.update({'User-Agent': 'programmatic access by '+email_for_ncbi,
                       'From':email_for_ncbi,
                       })
ncbi_efetch_base_url = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?api_key="+ncbi_api_key+"&email="+email_for_ncbi+"&db=nuccore&rettype=ft&conwithfeat=on"
#the &conwithfeat=on part is not documented by NCBI, but it's required to make certain large records download that otherwise would just be empty 

mibig_file_list=glob.glob("mibig_json_2.0/*.json")

results_df=pandas.DataFrame(columns=["ncbi_protein_accession","mibig_accession","protein_id_type","mibig_taxid"])
results_df.set_index(["ncbi_protein_accession","mibig_accession"],inplace=True)

gb_regex=re.compile("protein_id.+?(gb|dbj|ref|gp|emb|tpe|tpg|tpd)\|(.{7,13}?)\..{1,2}\|")
for filename in mibig_file_list:
#     if(int(filename[-10:-5])<500):   #useful if you need to start somewhere in the middle
#          continue
    with open(filename) as f:
        mibig_file_dict=json.load(f)
        print("retrieving ncbi protein IDs for {}".format(mibig_file_dict["cluster"]["mibig_accession"]))
        mibig_locus_dict=mibig_file_dict["cluster"]["loci"]
        mibig_accession=mibig_file_dict["cluster"]["mibig_accession"]
        if "start_coord" in mibig_locus_dict:
            request_url=ncbi_efetch_base_url+"&id={}&seq_start={}&seq_stop={}".format(mibig_locus_dict["accession"],
                                                                                      mibig_locus_dict["start_coord"],
                                                                                      mibig_locus_dict["end_coord"]
                                                                                     )
        else:
            request_url=ncbi_efetch_base_url+"&id="+mibig_locus_dict["accession"]
            
        print("requesting "+request_url)
        if request_url in cache:
            print("Grabbing from cache")
            response_text=cache[request_url]        
        else:
            response_text=None
            r=session.get(request_url)
            if(r.ok):
                print("response in "+str(r.elapsed.total_seconds())+" s")
                if(r.text):
                    response_text=r.text
                    cache[request_url]=response_text
                else:
                    print("no results")
            else:
                print("error: "+str(r.status_code))
            time.sleep(5)
        if(response_text):
            total_proteins=0
            for regex_match in gb_regex.finditer(response_text):
                total_proteins+=1
                protein_id_type=regex_match.group(1)
                protein_accession=regex_match.group(2)
                if (protein_accession,mibig_accession) in results_df.index:
                    print("duplicates!!")
                else:  
                    results_df.loc[(protein_accession,mibig_accession),
                                   ["protein_id_type","mibig_taxid"]
                                   ]=[protein_id_type,mibig_file_dict["cluster"]["ncbi_tax_id"]]
                                  
            print("{}: found {} protein accessions".format(mibig_accession,total_proteins))
            print("table now has {} entries".format(len(results_df.index)))

results_df.to_csv("ncbi_accessions_to_mibig.tsv",sep="\t")<jupyter_output><empty_output><jupyter_text>we will use the file that this created, "ncbi_accessions_to_mibig.tsv", later on, and will use the dataframe in-memory in this notebook<jupyter_code>#generate 2 files for NCBI protein accessions. One for refseq accessions and one for genbank/EMB/DDBJ accessions

only_ref_df=results_df[results_df["protein_id_type"]=="ref"].reset_index().set_index("ncbi_protein_accession")
print("length of submit_to_uniprot_refseq_protein.txt: {} entries".format(len(only_ref_df)))
only_ref_df.to_csv("submit_to_uniprot_refseq_protein.txt",columns=[],header=False)

not_ref_df=results_df[results_df["protein_id_type"]!="ref"].reset_index().set_index("ncbi_protein_accession")
print("length of submit_to_uniprot_embl_genbank_cds.txt: {} entries".format(len(not_ref_df)))
not_ref_df.to_csv("submit_to_uniprot_embl_genbank_cds.txt",columns=[],header=False)<jupyter_output><empty_output><jupyter_text>We submit these two files to the UniProt ID mapping service:- Use the Refseq protein input option for "submit_to_uniprot_refseq_protein.txt". The results of this were saved in "mibig_uniprot_refseq_proteins.tsv"- And use the EMBL/Genbank CDS input option for "submit_to_uniprot_embl_genbank_cds.txt". For this file there are too many accession numbers to all submit at once, so I submitted in three batches and concatenated the results (saved in "mibig_uniprot_embl_genbank_cds_proteins_partA.tsv", etc., as follows:<jupyter_code>%%bash
#combine the three separate requests to uniprot
cat mibig_uniprot_embl_genbank_cds_proteins_partA.tsv > mibig_uniprot_embl_genbank_cds_proteins_combined.tsv
tail -n+2 mibig_uniprot_embl_genbank_cds_proteins_partB.tsv >> mibig_uniprot_embl_genbank_cds_proteins_combined.tsv
tail -n+2 mibig_uniprot_embl_genbank_cds_proteins_partC.tsv >> mibig_uniprot_embl_genbank_cds_proteins_combined.tsv
tail -n+2 mibig_uniprot_embl_genbank_cds_proteins_partD.tsv >> mibig_uniprot_embl_genbank_cds_proteins_combined.tsv<jupyter_output><empty_output><jupyter_text>Lastly we combine the metadata from the two different Uniprot ID mapping types, and merge that data with the MIBiG to NCBI protein accession mapping, to finally give a mapping from MIBiG BGC to Uniprot accession in the "uniprot_to_mibig_merged_on_protein_id.tsv" file.<jupyter_code>print("Table matching protein ids to mibig has {} entries".format(len(results_df)))

uniprot_refseq_df=pandas.read_table("mibig_uniprot_refseq_proteins.tsv",dtype="str")
uniprot_refseq_df.rename(columns={uniprot_refseq_df.columns[1]:"protein_id"},inplace=True)
uniprot_refseq_df.set_index(["Entry","protein_id"],inplace=True)
print("Refseq table has {} entries".format(len(uniprot_refseq_df)))

uniprot_embl_genbank_df=pandas.read_table("mibig_uniprot_embl_genbank_cds_proteins_combined.tsv",dtype="str")
uniprot_embl_genbank_df.rename(columns={uniprot_embl_genbank_df.columns[0]:"protein_id"},inplace=True)
uniprot_embl_genbank_df.drop(columns=uniprot_embl_genbank_df.columns[2],inplace=True)
uniprot_embl_genbank_df.set_index(["Entry","protein_id"],inplace=True)
uniprot_embl_genbank_df.drop_duplicates(inplace=True)
print("EMBL/Genbank CDS table has {} entries".format(len(uniprot_embl_genbank_df)))

uniprot_df=pandas.concat([uniprot_refseq_df,uniprot_embl_genbank_df],verify_integrity=True)
print("Combined uniprot table has {} entries".format(len(uniprot_df)))

#Merging on only protein_id. Useful for quickly checking if a uniprot ID is in mibig
merged_df_on_protein_id=pandas.merge(uniprot_df.reset_index(),results_df.reset_index(),left_on=["protein_id"],right_on=["ncbi_protein_accession"])
merged_df_on_protein_id.set_index("Entry",inplace=True)
merged_df_on_protein_id.to_csv("uniprot_to_mibig_merged_on_protein_id.tsv",sep="\t")
print("Merging on protein ID only leaves {} entries".format(len(merged_df_on_protein_id)))


#Merging on both protein_id and taxid.
#Unfortunately we do lose some entries where NONE of the uniprot results match the taxid in mibig
merged_df_on_protein_id_and_taxid=pandas.merge(uniprot_df.reset_index(),results_df.reset_index(),left_on=["protein_id","Taxonomic lineage IDs"],right_on=["ncbi_protein_accession","mibig_taxid"])
merged_df_on_protein_id_and_taxid.set_index("Entry",inplace=True)
merged_df_on_protein_id_and_taxid.to_csv("uniprot_to_mibig_merged_on_protein_id_and_taxid.tsv",sep="\t")
print("Merging on protein ID and taxid leaves {} entries".format(len(merged_df_on_protein_id_and_taxid)))<jupyter_output><empty_output><jupyter_text>Uniprot entries curated to have catalytic activitySearch uniprot for 'annotation:(type:"catalytic activity" evidence:manual)', made sure the Pfam column was turned on, and exported it as "Uniprot_curated_catalytic_activity.tsv".Uniprot also has an older index of proteins involved in metabolic pathways, called pathway.txt: https://www.uniprot.org/docs/pathwayI extracted the uniprot accessions out of that file using:<jupyter_code>%%bash
egrep -o '(......)' pathway.txt | awk '{print substr($0,2,6)}'>pathway_accessions.txt
#This is not ideal and will miss a couple of entries; in the future I should use the same regular expression as in the "Mapping to Brenda to Uniprot" section<jupyter_output><empty_output><jupyter_text>then I queried the uniprot KB using those accessions in 3 batches of less than 50,000, and concatenated the files to generate "uniprot_pathway_complete.tsv"  Pre-process all pfam annotations in UniprotFor this we actually use the protein2ipr.dat.gz file from the Interpro FTP server, which is conveniently sorted by uniprot accession, speeding up downstream processing<jupyter_code>##parsing protein2ipr.dat.gz file. uniprot IDs are already sorted alphabetically in this file
#the output files require about 2 Gb

protein2ipr_filename="D:/domain_analysis/protein2ipr.dat.gz"
output_folder="D:/domain_analysis/pfam_to_uniprot_files/"




import gzip
import collections

def write_progress_and_clear_lists(list_dict):
    print("writing progress...")
    for k,v in list_dict.items():
        with gzip.open(output_folder+k+".gz", 'ta') as output_file: #make this gzip.open for gzipped files
            for l in v:
                output_file.write(l+"\n")
            v.clear()
    print("...done writing")

total_lines=872493862 # figured this out using "zcat protein2ipr.dat.gz | wc -l"
with gzip.open(protein2ipr_filename, "tr") as input_file:
    domains_to_uniprot_list_dict=collections.defaultdict(list)
    counter=0
    for line in input_file:
        fields=line.split("\t")
        domain_accession=fields[3]
        if domain_accession[0:2]=="PF":
            domains_to_uniprot_list_dict[domain_accession].append(fields[0])
        if(counter%5000000==0 and counter>0):
            print("On line #{} ({:.2f}%) now".format(counter,100.0*counter/total_lines))
            if(counter%25000000==0):
                write_progress_and_clear_lists(domains_to_uniprot_list_dict)
        counter+=1
write_progress_and_clear_lists(domains_to_uniprot_list_dict)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1015485407980941, 'fraction_numerical': 0.007444907683144729, 'mean_word_length': 4.850609756097561, 'pattern_counts': {'":': 2, '<': 35, '<?xml version=': 0, '>': 44, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '594295', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 4757, 'n_tokens_neox': 4540, 'n_words': 908, 'path': 'co-occurrence_of_enzyme_domains_(preparing_data).ipynb', 'repo_name': 'tderond/CO-ED', 'size': '13439'}