Standardized genome-wide function prediction enables comparative functional genomics: a new application area for Gene Ontologies in plants

Abstract

Background

Genome-wide gene function annotations are useful for hypothesis generation and for prioritizing candidate genes potentially responsible for phenotypes of interest. We functionally annotated the genes of 18 crop plant genomes across 14 species using the GOMAP pipeline.

Results

By comparison to existing GO annotation datasets, GOMAP-generated datasets cover more genes, contain more GO terms, and are similar in quality (based on precision and recall metrics using existing gold standards as the basis for comparison). From there, we sought to determine whether the datasets across multiple species could be used together to carry out comparative functional genomics analyses in plants. To test the idea and as a proof of concept, we created dendrograms of functional relatedness based on terms assigned for all 18 genomes. These dendrograms were compared to well-established species-level evolutionary phylogenies to determine whether trees derived were in agreement with known evolutionary relationships, which they largely are. Where discrepancies were observed, we determined branch support based on jackknifing then removed individual annotation sets by genome to identify the annotation sets causing unexpected relationships.

Conclusions

GOMAP-derived functional annotations used together across multiple species generally retain sufficient biological signal to recover known phylogenetic relationships based on genome-wide functional similarities, indicating that comparative functional genomics across species based on GO data holds promise for generating novel hypotheses about comparative gene function and traits.

Background

Phenotypes and traits have long been the primary inspiration for biological investigation. Phenotypes are the result of a complex interplay between functions of genes and environmental cues. In an effort to organize and model gene functions, various systems of classification have been developed including systems like KEGG, which is focused on protein function including gene activities superimposed on metabolic pathways [1]. Other such systems include the various Cyc databases, MapMan, and the Gene Ontologies (GO), a vocabulary of gene functions organized as a directed acyclic graph, which makes it innately tractable for computational analysis [2–4].

GO-based gene function annotation involves the association of GO terms to individual genes. Functions may be assigned to genes on the basis of different types of evidence for the association. For example, functional predictions can be inferred from experiments (EXP), expression patterns (IEP), and more [5]. Computational pipelines are often used to generate functional predictions for newly sequenced genomes, where the genome is first sequenced and assembled, then gene structures (gene models) are predicted, then functions are associated with those gene predictions. Genome-wide gene function prediction datasets are frequently used to analyze gene expression studies, to prioritize candidate genes linked to a phenotype of interest, to design experiments aimed at characterizing functions of genes, and more [6–8]. How well a gene function prediction set models reality is influenced by how complete and correct the underlying genome assembly and gene structure annotations are, as well as by how well the software used to predict functions performs.

GOMAP (the Gene Ontology Meta Annotator for Plants) is a gene function prediction pipeline for plants that generates high-coverage and reproducible functional annotations [9]. The system uses multiple functional prediction approaches, including sequence similarity, protein domain presence, and mixed-method pipelines developed to compete in the Critical Assessment of Function Annotation (CAFA) Challenge [10], a community challenge that has advanced the performance of gene function prediction pipelines over the course of 5 organized competitions [11].

We previously annotated gene functions for the maize B73 genome and demonstrated that GOMAP’s predicted functions were closer to curated gene-term associations from the literature than those of other community functional annotation datasets, including those produced by Gramene (Ensembl pipeline) and Phytozome (Interpro2GO pipeline) [12]. Using the newly containerized GOMAP system [9], we report here the functional annotation of 18 plant genomes across the 14 crop plant species listed in Table 1 and report comparisons of performance based on comparison to gold standard gene function datasets, where possible.

Given these multiple annotations across various plant species, we next considered whether these datasets could be used together for comparative functional genomics in plants. We describe here a simple and crude method by which we used gene function annotations to generate dendrograms of genome-level similarity in function. This idea is similar to that of Zhu et al., who determined the evolutionary relationships among microorganisms based on whole-genome functional similarity [50]. Here we expand on that approach, analyzing genome-wide GO assignments to generate parsimony and distance-based dendrograms (see Fig. 1 for process overview). We compared these with well-established species phylogenies (Fig. 2) to determine whether trees derived from gene function show any agreement with evolutionary histories, taking agreement between generated dendrograms and known evolutionary histories to be evidence that sufficient comparative biological signal exists to begin to use GO functional annotations across multiple plant genomes for comparative functional genomics investigations.

Results of Analyses

Overview

