Verify if all the components of a pathway are present in a genomes

Question

[dago]

I would like to verify if all the protein in specific pathways (e.g. secretion systems, two component regulatory system) are present in a genome.

Is there a way to do this automatically iterating the search for many genomes?

Answer 1

[umer.zeeshan.ijaz]

Cleaned up the post and moved to a Gist by Ram on 22-Apr-2022.

[dago]

Wow this looks like a really detailed work flow. I will give it a try but I am a bit dumb with the KEGG API. Any tutorial I could use to learn it?

[umer.zeeshan.ijaz]

Hi dago,

I figured it out by hit and trial method. Bioconductor's KEGGREST could be a good starting point to explore KEGG REST service further. Please also have a look at the HUManN workflow on how they incorporated MinPath. Here are a few more one-liners that you can use on a GENBANK file to understand the KEGG service better:

The API uses KEGG ENZYME database, which is an implementation of the Enzyme Nomenclature (EC Numbers) on the ExplorEnz database, and is maintained in the KEGG LIGAND relational database with additional annotation of reaction hierarchy, organism information, and sequence data links.

To use these one-liners on your GENBANK files, replace test.gbk with the name of the file you are using. First step is to extract a tab-separated list of only those contigs (once you have annotated them through PROKKA) which have enzymes in them. All other contigs are ignored

$ awk '/EC_number/{print gensub(" +/EC_number=\"(.*)\"","ec:\\1","g")}/^LOCUS/{print ";"$2}'  test.gbk | tr '\n' '\t' | tr ';' '\n' | awk 'NF>1{print $0}'
NODE_2925_length_923_cov_1.698808   ec:4.2.1.17 
NODE_4520_length_55408_cov_7.328635 ec:1.11.1.6 ec:1.2.1.22 ec:1.1.99.1 ec:2.5.1.18 
NODE_6043_length_823_cov_1.500607   ec:6.2.1.25 
NODE_918_length_270_cov_10.744445   ec:1.5.3.1  
NODE_312_length_2459_cov_13.763318  ec:4.3.1.17 
NODE_282_length_4935_cov_13.040932  ec:2.7.6.5  ec:3.1.7.2  ec:2.7.7.6  ec:2.7.4.8  
NODE_3153_length_49002_cov_8.110995 ec:1.4.3.19 ec:4.1.3.39 ec:5.1.1.8  ec:3.5.4.22 ec:1.2.1.26 ec:1.1.1.46 ec:1.1.1.1  ec:4.1.1.28 
NODE_241_length_4648_cov_13.109509  ec:4.2.1.20 
NODE_533_length_2973_cov_13.546249  ec:4.2.1.20
...

The following one-liner extracts the list of all enzymes found in the GENBANK file and uses rest-style KEGG API to generate names from EC numbers:

