The standard operating procedure of the DOE-JGI Microbial Genome Annotation Pipeline (MGAP v.4)
© Huntemann et al. 2016
Received: 18 April 2015
Accepted: 13 October 2015
Published: 26 October 2015
The Erratum to this article has been published in Standards in Genomic Sciences 2016 11:27
The DOE-JGI Microbial Genome Annotation Pipeline performs structural and functional annotation of microbial genomes that are further included into the Integrated Microbial Genome comparative analysis system. MGAP is applied to assembled nucleotide sequence datasets that are provided via the IMG submission site. Dataset submission for annotation first requires project and associated metadata description in GOLD. The MGAP sequence data processing consists of feature prediction including identification of protein-coding genes, non-coding RNAs and regulatory RNA features, as well as CRISPR elements. Structural annotation is followed by assignment of protein product names and functions.
Introduction and Requirements
The DOE-JGI Microbial Genome Annotation Pipeline performs structural and functional annotation of bacterial and archaeal genomes included into the Integrated Microbial Genome (IMG) system . Annotation consists of the identification of RNA and protein-coding genes and repeats, as well as the prediction of functions for each gene (product name assignment). The annotated microbial genomes datasets produced by the MGAP are integrated into IMG, where they can be analyzed or manually edited in the context of a comprehensive set of publicly available genomes.
Compared to the previous version, the pipeline contains several new steps (e.g. the addition of several modules that form a data preprocessing track), some tool replacements (e.g. Prodigal is now used for ab intio gene prediction instead of GeneMark and Metagene Annotator), as well as updated software and models (e.g. the rRNA Hidden Markov Models). So far this version of the pipeline has processed more than 800 microbial genomes with at least 10 new ones being processed on a daily basis. While no comprehensive reannotation of all genomes processed by the older version of the pipeline has been planned, point updates of the features affected by the pipeline changes are under way. Some of them, including CRISPR annotations and Pfam and TIGRfam assignments, have been finalized. The MGAP requires a multi-FASTA file of assembled nucleotide sequences as an input for gene calling. In addition, each sequence dataset submitted for annotation needs to be associated with an analysis project that has already been specified in the Genomes OnLine Database . Microbial genome annotation consists of three stages: sequence data preprocessing, feature prediction, and functional annotation. Feature prediction (which includes gene prediction and repeat identification) produces a GenBank file that does not have any functional information for the predicted genes. Subsequently, these genes are assigned functions and are integrated into IMG.
Procedure and Implementation
The MGAP stages and individual steps are further described below.
Sequence data preprocessing
As illustrated in Fig. 1(ii), feature prediction starts with the detection of CRISPR arrays, followed by non-coding RNA genes (tRNA, rRNA and other RNA genes), and finally the prediction of protein coding genes.
CRISPR elements are identified using the programs CRT  and PILER-CR v1.06 . For PILER-CR the maximum spacer length is set to 100 and the CRISPR element needs to have at least 5 repeats, which have at least 90 % identity to each other and at least 75 % identity to the consensus sequence of the repeat. For CRT the pipeline runs a modified version of the latest official CRT-CLI 1.2 version. Specifically, the modified CRT has the capability to read multi-FASTA files, detect truncated repeats at the ends of the contigs/scaffolds and deal with spacer artifacts and repeats that contain Ns. This version also executes checks for repeat and spacer length ratios, while the length and similarity checks are performed as part of “all vs. all” spacer and repeat comparisons. Additionally, the progression step of the sliding search window was reduced to 1, while threshold values and search ranges, which were strictly defined in the original software, can be changed from default values on the command line together with the new options and arguments.
In our modified CRT the default values for the minimum and maximum repeat lengths are set to 20 and 50 bp, respectively, while the minimum and maximum spacer lengths are set to 20 and 60 bp, respectively. The ratio of the spacer lengths to the repeat lengths have to be between 0.6 and 2.5. The search window is 7 bp long and an element needs to have at least 3 repeats that have a minimum of 70 % identity. The predictions from PILER-CR and CRT are concatenated and, in case of overlapping predictions, the CRT prediction is retained.
Protein-coding genes and non-coding RNA genes are identified using a combination of Hidden Markov Models and sequence similarity-based approaches. The first category of non-coding RNAs, tRNAs, are predicted using tRNAscan-SE 1.3.1 , and the best scoring predictions are selected. Ribosomal RNA genes (5S, 16S, 23S) are predicted using hmmsearch tool from the package HMMER 3.1b2  and a set of in-house curated HMMs derived from an alignment of full-length rRNA genes selected from IMG isolate genomes. Both tRNAscan-SE and hmmsearch use a domain-specific set of models for Bacteria and Archaea, which is selected based on the taxonomic information provided in the corresponding GOLD Analysis Project.
MGAP also predicts other non-coding RNAs and regulatory RNA features, such as riboswitches. With the exception of tRNAs and rRNAs, all models from Rfam 10.1  are used to search the genome sequences. For faster detection, sequences are first compared to a database containing all the ncRNA genes and other RNA features in the Rfam database using BLASTn, with a very loose e-value cutoff (1.0e −10). Subsequently, contigs/scaffolds with hits to this database are searched against Rfam covariance models using the program cmsearch from the INFERNAL 1.0.2 package .
The identification of protein-coding genes is performed using the Prodigal v2.6.2 ab initio gene prediction program . Overlaps between predicted features of different type (e. g. ncRNAs and protein-coding genes) get resolved based on an in-house curated set of rules. Every annotated gene is assigned a locus tag of the form PREFIX_SUFFIX where the prefix is the identifier of the GOLD Analysis Project associated with the genome dataset and the suffix is a number that identifies a certain gene on a particular sequence. This assignment scheme guarantees that each gene within a sequencing project gets a unique locus tag. The output of this stage is a GenBank format genome data file.
Every GenBank file by the pipeline or submitted by a user, must pass an additional validation step before it is forwarded to the next stage for functional annotation. The validation involves checking the file format, matching the gene coordinates to their translations and sequence lengths. The validation also removes phage PhiX sequences, which can be found as a common contaminant in isolate genomes sequenced with Illumina technology . These are identified by running BLASTn against the PhiX genome sequence with an e-value of 0.01 and 90 % identity. If a hit covers 80 % or more of a query contig/scaffold, the latter gets removed. As a final step the validation script also assesses the quality of a genome, which determines whether it will be included as a reference genome for taxonomic comparisons with other genomes and metagenomes. A genome could get marked as “low quality” and excluded from taxonomic reference database if (a) it lacks phylum-level taxonomic assignment or (b) its coding density (defined as total length of nucleotide sequence of predicted genes divided by the total length of nucleotide sequence) is less than 70 % or greater than 100 % or (c) there are more than 300 sequences per million base pair or (d) the number of genes per million base pair is less than 300 or greater than 1200. These values were set after manual analysis and benchmarking and are intended to prevent highly fragmented genomes and/or genomes with high rate of sequencing artifacts from being included in a taxonomic reference database.
After a genome dataset undergoes structural annotation, the resulting protein-coding genes are compared to protein families and the proteome of selected “core” genomes which are publicly available, and a protein product name is assigned to each gene as discussed below and illustrated in Fig. 1(iii).
COG & KOG assignment: protein sequences are compared to COG PSSMs obtained from the CDD database  using the program RPS-BLAST at an e-value cutoff of 1e–2, with the top hit retained. The alignment length needs to be at least 70 % of the consensus sequence length.
KEGG Orthology term assignment: Genes are associated with KO terms  as follows. First, the genes that can be unambiguously mapped to the entries in KEGG Genes database are assigned the KO terms associated with the corresponding KEGG gene. The gene to KEGG gene mapping is based on NCBI’s GI numbers and GeneIDs. For genes that are not mapped to KEGG genes, USEARCH is run against the database of KEGG genes by applying UBLAST . The results of this search are organized in a list of candidate KO assignments. KO terms are assigned to genes using a subset of this list, whereby the threshold is defined by an E-value cutoff of 1e–5, KO assignments are selected from the top 5 hits, with 30 % or better alignment sequence identity, and alignment percentage of at least 70 % over the length of the query gene and KEGG subject gene.
MetaCyc assignment: genes are associated with MetaCyc  reactions as follows. First, genes are mapped to KO terms as described above, whereby KO terms are associated with Enzyme Commission numbers using the KEGG KO term to Enzyme relationship provided by KEGG. Next, genes are associated with MetaCyc reactions via EC numbers.
Pfam & TIGRfam assignments: protein sequences are searched against Pfam  and TIGRfam  databases using HMMER 3.0. For TIGRfam, the noise cutoff (−−cug_nc) is used, with hits below the trusted cutoff and at/above the noise cutoff flagged as “marginal”. For Pfam, the gathering threshold (−−cut_ga) is used inside the pfam_scan.pl script. The script also helps resolving overlaps between hits to Pfam models from the same clan in order to generate final Pfam assignments.
InterPro Scan: Additional protein family annotations for SMART, PrositeProfiles, PrositePatterns, and SuperFamily are provided by InterPro Scan (run with default parameters) .
Method 1 is based on the BLAST-able database of gene sequences that have been associated manually with IMG terms. Gene sequences of new genomes are compared against this database: a gene g with hit to a sequence associated with IMG term t satisfying the following criteria will be assigned this IMG term: top hit, e-value < = 1e-5, > = 90 % identity, alignment > = 80 % on both query and subject sequences, smallest to largest sequence length ratio of query and subject sequence > = 70 %.
Method 2 is applied to genes that are not assigned an IMG terms using Method 1 and is based on a set of rules devised by domain experts at JGI for mapping functional annotations (COG, Pfam, TIGRfam) to IMG terms. An example of such a rule is: “assign IMG Term 6 (replicative DNA helicase loader DnaB) to a gene if the gene is annotated with COG3611 (replication initiation/membrane attachment protein)”.
- 3.Method 3 is applied to genes that are not assigned an IMG term using Methods 1 and 2, and is based on gene bi-directional best hits (BBHs) across all genomes available in IMG. For a gene g:
Get g’s top 5 BBH genes satisfying the following conditions: sequence alignment length > = 70 %, percent identity > = 25 %. No IMG terms can be assigned to gene g unless there are at least 5 BBH genes satisfying these conditions.
Let Set T be the set of all manually assigned IMG terms (i.e., not automatically populated terms) of any of the 5 BBH genes above. Check each term T1 in Set T: (i) if the 5 BBH genes have conflicting term assignments (e.g., some were assigned term T1, while others were assigned term T2), then no terms in Set T can be assigned to gene g; (ii) if there are no conflicting IMG term assignments and at least 2 of the 5 BBH genes are associated with term T1, then assign T1 to gene g; (iii) if there are no conflicting IMG term assignments but fewer than 2 of the 5 BBH genes are associated with the same IMG term, then repeat this step with top 10 BBH genes.
Bi-directional best hits are computed between protein sequences of pair of genomes using USEARCH by applying UBLAST with a nominal e-value cutoff of 1e–2. An effective database size (−ka_dbsize 700,000,000) is used in order to make the e-values comparable across pairs of genome computations.
Protein product names
ITERM:xxxxx: The source of the gene product name is its associated IMG term xxxxx;
TIGRxxxxx: The source of the gene product name is its associated TIGRfam xxxxx;
COGxxxx: The source of the gene product name is its associated COGxxxx;
PFAMxxxx: The source of the gene product name is the associated PFAMxxxx.
IMG Protein product names are assigned to genes in two stages: (i) product names are assigned based on IMG terms whenever they are available, (ii) if IMG terms are not available then protein family associations for genes in individual genome datasets are employed for assigning product names. Protein product names assignments based on IMG terms rely on protein sequence similarities between the genes of the new dataset and genes of all other genomes in the IMG data warehouse.
If gene g is associated with one or more IMG terms, then the IMG term(s) becomes the new IMG product name(s) of g.
For genes that are not associated with IMG terms based on Methods 1–3 described above, and therefore do not have a product name based on IMG terms, assignment of TIGRfam names as product names is attempted: the gene without a product name is assigned a name of a TIGRfam if it has a TIGRfam hit. If a gene has a hit to only one TIGRfam, the name of this TIGRfam is assigned; if more than 1 TIGRfam is assigned, the name of a TIGRfam of the type “equivalog” is assigned.
For genes that were not associated with a product name using IMG terms or TIGRfam names, product names are assigned based on the name of their COG hit. If the COG name is “uncharacterized conserved protein” or contains “predicted”, the name has the format “COG.cog_name, COG.cog_id”.
For genes that were not associated with a product name using IMG terms, TIGRfam or COG names, product names are assigned based on the name of their Pfam hit, where the product name is a concatenation of Pfam family description (attribute “description” in pfam_family) with “protein”. If a protein has hits to multiple Pfams, their descriptions are concatenated with “/” as a separator and a word “protein” added in the end.
A translation table for protein product names based on IMG terms, TIGRfam, COG and Pfam descriptions is employed in order to ensure that the product names are compatible with GenBank requirements when the datasets are submitted to GenBank.
Functional annotation sources
COG 2014 Version (November 2014)
KEGG Release 71.0, July 2014
MetaCyc Release 18.1, June 2014
PFAM 28.0, May 2015
TIGRfam Release 14.0, January, 2014
InterPro Data Release 48
Prediction of protein sequence features
Signal peptide feature prediction employs SignalP 4.1 . The model used is determined by the gram stain annotation field for the genome (gram+, gram-, Euk). If the gram stain field is not specified in the IMG data, assume it’s gram- (gram negative).
Transmembrane helices are predicted using TMHMM2.0c .
Both tools get executed with the default parameters.
The DOE-JGI MGAP performs annotation for bacterial and archaeal genomes. The pipeline consists of custom scripts and publicly available tools. Consistency and reproducibility of the results produced by the MGAP depend on the tools and annotation resources used in the pipeline. Thus, updated versions of resources such as Rfam, Pfam, and KEGG may improve the breadth and depth of functional annotations.
In order to apply the DOE-JGI MGAP on their datasets, users need to first specify their analysis projects in GOLD and then provide their genome datasets via IMG’s data submission site.
We will continue to extend the MGAP with the goal of improving the identification and characterization of genes in the genome datasets it processes.
This work is funded by Director, Office of Science, Office of Biological and Environmental Research, Life Sciences Division, U.S. Department of Energy (Contract No. DE-AC02-05CH11231).
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Markowitz VM, Chen IM, Palaniappan K, Chu K, Szeto E, Pillay M, et al. IMG 4 version of the integrated microbial genomes comparative analysis system. Nucleic Acids Res. 2014;42:D560–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Reddy TB, Thomas AD, Stamatis D, Bertsch J, Isbandi M, Jansson J, et al. The Genomes OnLine Database (GOLD) v. 5: a metadata management system based on a four level (meta)genome project classification. Nucleic Acids Res. 2015;43:D1099–106.View ArticlePubMedPubMed CentralGoogle Scholar
- Morgulis A, Gertz EM, Schäffer AA, Agarwala R. A fast and symmetric DUST implementation to mask low-complexity DNA sequences. J Comput Biol. 2006;5:1028–40.View ArticleGoogle Scholar
- Bland C, Ramsey TL, Sabree F, Lowe M, Brown K, Kyrpides NC, et al. CRISPR recognition tool (CRT): a tool for automatic detection of clustered regularly interspaced palindromic repeats. BMC Bioinformatics. 2007;8:209.View ArticlePubMedPubMed CentralGoogle Scholar
- Edgar RC. PILER-CR: fast and accurate identification of CRISPR repeats. BMC Bioinformatics. 2007;8:18.View ArticlePubMedPubMed CentralGoogle Scholar
- Lowe TM, Eddy SR. tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. 1997;25:955–64.View ArticlePubMedPubMed CentralGoogle Scholar
- Eddy SR. Accelerated profile HMM searches. PLoS Comput Biol. 2011;7, e1002195.View ArticlePubMedPubMed CentralGoogle Scholar
- Griffiths-Jones S, Moxon S, Marshall M, Khanna A, Eddy SR, Bateman A. Rfam: annotating non-coding RNAs in complete genomes. Nucleic Acids Res. 2005;33:D121–4.View ArticlePubMedPubMed CentralGoogle Scholar
- Nawrocki EP, Kolbe DL, Eddy SR. Infernal 1.0: inference of RNA alignments. Bioinformatics. 2009;25:1335–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Hyatt D, Chen GL, Locascio PF, Land ML, Larimer FW, Hauser LJ. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics. 2010;11:119.View ArticlePubMedPubMed CentralGoogle Scholar
- Mukherjee S, Huntemann M, Ivanova N, Kyrpides NC, Pati A. Large-scale contamination of microbial isolate genomes by Illumina PhiX control. Stand Genomic Sci. 2015;10:18.View ArticlePubMedPubMed CentralGoogle Scholar
- Marchler-Bauer A, Anderson JB, Derbyshire MK, DeWeese-Scott C, Gonzales NR, Gwadz M, et al. CDD: a conserved domain database for inter-active domain family analysis. Nucleic Acids Res. 2007;35:D237–40.View ArticlePubMedPubMed CentralGoogle Scholar
- Kanehisa M, Goto S, Sato Y, Kawashima M, Furumichi M, Tanabe M. Data, information, knowledge and principle: back to metabolism in KEGG. Nucleic Acids Res. 2014;42:D199–205.View ArticlePubMedPubMed CentralGoogle Scholar
- Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010;26:2460–1.View ArticlePubMedGoogle Scholar
- Caspi R, Altman T, Billington R, Dreher K, Foerster H, Fulcher CA, et al. The MetaCyc database of metabolic pathways and enzymes and the BioCyc collection of pathway/genome databases. Nucleic Acids Res. 2014;42:D459–71.View ArticlePubMedPubMed CentralGoogle Scholar
- Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate J, Boursnell C, et al. The Pfam protein families database. Nucleic Acids Res. 2012;40:D290–301.View ArticlePubMedPubMed CentralGoogle Scholar
- Selengut JD, Haft DH, Davidsen T, Ganapathy A, Gwinn-Giglio M, Nelson WC, et al. TIGRFAMs and Genome Properties: tools for the assignment of molecular function and biological process in prokaryotic genomes. Nucleic Acids Res. 2007;35:D260–4.View ArticlePubMedPubMed CentralGoogle Scholar
- Jones P, Binns D, Chang HY, Fraser M, Li W, McAnulla C, et al. InterProScan 5: genome-scale protein function classification. Bioinformatics. 2014;30:1236–40.View ArticlePubMedPubMed CentralGoogle Scholar
- Chen IM, Markowitz VM, Chu K, Anderson I, Mavromatis K, Kyrpides NC, et al. Improving microbial genome annotations in an integrated database context. PLoS One. 2013;8, e54859.View ArticlePubMedPubMed CentralGoogle Scholar
- Petersen TN, Brunak S, von Heijne G, Nielsen H. SignalP 4.0: discriminating signal peptides from transmembrane regions. Nat Methods. 2010;10:785–6.Google Scholar
- Krogh A, Larsson B, von Heijne G, Sonnhammer EL. Predicting transmembrane protein topology with a hidden Markov model: application to complete genomes. J Mol Biol. 2001;3:567–80.View ArticleGoogle Scholar