As shown in Fig. 1, gene function annotation sets were created and compiled for each genome. For those with existing annotation sets available on Gramene or Phytozome [51,52], the datasets were compared. From there, matrices that included genomes as rows and terms as columns were generated. These were used directly to build parsimony trees or to create distance matrices for NJ tree construction [53–55]. In subsequent analyses, jackknifing was used to remove terms (columns) or to remove genomes (rows) to map the source of signal for tree-building results [56].

Functional annotation sets produced

Table 2 reports quantitative attributes of each of the annotation sets. In summary, GOMAP covers all annotated genomes with ≥1 annotation per gene and provides between 3.8 and 12.1 times as many annotations as Gramene or Phytozome.

Quality evaluation of gene function predictions is not trivial and is approached by different research groups in different ways. Most often datasets are assessed by comparing the set of predicted functions for a given gene to a gold standard consisting of annotations that are assumed to be correct. This assumption of correctness can be based on any number of criteria. Here we used as our gold standard dataset all annotations present in Gramene63 that had a non-IEA (non-inferred by electronic annotation) evidence code; i.e., we used only annotations that had some manual curation. This enabled us to assess 10 of the genomes described in Table 2. It is perhaps noteworthy that the IEA and non-IEA annotation sets from Gramene63 frequently contain overlaps, indicating that some of the predicted annotations were manually confirmed afterwards by a curator and that in such cases, a new annotation was asserted with the new evidence code rather than simply upgrading the evidence code from IEA to some other code, thus preserving the IEA annotations in Gramene63 that are produced by the Ensembl analysis pipeline [57], a requirement for comparing GOMAP-produced IEA datasets to the IEA datasets produced by the Ensembl pipeline.

A general limitation of using gold standards for quality evaluation is that they can never be assumed to be complete, and therefore false-positive results in the prediction cannot be distinguished from false-negative results in the gold standard. In other words, is gene X, function Y truly a wrong prediction or has it simply not yet been discovered experimentally? This problem is laid out in more detail in [58]. As a consequence, the quality of larger prediction sets will be systematically underestimated compared to smaller ones, and this effect is exacerbated the more incomplete the gold standard is.

There are many different metrics that have been used to evaluate the quality of predicted functional annotations. For the maize B73 GOMAP annotation assessment in [12], we had used a modified version of the hierarchical evaluation metrics originally introduced in [59] because they were simple, clear, and part of an earlier attempt at unifying and standardizing GO annotation comparisons [60]. In the meantime, Plyusnin et al. published an approach for evaluating different metrics showing variation among the robustness of different approaches to quality assessment [61]. On the basis of their recommendations, we use here the SimGIC2 and term-centric area under precision-recall curve (TC-AUCPCR) metrics. We also evaluated with the F_max metric, simply because it is widely used (e.g., by [10]), even though according to Plyusnin et al., it is actually a flawed metric [61]. Results of the quality assessments for the 10 genomes where a gold standard was available are reported in Table 3 and Supplementary Fig. S2. While evaluation values differ between metrics and the scores are not directly comparable, a few consistent patterns emerge: GOMAP annotations are almost always better than Gramene and Phytozome annotations in the cellular component and molecular function aspect, with the only 3 exceptions being the molecular function aspect for T. aestivum using the TC-AUCPCR and the F_max metric and the cellular component aspect for M. truncatula A17 using the F_max metric. Conversely, GOMAP predictions achieve consistently lower quality scores in the biological process aspect with the exception of B. dystachion, O. sativa, and S. bicolor with the TC-AUCPR metric. Generally, annotations that are better in 1 aspect are also better in the other 2 aspects, but the ranking of annotations does not necessarily hold across metrics. The Phytozome annotation for O. sativa is an outlier in terms of its comparative quality, potentially because it is based on a modified structural annotation that differs substantially from the gold standard and the other annotations under comparison.

Phylogenetic tree analyses