$ for I in $(awk '/EC_number/{print gensub(" +/EC_number=\"(.*)\"","\\1","g")}'  test.gbk | sort | uniq); do echo $(curl -s http://rest.kegg.jp/find/enzyme/$i) | awk '{$1=$1"\t"}1'; done
ec:1.1.1.108     carnitine 3-dehydrogenase
ec:1.11.1.15     peroxiredoxin; thioredoxin peroxidase; tryparedoxin peroxidase; alkyl hydroperoxide reductase C22; AhpC; TrxPx; TXNPx; Prx; PRDX
ec:1.1.1.158     Transferred to 1.3.1.98
ec:1.11.1.6  catalase; equilase; caperase; optidase; catalase-peroxidase; CAT
ec:1.1.1.169     2-dehydropantoate 2-reductase; 2-oxopantoate reductase; 2-ketopantoate reductase; 2-ketopantoic acid reductase; ketopantoate reductase; ketopantoic acid reductase
ec:1.11.1.7  peroxidase; lactoperoxidase; guaiacol peroxidase; plant peroxidase; Japanese radish peroxidase; horseradish peroxidase (HRP); soybean peroxidase (SBP); extensin peroxidase; heme peroxidase; oxyperoxidase; protoheme peroxidase; pyrocatechol peroxidase; scopoletin peroxidase; Coprinus cinereus peroxidase; Arthromyces ramosus peroxidase
...

For the extracted enzymes, we can list all the KEGG Ortholog (KO) groups each enzyme is part of, along with their detailed description. The KO system is the basis for representation for all proteins and functional RNAs that correspond to KEGG pathway nodes, BRITE hierarchy nodes, and KEGG module nodes:

$ for I in $(awk '/EC_number/{print gensub(" +/EC_number=\"(.*)\"","\\1","g")}'  test.gbk | sort | uniq); do echo $(curl -s http://rest.kegg.jp/link/ko/ec:$i | grep -Po '(?<=ko:).*' | sed -e "s/K/ko:K/g") | xargs -n 1 | xargs -I {} curl -s http://rest.kegg.jp/find/ko/{} | awk -v k=$i '{print "ec:"k"\t"$0}' ; done
ec:1.10.2.2 ko:K00411   UQCRFS1, RIP1, petA; ubiquinol-cytochrome c reductase iron-sulfur subunit [EC:1.10.2.2]
ec:1.1.1.1  ko:K00001   E1.1.1.1, adh; alcohol dehydrogenase [EC:1.1.1.1]
ec:1.1.1.1  ko:K00121   frmA, ADH5, adhC; S-(hydroxymethyl)glutathione dehydrogenase / alcohol dehydrogenase [EC:1.1.1.284 1.1.1.1]
ec:1.1.1.1  ko:K04072   adhE; acetaldehyde dehydrogenase / alcohol dehydrogenase [EC:1.2.1.10 1.1.1.1]
ec:1.1.1.1  ko:K11440   gbsB; choline dehydrogenase [EC:1.1.1.1]
ec:1.1.1.1  ko:K13951   ADH1_7; alcohol dehydrogenase 1/7 [EC:1.1.1.1]
ec:1.1.1.1  ko:K13952   ADH6; alcohol dehydrogenase 6 [EC:1.1.1.1]
ec:1.1.1.1  ko:K13953   adhP; alcohol dehydrogenase, propanol-preferring [EC:1.1.1.1]
ec:1.1.1.1  ko:K13954   yiaY; alcohol dehydrogenase [EC:1.1.1.1]
ec:1.1.1.1  ko:K13980   ADH4; alcohol dehydrogenase 4 [EC:1.1.1.1]
ec:1.1.1.100    ko:K00059   fabG; 3-oxoacyl-[acyl-carrier protein] reductase [EC:1.1.1.100]
ec:1.1.1.100    ko:K11610   mabA; beta-ketoacyl ACP reductase [EC:1.1.1.100]
ec:1.1.1.108    ko:K17735   lcdH, cdhA; carnitine 3-dehydrogenase [EC:1.1.1.108]
...

Similarly, for the extracted enzymes, we can also list all the known reactions these enzymes are a part of. Each reaction is identified by the R number and is linked to ortholog groups of enzymes enabling integrated analysis of genomic and chemical information.

$ for I in $(awk '/EC_number/{print gensub(" +/EC_number=\"(.*)\"","\\1","g")}'  test.gbk | sort | uniq); do echo $(curl -s http://rest.kegg.jp/link/reaction/ec:$i | grep -Po '(?<=rn:).*') | xargs -n 1 | xargs -I {} curl -s http://rest.kegg.jp/find/reaction/rn:{} | awk -v k=$i '{print "ec:"k"\t"$0}' ; done
ec:1.10.2.2 rn:R02161   Ubiquinol:ferricytochrome-c oxidoreductase; Ubiquinol + 2 Ferricytochrome c <=> Ubiquinone + 2 Ferrocytochrome c + 2 H+
ec:1.1.1.1  rn:R00228   acetaldehyde:NAD+ oxidoreductase (CoA-acetylating); Acetaldehyde + CoA + NAD+ <=> Acetyl-CoA + NADH + H+
ec:1.1.1.1  rn:R00623   primary_alcohol:NAD+ oxidoreductase; Primary alcohol + NAD+ <=> Aldehyde + NADH + H+
ec:1.1.1.1  rn:R00624   Secondary_alcohol:NAD+ oxidoreductase; Secondary alcohol + NAD+ <=> Ketone + NADH + H+
ec:1.1.1.1  rn:R00754   ethanol:NAD+ oxidoreductase; Ethanol + NAD+ <=> Acetaldehyde + NADH + H+
ec:1.1.1.1  rn:R01172   butanal:NAD+ oxidoreductase (CoA-acylating); Butanal + CoA + NAD+ <=> Butanoyl-CoA + NADH + H+
ec:1.1.1.1  rn:R02124   retinol:NAD+ oxidoreductase; Retinol + NAD+ <=> Retinal + NADH + H+
ec:1.1.1.1  rn:R02878   1-Octanol:NAD+ oxidoreductase; 1-Octanol + NAD+ <=> 1-Octanal + NADH + H+
...

We may also be interested in knowing which other organisms contain the same enzymes:

$ for I in $(awk '/EC_number/{print gensub(" +/EC_number=\"(.*)\"","\\1","g")}'  test.gbk | sort | uniq); do echo $(curl -s http://rest.kegg.jp/link/genome/ec:$i | grep -Po '(?<=genome:).*') | xargs -n 1 | xargs -I {} curl -s http://rest.kegg.jp/find/genome/genome:{} | awk -v k=$i '{print "ec:"k"\t"$0}' ; done
ec:1.10.2.2 genome:T01003   rno, RAT, 10116; Rattus norvegicus (Norway rat)
ec:1.10.2.2 genome:T01008   bta, BOVIN, 9913; Bos taurus (cow)
ec:1.1.1.1  genome:T00030   dme, DROME, 7227; Drosophila melanogaster (fruit fly)
ec:1.1.1.1  genome:T01001   hsa, HUMAN, 9606; Homo sapiens (human)
ec:1.1.1.1  genome:T01003   rno, RAT, 10116; Rattus norvegicus (Norway rat)
ec:1.1.1.1  genome:T01058   ecb, HORSE, 9796; Equus caballus (horse)
ec:1.1.1.100    genome:T00007   eco, ECOLI, 511145; Escherichia coli K-12 MG1655
ec:1.1.1.108    genome:T00035   pae, PSEAE, 208964; Pseudomonas aeruginosa PAO1
ec:1.11.1.15    genome:T01001   hsa, HUMAN, 9606; Homo sapiens (human)
ec:1.11.1.15    genome:T01003   rno, RAT, 10116; Rattus norvegicus (Norway rat)
ec:1.11.1.5 genome:T00005   sce, YEAST, 559292; Saccharomyces cerevisiae S288c
ec:1.1.1.157    genome:T00564   ckl, CLOK5, 431943; Clostridium kluyveri DSM 555
ec:1.11.1.6 genome:T00115   lpl, LACPL, 220668; Lactobacillus plantarum WCFS1

Best Wishes,
Umer

Answer 2

[5heikki]

A combination of KEGG's API, blastp and scripting should do it. If you're doing more small scale stuff, it might be feasible to just use http://www.kegg.jp/blastkoala/ and then the mapping tools provided at the site..

[dago]

Thanks for the answer. I tried that as well. But I was just wondering how sensitive this is. In blast KOALA the search is done using as reference the information for a Genus, and the information available for the genus I am working with are quite limited. Therefore, will I miss some protein known in other taxon to be involved in a specific process but not described in my genus?

[5heikki]

I don't think so:

The database files are generated from KEGG GENES as a collection of representative genomes by removing similar organisms at the species, genus or family level. When multiple members are present in each species/genus/family group, the first genome is taken as a representative genome. When the other members in the group contain different K numbers that are not present in the representative genome, those genes are added as if they are present in additional chromosomes or plasmids.

Stated right there in the blast koala site..

[dago]

Ok, got it, but for my Genus there are only two genomes available, will it change?

Answer 3

[Ram]

I'd take in a list of proteins as query and BLASTP against organism-filtered NR, verify identity and similarity in top 3 hits per query to ensure they fall under a threshold.

Answer 4

[Asaf]

You might want to check the methods of this paper: Long-term phenotypic evolution of bacteria : Nature : Nature Publishing Group

They used Flux Balance analysis to determine the substrates a bacterial strain can grow on.

[dago]

This is really cool, thanks for sharing. I am afraid this is to complex for my capacity..:P