With the comparative quality of gene function predictions in hand, we approached the question of whether the datasets could be used together for comparative functional analysis across all genomes. As a simple first step, we began to work toward understanding the degree to which trees built on the basis of gene functions agree with known, well-documented evolutionary relatedness. We constructed neighbor-joining (NJ) and parsimony trees of the 18 plant genomes and visulized them using iTOL [62]. The 2 tree topologies, rooted at P. lambertiana, were compared to one another and to the topology of the expected tree (Fig. 2). For both the NJ (Fig. 3A) and parsimony trees (Fig. 3B), 1 common difference is noted: S. bicolor is not at the base of the Z. mays clade as expected and is clustered with B. distachyon instead. Notable differences between the NJ and parsimony tree are the following: C. sativa appears at the base of the eudicots instead of G. raimondii in the NJ tree, while G. raimondii is grouped with C. sativa and A. hypogaea is grouped with G. max in the parsimony tree. Second, O. sativa was expected to be at the base of the Bambusoideae, Oryzoideae, and Pooideae (BOP) clade but appears at the base of Z. mays in the NJ tree, and at the base of all angiosperms in the parsimony tree. Differences among relationships within the Z. mays clade constaining B73, PH207, W22, and Mo17 were disregarded given the high degree of similarity across annotation sets and the fact that these relationships are not clear given the complex nature of within-species relationships.

Owing to differences between the function-based dendrograms and the expected tree, jackknifing analysis was carried out by removing terms (columns in underlying datasets) to determine the degree to which the underlying datasets support specific groupings based on functional term assignments. This analysis was carried out for both NJ and parsimony trees. First, trees were generated by omitting 5%–95% of the dataset in increments of 5 to determine the threshold at which the tree topologies deviated from those generated using the full dataset. That threshold was reached at 45% for both NJ and parsimony; therefore, we used trees generated with 40% of the data removed for reporting branch support values for the topology (Fig. 3). Comparing the 2 trees, the parsimony topology was not as solid as that of the NJ at jackknife values ≤40%. On the basis of this robustness for NJ tree-building in general, we carried out all subsequent analyses using NJ tree-building methods.

We considered investigating the effect of using 1 GO aspect to generate our NJ tree. In other words, we generated the NJ trees using cellular component GO terms, molecular function GO terms, and biological process GO terms separately (Supplementary Fig. S3). Of the 14,303 total GO terms, 1,524 are cellular component terms, 3,926 are molecular function terms, and 8,853 are biological process terms. Of the 3 single-aspect phylogenetic trees, the one built using molecular function terms is the closest to our NJ tree obtained using all GO terms in our datasets (Fig. 3A). The only difference is that A. hypogaea and G. max are clustered in the molecular function tree, while they are not in our NJ tree Fig. 3A. In the cellular component tree, G. raimondii and C. sativa are clustered together when they are not in the NJ tree with all GO aspects (Fig. 3A). Also, O. sativa is at the base of the monocots just like in the expected tree, but not in the NJ tree Fig. 3A. In the biological process tree, O. sativa is at the base of the angiosperms and there is no clear separation between monocots and dicots. In all 3 single-aspect phylogenetic trees and our all-aspect NJ tree, A. hypogaea is never placed at the base of the NPAAA clade. Also, B. distachyon and S. bicolor are always clustered together. Overall, the topologies constructed using 1 GO term aspect at a time are close to that of our NJ tree, such that not one GO term aspect alone restored the topology of the expected tree.

To map the source of discrepancies to specific gene annotation sets, we generated various NJ trees excluding 1 genome each time, an additional tree with both Medicago genomes excluded simultaneously, and another with all Z. mays genomes excluded simultaneously. To exemplify this, see the monocot clade in Fig. 2 and the lower (monocot) clade in Fig. 3A. When the NJ tree was generated, 2 species are misplaced: S. bicolor and O. sativa. As shown in Fig. 4a, removal of O. sativa corrects 1 error (itself) but does not correct the errant grouping of S. bicolor with B. distachyon. In Fig. 4b, it is shown that the removal of S. bicolor corrects the errant grouping of itself and B. distachyon, but O. sativa placement remains incorrect. However, as shown in Fig. 4c, the removal of B. distachyon generates a tree where all relationships are in agreement with known species-level relationships. (Note well: all individual annotation sets were progressively removed, not just the 3 shown in the example.)

With this observation in hand, we sought to determine the minimum number of genomes that could be removed to create a tree that matched the expected tree topology. All possible combinations of removing 0–5 genomes to restore the topology were tested, and 10 combinations of minimum amount of genomes to be removed were obtained. The removal of 4 genomes was required to generate function-based trees consistent with known phylogenetic relationships. Of the 10, we selected the 1 that had the genomes that were most frequently part of a solution (O. sativa, 8; B. distachyon, 7; C. sativa, 6; A. hypogaea, 5; S. bicolor, 4; G. raimondii, 4; G. max, 4; T. aestivum) to show in this article (the other combinations can be found in our publicly available dataset). To elaborate, the genomes removed here are O. sativa, B. distachyon, C. sativa, and A. hypogaea (Fig. 5). Jackknifing analysis was also carried out for this dataset with support shown. Branch support is generally higher than that for the full dataset (i.e., branch support is higher in Fig. 5 than in Fig. 3A), and removing genomes that are causing variations seems to stabilize the tree.

Potential causes of unexpected groupings

As a first step toward explaining discrepancies between known evolutionary relationships and those resulting from comparative analysis of genome-wide gene function predictions, we assessed the quality of each genome assembly and structural annotation set using GenomeQC [63]. Tables 4 and 5 and Figs. 6 and 7 represent the resulting assembly quality, structural annotation measures of quality, and proportion of single-copy BUSCOs [64] that were generated. Although these analyses make evident that the species annotated are comparatively different in both natural genome characteristics and in assembly and annotation quality aspects, it is not the case that the 4 species responsible for deviations between the functional annotation dendrograms and known phylogenetic relationships (i.e., C. sativa, A. hypogaea, O. sativa, and B. distachyon) create these discrepancies owing to issues of genome assembly and/or annotation quality. One potential for some explanation is in relation to C. sativa, which is the only genome that has an assembly length larger than the expected (see Table 4), and a comparatively large proportion of missing BUSCOs in the assembly (see Fig. 6). Similarly, for A. hypogaea and O. sativa, there is a large proportion of missing BUSCOs in the annotations (see Fig. 7).

Discussion

In this study, we used the GOMAP pipeline to produce whole-genome GO annotations for 18 genome assembly and annotation sets from 14 plant species [9]. Assessments of the number of terms predicted, as well as the quality of predictions, indicate that GOMAP functional prediction datasets cover more genes, contain more predictions per gene, and are of similar quality to prediction datasets produced by other systems, thus supporting the notion that these high-coverage datasets are a useful addition for researchers who are interested in genome-level analyses, including efforts aimed at prioritizing candidate genes for downstream analyses. Given that we can now produce high-quality, whole-genome functional annotations for plants in a straightforward way, we intend to produce more of these over time (indeed we recently annotated Vitis vinifera [76], Brassica rapa [77], Musa acuminata[78], Theobroma cacao [79], Coffea canephora [80], Vaccinium corymbosum [81], Solanum lycopersicum [82], and Solanum pennellii [83]).

With 18 genome functional annotations in hand, we sought to determine whether and how researchers could use multispecies GO annotation datasets to perform comparative functional genomics analyses. As a proof of concept, we adapted phylogenetic tree-building methods to use the gene function terms assigned to genes represented by the genomes to build dendrograms of functional relatedness and hypothesized that if the functions were comparable across species, the resulting trees would closely match evolutionary relationships. To our delight and surprise, the NJ and parsimony trees (Fig. 3) did resemble known phylogenies, but were not exact matches to broadly accepted phylogenetic relationships.

After removing the minimum number of genomes that resulted in restoration of the expected evolutionary relationships, we found that the individual species that may be responsible for the discrepancies observed in Fig. 3 were C. sativa, A. hypogaea, O. sativa, and B. distachyon. We hypothesize that the following could account for such errant relationships:

1. Quality of sequencing and coverage assembly: genomes of similarly high sequence coverage that have excellent gene calling would be anticipated to create the best source for functional annotation. Genomes of comparatively lower, or different, character would be anticipated to mislead tree-building and other comparative genomics approaches.

2. Shared selected or natural traits: species that have been selected for, e.g., oilseeds may share genes involved in synthesis of various oils. Other shared traits would be anticipated to cause similarities for species with those shared traits.

3. Lack of good representation of diverse plant biology aspects in the GO graph: most plant-specific GO terms were derived from functional analysis of 1 model species, Arabidopsis thaliana. This single source for presence of plant-specific functions limits the graph from containing unique functional aspects of plant biology represented in other species’ genomes. In addition, this limitation could lead to the assignment of unknown or errant functions based on a lack of closely related and/or plant-specific terms.

4. Use of a simple method of tree-building based on the presence or absence of gene function terms: the method that we devised and describe here is not sophisticated enough to make full use of information in the GO graphs such that we recover the full detail of the species’ evolutionary histories from the simple method.

To consider the first of these, we looked at genome assembly and annotation quality metrics (see Tables 4 and 5 and Figs. 6 and 7). For B. distachyon we could find no compelling evidence that assembly structural or functional annotation quality differed significantly from all others, except in the case of C. sativa, where we noted that the assembly length exceeded the predicted genome size based on C-values for genome sizes reported previously [67]. In this case, the fact that the C. sativa line sequenced is not inbred [84] may be responsible for the inflated assembly size relative to what is expected. This means that in the assembly, there are likely regions where alleles between chromosomes do not align, which would inflate the overall length of the assembly. In addition, the assembly misses a large proportion of BUSCO genes compared to most other genomes included in this analysis. Indeed, the comparatively low-quality assembly for Cannabis genome has been noted by others [85], and our preliminary investigations indicate that the assembly length is in fact longer than expected.

In an attempt to better understand conflicting phylogenetic signals that could be caused by the second potential cause, i.e., shared selected or natural traits, we mapped all GO terms that exist in our binary matrix and traced character history (presence/absence) on the nodes and leaves of our expected evolutionary tree using the software Mesquite version 3.61 [86]. These data summarize the gain or loss of each GO term across the species described in this article and can be found in our GitHub repository. We carried out a number of simple experiments to reveal which terms could be causal for errant relationships (e.g., dropping all unique terms from the B. distachyon dataset, reconstructing the term states at nodes that should be where B. distachyon should occur) and could not identify any biologically compelling patterns. (Because these analyses were not fruitful, they were not specifically included in our Methods section, although we do include the input datasets here in section “Availability of Source Code and Supporting Data,” for others to consider and peruse independently.)

A likely limitation in this analysis is the effect of the third potential cause of discrepancies between our generated phylogenetic trees and the expected topology: a deficiency of terms that describe diverse gene functions across the diversity of plant species. Because most GO terms specific to plant biology are likely derived from Arabidopsis, a model dicot species, gene functions unique to other species are expected to be missing from the GO graphs [87]. This source of error will only be corrected over time as gene functions unique to diverse plant species are populated into the GO graph.

We consider the most likely explanation for observed discrepancies between the known evolutionary phylogenies and dendrograms created on the basis of GO terms describing gene function to be a result of the fourth explanation: the simplicity of the tree-building models and methods we used for these analyses. Because the tree-building and analytics described in this article were based on the presence/absence of GO terms, novel terms are highly influential on the outcomes of the analysis and the number of times a term is used does not influence the outcome at all. In contrast, plant genomes are notable for having many duplicated genes as a result of whole-genome and segmental duplications over evolutionary history, so these duplications are in fact a feature of and marker for what happened to that genome over time. Therefore, using presence/absence of GO terms where the counts of term occurrences are not weighted may be too simple to get at the genuine biological complexities represented in any given plant genome. Our simplistic demonstration of the utility of GO datasets for comparative functional genomics shows that more sophisticated methods are promising for comparative functional genomics analyses.

It should be noted that comparative analyses using gene functions are not completely absent from the literature—although they are absent for large genome comparisons in plants. An example of such existing comparative use of GO is one where a given tree topology was used to look for gains and losses of functions mapped to independently derived trees, which was reported by Schwacke et al. [88]. They report, as an example of their method, an analysis of gene loss in Cuscuta, a parasitic plant, based on analysis using the Mapman ontology. In their work, they showed considerable loss of genes, which is a hallmark of the parasitic lifestyle. Our efforts differ in that we used the functions directly to infer tree structures as a demonstration that sufficient biological signal is present in GO-based datasets of genome-wide function prediction to reproduce known biological relationships. The method that we used was quick and dirty, and we anticipate that refinements in approach that consider multiple copies of genes, as well as using different types of graph and network representations beyond tree structures, are logical next steps for refining the use of GO terms for comparative functional genomics analyses in plants. With that in mind, we look forward not only to developing systems to support GO-based comparative functional genomics tools but also to seeing the tools that other research groups will develop to approach the use of these datasets to formulate novel comparative functional genomics hypotheses.

Methods

Acquiring input datasets

For each of the 18 genomes listed in Table 1, information on how to access input annotation products is listed by DOI. For each, 1 representative translated peptide sequence per protein-coding gene was selected and used as the input for GOMAP, a gene function prediction tool for plants that is actively maintained, updated, and versioned. Details of how GOMAP annotations are derived including the specificity of component datasets and which terms are retained are described elsewhere [9, 12]. In brief, GOMAP annotations are a combination of the annotations from multiple sources. GOMAP combines the annotations from all the sources and removes the less specific annotations that could be inferred from the more specific annotations, keeping only the most specific terms for each gene that cannot be inferred from other terms (i.e., only leaf terms). Unless the authors of the genome provided a set of representative sequences designated as canonical, we chose the longest translated peptide sequence as the representative for each gene model. In general, non-IUPAC characters and trailing asterisks were removed from the sequences, and headers were simplified to contain only non-special characters. The corresponding script for each dataset can be found at the respective DOI. On the basis of this input, GOMAP yielded a functional annotation set spanning all protein-coding genes in the genome. Using the Gene Ontology version releases/2020-10-09, this functional annotation set was cleaned up by removing duplicates, annotations with qualifiers (NOT, contributes_to, colocalizes_with; column 4 in the GAF 2.1 format), and obsolete GO terms. Any terms containing alternative identifiers were merged to their respective main identifier, uncovering a few additional duplicates, which were also removed. Supplementary Table S2 shows the number of annotations removed from each dataset produced.

To compare the quality of GOMAP predictions to currently available functional predictions from Gramene and Phytozome, we downloaded IEA annotations from Gramene (version 63) [51, 89] and Phytozome (version 12) [52, 90] for each species with functional annotations of the same genome version. These datasets were cleaned as above. Similarly cleaned non-IEA annotations from Gramene63 served as the gold standard wherever they were available. More detailed information on how these datasets were accessed can be found at [91].

Quantitative and qualitative evaluation

The number of annotations in each clean dataset was determined and related to the number of protein-coding genes (based on transcripts in the input FASTA file). This was done separately for each GO aspect as well as in total.

The ADS software version published in [61] is available from [92]. We used version b6309cb (also included in our code as a submodule) to calculate SimGIC2, TC-AUCPCR, and F_max quality scores. To provide the information content required for the SimGIC2 metric, the Arabidopsis GOA from [93] was used in version 2021-02-16.

Cladogram construction

An NJ tree was constructed on the basis of the generated pairwise distance matrix using PHYLIP (PHYLIP, RRID:SCR_006244) [53]. Additionally, term sets S of all genomes were combined into a binary matrix (with rows corresponding to genomes and columns corresponding to GO terms, values of 0 or 1 indicating whether a term is present or absent in the given set). PHYLIP pars was used to construct a parsimony tree from this binary matrix.

P. lambertiana, a gymnosperm, was included in the dataset as an outgroup to the angiosperms to separate between the monocot and eudicot clades. iTOL (iTOL, RRID:SCR_018174) [94] was used to visualize the trees using their Newick format, and root them at P. lambertiana. Moreover, a cladogram representing the known phylogeny of the included taxa was created by hand based on known evolutionary relationships [95–99]. This was used to compare the generated phylogenetic relationship based on functional similarity with the evolutionary relationships of the plant genomes.

Jackknifing analysis was carried out for both parsimony and NJ trees to assess the support for each clade on the basis of the proportion of jackknife trees showing the same clade. To this end, 40% of the terms in T were randomly removed, ancestors of the remaining terms were added, and trees constructed as above. The majority rule consensus tree of 100 individual trees was calculated with the jackknife values represented on each branch. The tree was then visualized using iTOL using its Newick format, and rooted again at P. lambertiana.

Genome quality evaluation

Data Availability

All data and source code generated are freely available at [105] under the terms of the CC0 license (also archived at [106]). All software requirements and dependencies are packaged into a Singularity container (now renamed as Apptainer) so no other set-up is required to reproduce our results (container available at [107]; download and use instructions are included in the repository README).

An up-to-date list of all available annotation sets can be found at [49], and the GOMAP software used to generate them is available at [108].

Additional Files

Supplementary Table S1: Additional assembly statistics from GenomeQC.

Supplementary Table S2: Number of removed annotations during clean-up.

Supplementary Figure S1: Number of total annotations in each GO IEA dataset analyzed, colored by GO aspect (cellular component in green, molecular function in orange, and biological process in blue). Species are ordered on the basis of the number of GO terms in the GOMAP dataset, with the species producing the most GO terms (Triticum) on top and the least (Pinus) on bottom.

Supplementary Figure S2: Quality scores of the predicted annotation visualized as a gray-scale heat map. Each row represents a single species from a single data source. Input gene model dataset origin is indicated: blue for GOMAP, yellow for Gramene63, and orange for Phytozome12. Across the top are the 3 metric types used to assess the annotation quality. Across the bottom are subgraph indicators: C is cellular component, F is molecular function, and P is biological process. Darker cells indicate a higher (better) score whereas lighter cells indicate a lower score. Note that scales are different for each metric type (meaning that comparisons across the 3 metric types are not meaningful). Rows are clustered by pairwise correlation/similarity across all metrics. The dendrogram at left retraces clustering.

Supplementary Figure S3: Neighbor-joining trees built on annotation subsets per GO aspect. Phylograms are colored and rooted as described in Fig. 2. For each tree, only the annotations from the respective aspect of the Gene Ontology were used.

Abbreviations

ADS: artificial dilution series; BOP clade: Bambusoideae, Oryzoideae, and Pooideae; bp: base pairs; BUSCO: Benchmarking Universal Single-Copy Orthologs; CAFA: Critical Assessment of Functional Annotation; CDS: coding sequence; DOI: Digital Object Identifier; EXP evidence code: Inferred from Experiment; GAF: GO Annotation File; GO: Gene Ontology; GOA: Gene Ontology Annotation; GOMAP: Gene Ontology Meta Annotator for Plants; IEA evidence code: Inferred from Electronic Annotation; IEP evidence code: Inferred from Expression Pattern; iTOL: Interactive Tree of Life; IUPAC: International Union of Pure and Applied Chemistry; KEGG: Kyoto Encyclopedia of Genes and Genomes; Mb: megabase pairs; NJ: neighbor joining; NPAAA: non-protein amino acid-accumulating clade; NSF: National Science Foundation; PACMAD clade: Panicoideae, Arundinoideae, Chloridoideae, Micrairoideae, Aristidoideae, and Danthonioideae; PHYLIP: PHYLogeny Inference Package; TC-AUCPCR: term-centric area under precision-recall curve.

Competing Interests

The authors declare that they have no competing interests.

Funding

This work has been supported by the Iowa State University Plant Sciences Institute Faculty Scholars Program to C.J.L.D., the Predictive Plant Phenomics NSF Research Traineeship (No. DGE-1545453) to C.J.L.D. (C.F.Y. and K.O.C. are trainees), the AI Research Institutes program supported by NSF and USDA-NIFA under AI Institute: for Resilient Agriculture (Award No. 2021-67021-35329) to C.J.L.D. (supporting L.F.), and IOW04714 Hatch funding to Iowa State University.

Authors' Contributions

L.F., D.P., C.F.Y., K.O.C., H.A.D., P.J., D.C.S., H.V., and K.W. generated annotations for plants as described in this article. D.P. and C.J.L.D. co-conceived the idea for phylogenetic analysis. D.P. wrote the code for the analyses in this article. L.F. worked with D.P. to create dendrograms and compare those to phylogenetic trees. L.F. carried out assembly and annotation metric comparisons. L.F., D.P., and C.J.L.D. wrote the manuscript. All authors read, offered suggestions to improve, and approved the final copy of the manuscript.

Authors’ Information

K.W. created the GOMAP system during his time as a graduate student at Iowa State University. L.F., D.P., C.F.Y., K.O.C., H.V., and P.J. are currently graduate students. D.P. attended Iowa State University for 2 semesters on a Study Grant from the German-American Fulbright Commission. H.A.D. and D.C.S. are undergraduate students. Each graduate and undergraduate student annotated ≥1 genome over the course of a research rotation lasting ≤1 semester. C.J.L.D. coordinated research activities and manuscript preparation.

ACKNOWLEDGEMENTS

Thanks to Steven Cannon for help to understand phylogenetic relationships among eudicots and helpful discussions and to Toby Kellogg and Jedrzej Szymanski for discussions and ideas on how to consider the data. Thanks to Nancy Manchanda for reviewing documentation and for checking genome versions used as input for GenomeQC. Thanks to Darwin Campbell who guided the deposition of datasets with CyVerse for release and DOI assignment. Thanks to the reviewers for their suggestions which have improved this manuscript.

References

Author notes