OrthoFinder: solving fundamental biases in whole genome comparisons dramatically improves orthogroup inference accuracy — TEI-tagged view

David M. EmmspersName, Steven KellypersName (2015date). “OrthoFinder: solving fundamental biases in whole genome comparisons dramatically improves orthogroup inference accuracy.” Genome Biologytitle·j.

Object of study — accurate and scalable inference of multi-species orthogroupsterm from protein sequences despite gene-length and phylogenetic-distance biases

Research question — whether normalised similarity scores and query-specific cognate thresholds improve complete orthogroupterm recovery over existing pairwise and clustering approaches

Contribution — OrthoFinderterm, a five-phase orthogroupterm inference algorithmterm that normalises BLAST scoresterm, builds a cognate-gene graph, clusters it with MCL, and is benchmarked for accuracy, missing-data robustness, scale, and plant transcription-factor discovery

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

OrthoFinderterm: solving fundamental biases in whole genome comparisons dramatically improves orthogroupterm inference accuracy Abstract Identifying homology relationshipsterm between sequences is fundamental to biological research. Here we provide a novel orthogroupterm inference algorithmterm called OrthoFinderterm that solves a previously undetected gene length biasterm in orthogroupterm inference, resulting in significant improvements in accuracy. Using real benchmark datasets we demonstrate that OrthoFinderterm is more accurate than other orthogroupterm inference methods by between 8 % and 33 %. Furthermore, we demonstrate the utility of OrthoFinderterm by providing a complete classification of transcription factorterm gene families in plants revealing 6.9 million previously unobserved relationships.

The online version of this article (doi:10.1186/s13059-015-0721-2) contains supplementary material, which is available to authorized users.

Background and rationale Identifying homology relationshipsterm between sequences is fundamental to all aspects of biological research. In addition to the pivotal role these inferences play in furthering our understanding of the evolution and diversity of life, they also provide a coherent framework for the extrapolation of biological knowledge between organisms. In this context, orthology inferenceterm underpins genome and transcriptome annotationterm and provides the foundation on which synthetic and systems biology is built. Given the importance of this process to biological research there has been a rich heritage of methodology development in this area with the production of several effective orthologyterm databases and algorithmsterm.

The most widely used methods for orthology inferenceterm can be classified into two distinct groups. One group of methods approaches the problem by inferring pairwise relationshipsterm between genes in two species, and then extending orthologyterm to multiple species by identifying sets of genes spanning these species in which each gene-pair is an orthologueterm. Popular methods that adopt this approach include MultiParanoid [ 1 ] and OMA [ 2 ]. A confounding factor to these approaches is that gene duplicationsterm cause orthologyterm relationships that are not one-to-one [ 3 ] and so orthologyterm is not a transitive relationship (for example, if gene A is an orthologueterm of gene B, and gene B is an orthologueterm of gene C, it is not necessarily true that gene A is an orthologueterm of gene C) [ 4 ]. This lack of transitivity means that to capture all pairwise orthologyterm relationships individual genes must be allowed to be members of more than one set [ 2 ], or the gene sets must be restricted to subsets of species that share the same last common ancestor [ 1 ]. Methods that adopt these pairwise approaches have high levels of precisionterm in recovering orthologuesterm, however, they suffer from low rates of recallterm in discovering the complete orthogroupterm due to these complications arising from gene duplicationsterm.

The second group of methods do not adopt this pairwise strategy but rather attempt to identify complete orthogroupsterm; an orthogroupterm is the set of genes that are descended from a single gene in the last common ancestor of all the species being considered [ 2 , 5 – 9 ]. Here an orthogroupterm by definition contains both orthologuesterm and paraloguesterm, and in this context is frequently used as a unit of comparison for comparative genomics [ 10 – 12 ]. In this work we follow this latter approach as it is a logical extension of orthologyterm to multiple species. The most widely used orthogroupterm inference method is OrthoMCLterm [ 13 ] (usage assessed by citations n = 870 Scopus citations at the time of writing this article). OrthoMCLterm uses BLASTterm [ 14 ] to compute sequence similarity scores between sequences in multiple species and then uses the MCL clustering algorithmterm [ 15 ] to identify highly-connected clusters (groups of highly similar sequences) within this dataset.

In addition to the approaches discussed above, several methods have also been developed that incorporate gene syntenyterm/co-linearity information to assist in the inference of orthogroupsterm [ 16 , 17 ]. For groups of organisms such as the Kinetoplastids, where gene syntenyterm/co-linearity is well conserved [ 18 ] it can provide valuable additional information. However, synteny is not conserved over large evolutionary distances and thus can provide little assistance to the identification of related genes between distantly related groups such as plants and metazoa. Moreover, synteny is unavailable for de novo assembled transcriptomes and for fragmented, low-coverage genome assemblies. Thus there is a need to have accurate methods of orthogroupterm inference that do not require gene syntenyterm information.

Here we present OrthoFinderterm, a novel method that infers orthogroupsterm of protein coding genes. It is fast, easy to use and scalable to thousands of genomes. In tests using real benchmark datasets OrthoFinderterm outperforms all other commonly used orthogroupterm inference methods by between 8 % and 33 %. We further demonstrate the utility of OrthoFinderterm through the inference and analysis of plant transcription factor orthogroupsterm. Here we use phylogenetic methods to validate the orthogroupsterm and show that using OrthoFinderterm to infer orthogroupsterm identifies millions of previously unobserved relationships. Further information about the algorithmterm can be found at [ 19 ] and a standalone implementation of the algorithmterm is available under the GPLv3 licence at [ 20 ].

Problem definition, method evaluation and comparison to other approaches Gene length biasterm in BLAST scoresterm affects the accuracy of orthogroupterm detection The inference of orthogroupsterm across multiple species requires a fast method to measure pairwise sequence similarity between all sequences in the species being considered. BLASTterm [ 14 ] is the most widely used method to identify and measure similarity between sequences and thus it underpins the majority of orthologueterm identification methods [ 9 , 13 , 21 – 23 ]. Analysis of the pairwise BLAST scoresterm that are produced when the full set of protein sequences from one species is BLASTterm searched against those from another species revealed that there is a clear length dependency in the scores that are obtained (Fig. 1a and b ). Short sequences cannot produce large bit scoresterm or low e-valuesterm (Fig. 1a and b , respectively), whereas long sequences produce many hits with scores better than those for the best hits of short sequences (Fig. 1a and b ). Thus, methods that construct orthogroupsterm by evaluation of BLAST scoresterm in the absence of gene lengthterm·var information should result in a large number of missing genesterm·var (low recallterm) from orthogroupsterm that contain short genes and a large number of incorrectly clustered genes (low precisionterm) in orthogroupsterm that contain long genes. Fig. 1 Analysis of gene lengthterm·var dependency of BLASTp scores. a BLASTterm log 10 (bit scoreterm) for all hits between Homo sapiens (Homo_sapiens.GRCh37.60.pep.all, 21,841 sequences) and Mus musculus (Mus_musculus.NCBIM37.60.pep.all, 23,111 sequences) . b –log 10 (e-valueterm) for all hits between and Homo sapiens and Mus musculus . To avoid infinite values, BLAST scoresterm of zero have been replaced with the lowest obtainable value 10 −180 . The heat map in both cases goes from blue (lowest density of hits) to red (highest). c The F-scoreterm (red), recallterm (blue) and precisionterm (green) of orthogroupterm inference using OrthoMCLterm plotted as a function of sequence length. The sequences were sorted according to length and divided into four bins with the same number of sequences in each. The F-scoreterm, recallterm and precisionterm were calculated for each bin and the scores plotted against the geometric mean of the length of the sequences in each bin. The error bars show the lower and upper limits of sequence lengths for the shortest and longest sequences in each bin and the dot shows the geometric mean of these lengths. d Histogram of all protein-coding gene lengthsterm·var in Homo sapiens is provided for reference To determine if this was the case we assessed the performance of OrthoMCLterm using the OrthoBenchterm dataset [ 5 ]. OrthoBenchterm is the only publicly available benchmark dataset of manually curated orthogroupsterm. The dataset consists of 70 orthogroupsterm of protein coding genes covering 12 species within the Metazoa where each orthogroupterm contains all the genes derived from a single gene in the last common ancestor of the 12 species considered. For further details concerning the construction, species range and complexity of each orthogroupterm see [ 5 ]. The recallterm and precisionterm of OrthoMCLterm was assessed as a function of gene lengthterm·var in this dataset. This revealed that there were strong dependencies between the performance characteristics of OrthoMCLterm and the length of the gene that was being clustered (Fig. 1c , Additional file 1 : Table S1). Specifically, short sequences suffer from low recallterm rate (that is, many short sequences fail to be assigned to an orthogroupterm) and long sequences suffer from low precisionterm (that is, many long sequences are assigned to the incorrect orthogroupterm) as predicted from the analysis of BLAST scoresterm above. To put these results in perspective the distribution of protein lengths in Homo sapiens is provided in Fig. 1d .

A novel score transformterm eliminates gene length biasterm in orthogroupterm detection Given that orthogroupterm inference shows a clear gene lengthterm·var dependency, we sought to develop a transform of the BLAST scoresterm that would reduce the impact of gene lengthterm·var on clustering accuracy. To do this we developed a novel method that determines the gene lengthterm·var dependency of a given pairwise species comparison from an analysis of the bit scoresterm from an all-versus-all BLASTterm search between the two species. Bit scoresterm were used in place of e-valuesterm as the e-valueterm calculation enforces a limit of 1×10 −180 and thus all scores below this floor are given the same value (that is, 0) (Fig. 1b ) and thus length bias in e-valuesterm is non-uniform and irreversible. As bit scoresterm do not have a threshold value, and they have been previously shown to be capable of facilitating accurate inference of phylogenetic trees [ 24 ], they were selected as the raw data for the development of a novel score transformterm.

In brief, for each species-pair in turn, the all-vs-all BLASTterm hits (Fig. 2a ) were divided into equal sized bins of increasing sequence length according to the product of the query and hit sequence lengths. The top 5 % of hitsterm in each bin (ranked according to BLAST bit scoreterm·var) were used to represent ‘good’ hits for sequences of that length bin between the given species pair (Fig. 2b ). A linear model in log-log space was used to fit a line to these scores using least squares fitting (Fig. 2b ). All of the BLAST bit scoresterm·var that were obtained from each species-pair all-vs-all BLASTterm search are then transformed using this model so that the best hits between sequences in this species pair have equivalent scores that are independent of sequence length (Fig. 2c and d ). Following the transform the poor quality hits for longer sequences were no longer better than the best quality hits for short sequences (Fig. 2c ). This normalisation procedure is applied to each pairwise species comparison independently as the behaviour of the BLAST scoresterm is different for each pairwise species comparison (Additional file 2 : Figure S1). Importantly, this pairwise length normalisation between species also normalises for phylogenetic distanceterm between species (See ‘ Gene length and phylogenetic distance normalisationterm ’ & Additional file 2 : Figure S1). Specifically, the normalisation ensures that the best scoring hits between distantly related species achieve the same scores (on average) to the best scoring hits between closely related species (Additional file 2 : Figure S1). These length and phylogenetic distanceterm normalised scoresterm were then used as the measure of sequence similarity on which all subsequent analysis and clustering were performed. Fig. 2 The gene length and phylogenetic distance normalisationterm procedure for a single species pair. a BLAST bit scoresterm·var for all hits between Homo sapiens and Mus musculus . b BLAST bit scoresterm·var for the top 5 % of BLASTterm hits with least-squares fit of the equation log 10 B qh = a log 10 L qh + b ., where B qh is the bit scoreterm for the hit between sequence q and sequence h and L qh is the product of the gene lengthsterm·var (measured in amino acids). c Gene lengthterm·var and phylogenetic distanceterm normalised BLAST bit scoresterm·var. Note that there are a large number of poor scoring hits for long sequences due to these hits exceeding the BLASTterm search e-valueterm cutoff. d The same top 5 % of BLASTterm hits as shown in (b) after normalisation for reference Application of this novel score transformterm prior to clustering of the OrthoBenchterm dataset resulted in a dramatic reduction in the length dependency of the clustering results (Fig. 3 ). Unlike OrthoMCLterm (Fig. 3a ), neither precisionterm, recallterm nor F-scoreterm displayed any dependency on gene lengthterm·var (Fig. 3b ). Moreover, precisionterm was substantially increased over the entire range of sequence lengths (Fig. 3b ). Fig. 3 Comparison of OrthoFinderterm to other orthogroupterm inference methods. a The length dependency of OrthoMCLterm. b The length dependency of OrthoMCLterm using our normalised similarity scores. c The length dependency of the complete OrthoFinderterm algorithmterm. For A-C scores were calculated as in Fig. 1c . d Comparison of the results of OrthoFinderterm F-scoreterm with all other methods tested in OrthoBenchterm. e As in (d) but for recallterm. f As in (d) but for precisionterm. The error bars show the lower and upper limits of sequence lengths corresponding to the shortest and longest sequences in each bin and the dot shows the geometric mean of these lengths An improved method for orthogroupterm delimitation improves overall accuracy Given that we had reduced gene length biasterm and that precisionterm was high but recallterm was low, we assessed whether a method that could identify a higher proportion of cognate gene-pairs prior to clustering could produce an overall increase in clustering accuracy. Many orthologyterm assignment methods make use of reciprocal best BLAST hitterm (RBH) as it is widely regarded as a high precisionterm method for the identification of orthologuesterm gene-pairs [ 25 – 27 ]. Therefore we also sought to use reciprocal best BLAST hitsterm using our new length-normalised score to assist in construction of the orthogroup graphterm. Henceforth, we refer to a reciprocal best hit that is obtained using the length-normalised score as an RBNHterm (reciprocal best normalised hit).

In brief, for each gene that had successfully identified one or more RBNHsterm, the scores for these RBNHsterm were used to delimit an inclusion thresholdterm (see methods). As all scores are normalised for gene lengthterm·var and phylogenetic distanceterm, hits to other genes (in any species) that had scores above this inclusion thresholdterm were included as putative cognate gene-pairsterm and added to the orthogroup graphterm that was subjected to MCL clustering (for further details see methods). This novel data selection criterion resulted in a dramatic improvement in overall clustering accuracy while maintaining gene lengthterm·var independence (Fig. 3c ). The overall results for OrthoFinderterm, were 0.85 precisionterm, 0.81 recallterm and 0.83 F-scoreterm.

OrthoFinderterm outperforms all other methods from the OrthoBenchterm analysis Given that OrthoFinderterm exhibited high accuracy on the benchmark dataset we sought to determine the relative performance to other commonly used methods for orthogroupterm inference. OrthoFinderterm outperformed all other methods that have been applied to OrthoBenchterm [ 5 ] as measured by F-scoreterm (Fig. 3d ), performing 8 % better than TreeFam (the next best method) 25 % better than OrthoMCLterm (the most widely used method), and 33 % better than OMA (the lowest scoring method in this test). Importantly, the precisionterm and recallterm of OrthoFinderterm were balanced, demonstrating that the method is not biased towards over- or under-clustering of sequences. It should be noted that OMA exhibits a low recallterm in this test as its goal is to identify orthologuesterm instead of complete orthogroupsterm and thus paraloguesterm will be absent from the orthologueterm groups identified by this method. OMA is included here for completeness as it was included in the original OrthoBenchterm analysis [ 5 ].

In addition to accuracy, a number of other criteria were used to compare the performance of the different inference methods in the OrthoBenchterm paper. These criteria included the percentage of orthogroupsterm predicted without any errors, the number of erroneously assigned genes (that is, false positives, and thus also captured by the precisionterm) and missing genesterm·var (that is, false negatives, and thus also captured by recallterm) in the assignment of genes to orthogroupsterm and the proportion of orthogroupsterm affected by these false positive and false negatives. The results for OrthoFinderterm according to these criteria are reported in Additional file 3 : Figure S2 and are consistent with the increased accuracy of OrthoFinderterm compared to other methods. Additionally, the 70 orthogroupsterm that make up the OrthoBenchterm dataset comprise 40 that represent particular biological or technical challenges and 30 randomly chosen orthogroupsterm. Additional file 4 : Figure S3 shows the F-scoresterm for these two categories separately to illustrate the difference in performance of the method for ‘randomly selected’ and ‘difficult’ orthogroupsterm. OrthoFinderterm outperformed all other methods in both categories and achieved an F-scoreterm of 81 % and 90 % on the difficult and randomly selected orthogroupsterm, respectively.

OrthoFinderterm is suitable for the analysis of incomplete datasetsterm As many research groups are producing partial genome assemblies and transcriptome resources it is to be expected that sequence datasets will be missing genesterm·var due to incomplete assembly, low expression or errors in gene prediction. To demonstrate the suitability of OrthoFinderterm for analysing these incomplete datasetsterm we assessed the performance of OrthoFinderterm with between 5 % and 60 % of genes deleted at random from the OrthoBenchterm input sequences. This revealed that the accuracy of OrthoFinderterm is robust to missing data and that it achieved an F-scoreterm of over of 0.75 even when 60 % of the genes were missing from the input dataset (Additional file 5 : Figure S4). Thus OrthoFinderterm is suitable for orthogroupterm inference from partial and incomplete datasetsterm.

OrthoFinderterm is fast and scalable The number of speciesterm·var for which genome or transcriptome sequence resources are available is increasing rapidly and there is a corresponding need to be able to infer orthogroupsterm using these datasets as they emerge. To keep pace with these increasing demands the algorithmterm utilises sparse matricesterm as the central data structure and performs many steps using matrix operationsterm. For example, starting from pre-computed raw BLAST scoresterm the identification of orthogroupsterm for the OrthoBenchterm dataset (12 species, 235,033 sequences) takes 14 min 20 s using OrthoFinderterm on a single core of an Intel Core i7-4770 3.4GHz CPU. For comparison, OrthoMCLterm takes 20 h 12 min to perform the same operation using the same CPU and MySQL for its relational database management system. As the number of genomes that must be analysed increases, the scalability of the methods used becomes increasingly important. To demonstrate the scalability performance of OrthoFinderterm, the full set of sequenced plant genomes from Phytozome version 9.0 (n = 41 [ 28 ]) were clustered and the results are shown in Fig. 4 . Plant genomes were selected for this test as they are large with an average of 30,731 protein coding genes per species in Phytozome version 9.0 and thus they represent a stringent assessment of the scalability of OrthoFinderterm. The memoryterm·var (RAM) requirements increase linearly with the number of speciesterm·var clustered (Fig. 4a ). This is despite the fact that the number of BLASTterm hits increases quadratically with the number of speciesterm·var (Fig. 4c ). This linear scaling is achieved by processing the BLASTterm hits for each species sequentially and independently within OrthoFinderterm. Though the memoryterm·var requirements increase linearly, the time requirementsterm·var starting from pre-computed raw BLAST scoresterm increases quadratically with the number of speciesterm·var (Fig. 4b ). This is to be expected as the number of BLASTterm hits that must be processed also increases quadratically. For example, identifying the orthogroupsterm for all 41 plant species from Phytozome requires approximately 4 GB of RAM and took approximately 3 h on a single CPU core. Fitting the data to a line and extrapolating we estimate that approximately 450 plant sized genomes can be clustered on a linux computer with 64GB of RAM (Fig. 4a ). Thus OrthoFinderterm is capable of large analyses on conventional computing resources. It should be noted here that the BLASTterm searches incur the largest computational cost in any orthogroupterm inference analysis and that this cost is the same for all inference methods that use BLASTterm. In summary OrthoFinderterm is fast and scalable to hundreds of species on conventional computing resources. Fig. 4 Memoryterm·var and time requirementsterm·var of OrthoFinderterm. Sub-samples of between two and 41 plant genomes from Phytozome version 9.0 given pre-calculated BLASTterm results. The average number of genes per species was 30,731. a Maximum RAM requirements. b Time requirementsterm·var. c The number of BLASTterm hits that must be processed for a given number of speciesterm·var (provided to place the time and RAM requirements into context) Inference of high accuracy plant transcription factor orthogroupsterm Given that OrthoFinderterm has increased accuracy over other methods and that gene length biasterm has been eliminated from orthogroupterm inference, we sought to provide an additional demonstration of the utility of OrthoFinderterm for the inference of orthogroupsterm. To do this we selected plant transcription factorsterm as they are short genes and will thus suffer from low rates of recallterm in assignment to orthogroupsterm in the absence of gene length biasterm correction. Moreover transcription factorterm genes are preferentially retained following whole genome duplication eventsterm [ 29 , 30 ] and thus transcription factor orthogroupsterm are larger than average and contain multiple independent duplication events in multiple independent lineages that can cause some inference methods to fail. Finally, previous efforts to define transcription factor orthogroupsterm have utilised OrthoMCLterm [ 31 ]. Thus current transcription factor orthogroupsterm will have low recallterm resulting in fragmented orthogroupsterm spanning few species.

Using established rules for the identification and classification of transcription factorsterm [ 31 ] we identified and typed all of the transcription factorsterm present in the 41 genomes present in Phytozome v9. The complete predicted proteomes from these 41 genomes were then subject to clustering using OrthoFinderterm and OrthoMCLterm and the distribution of transcription factorsterm in the resultant orthogroupsterm were analysed. OrthoMCLterm was used here as it is the method by which all transcription factorterm families are currently classified [ 31 ]. Consistent with the increased recallterm rate for OrthoFinderterm, analysis of the resulting orthogroupsterm revealed that 8.5 % more transcription factorsterm were placed in orthogroupsterm using OrthoFinderterm than OrthoMCLterm (Fig. 5a , 97.6 % and 89.1 %, respectively, n = 52,744). Also consistent with the increased recallterm rate is that these orthogroupsterm were less fragmented than those that were produced by OrthoMCLterm (Fig. 5b , 897 and 3,024 orthogroupsterm, respectively). Importantly, the orthogroupsterm inferred using OrthoFinderterm were missing fewer RBHs (Fig. 5c , 2.15 % and 5.77 %, respectively) and clustered more of the same type of transcription factorterm together (Fig. 5d and e ). A major difference between those orthogroupsterm inferred using OrthoFinderterm and OrthoMCLterm is that those produced by OrthoFinderterm encompass a larger number of speciesterm·var than those recovered by OrthoMCLterm (Fig. 5f ), thus orthogroupsterm produced by OrthoFinderterm encompass greater phylogenetic distancesterm. Fig. 5 Inference of orthogroupsterm of plant transcription factorsterm. In all cases dark grey bars indicate the results for OrthoFinderterm and light grey bars indicate the results for OrthoMCLterm. a Comparison of the fraction of transcription factorsterm that are assigned to orthogroupsterm by OrthoFinderterm and by OrthoMCLterm. b Comparison of the number of transcription factor orthogroupsterm identified using each method. c The percentage of RBNHterm/RBH (for OrthoFinderterm/OrthoMCLterm) hits that are not contained in orthogroupsterm identified using each method. d The number of transcription factorsterm of the same type that each transcription factorterm is connected to in the orthogroupsterm produced by OrthoFinderterm. e as in (d) but for OrthoMCLterm. f Comparison of species coverage for transcription factor orthogroupsterm identified by each method. g The number of orthogroupsterm for each transcription factorterm type identified by OrthoFinderterm As OrthoFinderterm clustered the transcription factorsterm together into far fewer orthogroupsterm than OrthoMCLterm (897 versus 3024) we sought to demonstrate that it was correct in doing so. To do this we used gene-tree/species-tree reconciliationterm to determine if the orthogroupsterm were true orthogroupsterm·def if they incorrectly clustered sequences that are separated by a gene duplicationterm event that occurred before the last common ancestor of the species in the analysis. Overall, 858 of the 897 OrthoFinderterm orthogroupsterm (96 %) consisted entirely of genes that were correctly clustered together and only 39 contained some genes that were separated by a duplication prior to the last common ancestor (Additional file 6 : Table S2 and Additional file 7 : Table S3). Of the 897 OrthoFinderterm orthogroupsterm, 210 were identical to ones from OrthoMCLterm and 471 OrthoFinderterm orthogroupsterm were strict supersets of 2,271 OrthoMCLterm orthogroupsterm (Additional file 6 : Table S2 and Additional file 7 : Table S3). Of these, 90 % (425) were true orthogroupsterm·def that each encompassed on average four OrthoMCLterm orthogroupsterm (1,709 in total).

An illustrated example showing an OrthoFinderterm orthogroupterm and its constituent OrthoMCLterm orthogroupsterm is provided in Fig. 6 . Here the OrthoFinderterm orthogroupterm (labelled bHLH 8 in Additional file 6 : Table S2) contains all known type IVc bHLH transcription factorsterm [ 32 ]. Type IVc bHLH transcription factorsterm have previously been shown to be conserved from green algae to land plants and thus span the complete taxonomic range contained in this analysis [ 32 ]. The OrthoFinderterm orthogroupterm correctly united eight paraphyletic OrthoMCLterm orthogroupsterm and included 36 transcription factorsterm (highlighted in grey) that were not clustered into any orthogroupsterm by OrthoMCLterm (Fig. 6 ). The phylogenetic tree shows that there are no genes present in this OrthoFinderterm orthogroupterm that were the product of a gene duplicationterm event prior to the divergence of the last common ancestor of all species in the analysis. This is only one example and the complete set of phylogenetic trees for each OrthoFinderterm transcription factor orthogroupterm·def are provided in Additional file 6 : Table S2 along with the OrthoMCLterm subsets that comprise these groups where appropriate. Also contained in this table are the results of the gene-tree/species-tree reconciliationterm for each tree inferred from an OrthoFinderterm orthogroupterm. Fig. 6 A bootstrapped maximum likelihood phylogenetic tree of the OrthoFinderterm orthogroupterm containing the type IVc bHLH transcription factorsterm (bHLH 8). The OrthoMCLterm orthogroupsterm that are subsets of the OrthoFinderterm orthogroupterm are indicated by different coloured fonts. Thirty-six of the OrthoFinderterm clustered genes (coloured grey) failed to be clustered in any OrthoMCLterm orthogroupterm. The tree was inferred using RAxML using the PROTGAMMAAUTO model (the JTT was model was selected as having the highest likelihood) with 100 bootstrap replicates. Scale bar indicates the number of substitutions per site. Percentage bootstrap support values are indicated by coloured circles shown at internal nodes Taken together, using OrthoFinderterm to cluster transcription factorterm genes resulted in the identification 687 (897 less the 210 that were the same) novel orthogroupsterm of transcription factorsterm across 41 different species comprising 7.7 million pairwise relationshipsterm (of which 6.9 million are not detected by OrthoMCLterm). Thus using OrthoFinderterm to cluster transcription factorsterm has provided significant new insight into the relationship of transcription factorterm genes across plants. The number of orthogroupsterm for each transcription factorterm type is provided in Fig. 5g and the full classification including all constituent accession numbers is provided in Additional file 6 : Table S2.

Algorithmterm implementation and evaluation criteria OrthoFinderterm is an algorithmterm that infers orthogroupsterm across multiple species. The method does not classify the pairwise relationshipsterm that exist between genes within these orthogroupsterm. The method does not require synteny information and is thus equally suitable for clustering protein sequences predicted from genome or transcriptome datasets. OrthoFinderterm is run with a single commandterm and requires as input a directory containing one protein sequence FASTA fileterm per species to be clustered. OrthoFinderterm does not require preprocessing of FASTA files (such as filtering of sequences) and does not require knowledge or use of any relational database management system such as MySQL. It outputs orthogroupsterm in two file formats: the Quest for Orthologs community standard OrthoXMLterm [ 33 ] and in plain text format with one orthogroupterm per line.

There are two common problem definitions used by the majority of homology inference algorithmsterm. One is to predict pairs of orthologuesterm (pairs of genes from two different species descendent from a single gene in the last common ancestor of the two species) and pairs of recent, within-species paraloguesterm·def (genes-pairs arising from a duplication event since the last speciation event for that species). The other approach, and the one used here for OrthoFinderterm, is to predict orthogroupsterm. An orthogroupterm is the set of genes derived from a single gene in the last common ancestorterm of all the species under consideration. This is the approach used by OrthoMCLterm [ 13 ] and eggNOG [ 34 ]. OrthoFinderterm follows this second approach to produce orthogroupsterm of protein coding genes as this is a logical extension of orthologyterm to multiple species as it groups all genes descended from a single gene in the last common ancestor of all species being considered.

Methodological overview of the OrthoFinderterm algorithmterm An overview of the algorithmterm in shown in Fig. 7 , it proceeds in five phases corresponding to sections b-f in the figure: Fig. 7 Overview of the steps in the OrthoFinderterm algorithmterm for two example orthogroupsterm of genes from three species. a The unknown orthogroupsterm that the algorithmterm must recover, shown as a gene tree. b BLASTterm search of all genes against all genes. c Gene length and phylogenetic distance normalisationterm of BLAST bit scoresterm·var to give the scores to be used for orthogroupterm inference. d Selection of putative cognate gene-pairsterm from normalised BLAST scoresterm. e Construction of orthogroup graphterm, genes are nodes in the graph and pairs of genes are connected by an edge with edge weights given by the normalised bit scoreterm. f Clustering of genes into discrete orthogroupsterm using MCL BLASTterm all-versus-all search (Fig. 7b ). Protein BLASTterm (blastp) with an e-valueterm threshold of 10 −3 is used so as to avoid discarding putative good hits for very short sequences. A relaxed threshold is used at this stage of the method as subsequent steps filter out false positive hits using stringent, orthogroup-specific criteria for inclusion (described below). Gene length and phylogenetic distance normalisationterm of the BLAST bit scoresterm·var (Fig. 7c ). This step models the all-vs-all BLASTterm hits for each pairwise comparison between species to identify and remove the gene similarity dependency on gene lengthterm·var and phylogenetic distanceterm. This is done so that the best hits between all species achieve the same scores regardless of sequence length or phylogenetic distanceterm. Delimitation of orthogroupterm sequence similarity thresholds using RBNHsterm (Fig. 7d ). This step uses information from RBNHsterm (Reciprocal Best length-Normalised hit) to define the lower limit of sequence similarity for putative cognate genes of each query gene. To be included in the orthogroup graphterm a gene-pair must be an RBNHterm or produce a hit that is better scoring than the lowest scoring RBNHterm (irrespective of species) for either gene. Constructing an orthogroup graphterm for input into MCL (Fig. 7e ). Putative cognate gene-pairsterm are identified as above and are connected in the orthogroup graphterm with weights given by the normalised BLAST bit scoresterm·var. Clustering of genes into orthogroupsterm using MCL (Fig. 7f ) .

The steps 2 to 4 are the novel parts of our algorithmterm and are described in detail below.

Gene length and phylogenetic distance normalisationterm The aim of this normalisation procedure is to remove gene length biasterm from BLAST bit scoresterm·var and to normalise for phylogenetic distanceterm between species. MCL converts sets of similarity scores into clusters by breaking apart clusters of genes that have low similarity scores (and therefore are unlikely to be orthogroupsterm) and preserving clusters of sequences that have high similarity scores. If the similarity scores between long sequences are inherently larger than the similarity scores between short sequences then the clustering will preferentially break apart clusters of short sequences while preserving clusters of long sequences. This effect can be clearly seen in the results of a typical OrthoMCLterm cluster. Here, long sequences are placed in overly large clusters leading to low precisionterm, and short sequences remain un-clustered leading to low recallterm (Fig. 3a ). The species-wise normalisation implemented by OrthoFinderterm similarly ensures that orthologuesterm from more distant species (that have inherently lower similarity scores due to phylogenetic distanceterm) are not preferentially discarded and is similar to a step that is performed in OrthoMCLterm wherein all scores are divided by the average score between a given species pair [ 13 ].

Previous methods have exploited BLASTterm e-valuesterm (rather than bit scoresterm) as a measure of similarity between sequences. However, as can be seen in Fig. 1b the use of e-valuesterm for assessment of similarity between sequences is flawed. Here, the minimum e-valueterm that can be obtained for a given query sequence decreases with increasing sequence length until, at a certain length, the lower bound for e-valuesterm is reached and BLASTterm returns an e-valueterm of 0. This creates two problems: (1) long sequences will frequently have low quality hits with better e-valuesterm than the best possible hits of short sequences; and (2) the floor value for the e-valueterm calculation means that length bias is non-uniform and thus irreversible. Specifically, beyond the floor-value e-valuesterm cannot be used to distinguish between the qualities of hits as they are all assigned the same e-valueterm. As can be seen in the heat map shown in Fig. 1b , many hits obtain this floor-value for a given pairwise species comparison and thus their similarities are indistinguishable. This length-bias must be removed to prevent biasing downstream clustering applications.

In this method we construct a similarity measureterm between sequences based on the bit-score normalised to take into account the query and hit sequences lengths and the phylogenetic distanceterm between species. Unlike e-valuesterm, the bit-scores do not suffer from the presence of a threshold limit and thus different amounts of sequence similarity can be distinguished regardless of the lengths of the sequences involved. Let L q be the length of the query sequence and L h be the length of the hit sequence. In an analogous manner to the e-valueterm calculation made by BLASTterm and other sequence comparison methods, we use the variable L qh = L q L h to quantify the lengths of a pair of sequences that are being compared.

The length normalisation procedure is shown in Fig. 2 . For each species pair, we: Sort all BLASTterm hits according to L qh = L q L h . Put the hits into equal sized bins of 1,000 hits (put the ‘shortest’ 1,000 hits according to L qh into one bin, the next 1,000 hits into the next bin and so on for all the hits). If there are fewer than 5,000 hits then we put the hits into bins of 200. Using fixed sized bins means that it is not necessary for the algorithmterm to specify the location of the bins in advance. Sort the hits in each bin according to BLAST bit scoreterm·var and select the top 5 % of hitsterm from each bin. Find the parameters a and b that best describe the fit between sequence similarity scores and gene lengthterm·var for the selected hits using the equation log 10 B qh = a log 10 L qh + b where B qh is the BLAST bit scoreterm·var between sequences q and h. Normalise all obtained BLAST bit scoresterm·var (not just the top 5 %) between the given species pair according to, B qh ' = B qh /10 b L qh a , so that B ' qh , (the normalised scoreterm) is the BLAST bit scoreterm·var for a hit divided by the BLAST bit scoreterm·var that would be expected for the best hits between sequences of that length for the species pair under consideration.

The top 5 % of hitsterm are used rather than RBHs as selection of RBHs will be affected by the gene length-bias that we wish to correct. Moreover, gene duplicationterm events can frequently cause RBHs to fail (Additional file 8 : Figure S5) and thus reduce the number of data points that are available for fitting. The normalisation procedure ensures that the best hits between a given species pair achieve (on average) the same scores irrespective of their gene lengthterm·var.

OrthoFinderterm also normalises for phylogenetic distanceterm, this is done so that the similarity scores between orthologuesterm will be independent of phylogenetic distanceterm (that is, the true orthologuesterm in distantly related species will obtain similar scores to the true orthologuesterm in closely related species). If this step is not done then true orthologuesterm in distantly related species will always obtain lower scores than true orthologuesterm in closely related species. Thus during graph clustering (which is unaware of phylogenetic relationship between species) distantly related true orthologuesterm (and cognates) will become disconnected from each other more easily than closely related true orthologuesterm (and cognates) in the orthogroup graphterm. Previous efforts to prevent this phylogenetic bias include dividing the observed similarity score for any given gene-pair by the mean similarity score observed for all reciprocal best hits between genes in that species pair [ 13 ]. However, in the absence of gene lengthterm·var information this means that short genes will always be penalised more than long genes.

Though there is precedent for the use of L qh = L q L h to quantify the lengths of a pair of sequences that are being compared [ 14 ], different functions for gene lengthterm·var normalisation were also assessed. All other functions, including for example the use of the variable \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\tilde{L}}_{qh}={L}_q+{L}_h $$\end{document} L ˜ q h = L q + L h in place of L qh , gave a lower overall clustering accuracy.

Identification of putative cognate gene-pairsterm for inclusion in the orthogroup graphterm Once scores are normalised OrthoFinderterm exploits RBNHsterm to identify putative cognate gene-pairsterm. RBHs are a high precisionterm method to identify putative orthologuesterm [ 25 – 27 ] and OrthoFinderterm uses the reciprocal requirement exploiting its length and phylogenetic distanceterm normalised BLAST scoresterm. For each gene the scores for its RBNHsterm are used to delimit the extent of sequence similarity of that gene’s orthogroupterm. Specifically, for each query sequence, q, any hit, h, with a normalised scoreterm, B ' qh , greater than or equal to the score for the lowest scoring RBNHterm of q is selected as a putative cognate gene-pairterm·def of q and therefore is connected to q in the orthogroup graphterm that is subsequently subjected to MCL clustering.

The rationale for this approach is that the level of normalised similarity of a query gene and its RBNHsterm can be used to estimate the extent of similarity of other genes within the same orthogroupterm. All genes more similar to a query gene than any of the query gene’s RBNHsterm (irrespective of species) are likely members of the same orthogroupterm. Therefore, the normalised similarity score for the most dissimilar RBNHterm of a gene is used as a cutoff for inclusion of additional cognate gene-pairs from all species. That is q is connected to h in the orthogroup graphterm if B ' qh > B ' qR where R is an RBNHterm of q. This provides a simple and robust method for recovering cognate gene-pairs that may otherwise be difficult to detect due to duplication events that can cause the RBNHterm method to fail. Further details, explanation and worked examples are provided in Additional file 8 : Figure S5.

In summary, the novel method presented here generates, for each query gene, an independent prediction of all the genes in its orthogroupterm. This orthogroup graphterm is then clustered using MCL with its default inflation parameter of 1.5. The effect of varying the MCL inflation parameterterm on the OrthoFinderterm result is shown in Additional file 9 : Figure S6. The F-scoreterm of OrthoFinderterm is relatively stable to variation in MCL inflation parameterterm, however it is possible to trade precisionterm against recallterm by varying this parameter (Additional file 9 : Figure S6). For comparison the analogous analysis is also presented for OrhtoMCL (Additional file 10 : Figure S7).

Implementation OrthoFinderterm is written in python. It requires python together with the numpy and scipy libraries [ 35 ] to be installed. OrthoFinderterm requires the standalone BLASTterm+ and MCL algorithmsterm that are freely available. These standalone applications must be installed separately to OrthoFinderterm and are not included in the OrthoFinderterm package. The implementation makes use of sparse matricesterm to store hits between sequences. This provides a memoryterm·var efficient method of storing the data and allows key parts of the algorithmterm to be expressed using scipy’s highly optimised C++ implementations of sparse matrix operationsterm. OrthoFinderterm can either run the BLASTterm searches for you or can be run on pre-computed BLASTterm searches. If you chose to run BLASTterm searches independently then instructions are provided in the documentation for how to process your sequence names in the pre-computed BLASTterm output. Similarly OrthoFinderterm will also automatically run MCL for you. However if you wish to run MCL separately using different parameter settings then the MCL input files are stored for this purpose in a working directory.

Evaluation OrthoBenchterm [ 5 ] is the only manually curated dataset of orthogroupsterm for the assessment of orthogroupterm prediction algorithmsterm. It was used in this work for assessing OrthoFinderterm as it has been independently evaluated, it underpins the testing of multiple different methods and it is a well-defined and stringent test of the problem that OrthoFinderterm was designed to solve. Criteria such as functional similarity within orthogroupsterm, expressed for example using enzyme classification numbers [ 36 ], were not used in this work since not all proteins with the same function are members of the same orthogroupterm and members of the same orthogroupterm do not necessarily all have the same function. As we are using real benchmark datasets for which only a subset of sequences have been assigned to ‘true’ gene families the extent of true negative orthologue assignmentsterm is unknown (as is the case for all methods tested on this dataset). Thus we cannot use the Matthews correlation coefficient to assess the performance of the orthogroupterm inference methods. In the absence of this information the simplest and most transparent evaluation of the accuracy of any prediction method is to measure its precisionterm and recallterm. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ precisionterm=\frac{TP}{TP+FP} $$\end{document} precisionterm = T P T P + F P \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ recallterm=\frac{TP}{TP+FN} $$\end{document} recallterm = T P T P + F N Where TP is the number of true positive orthogroupterm assignments (that is, correct assignments), FP is the number of false positive orthogroupterm assignments (that is, incorrect assignments) and FN is the number of false negative orthologueterm assignments (that is, missing assignments). The F-scoreterm is the harmonic meanterm of these two measures, where the harmonic meanterm weights towards the worst performing measure. We also provide other evaluation measures from the original OrthoBenchterm analysis in Additional file 3 : Figure S2.

Inference of transcription factor orthogroupsterm To infer transcription factor orthogroupsterm we first identified the set of transcription factorsterm present in all genomes present in Phytozome V9. This identification was performed using the same rules for the presence and absence of PFAM domains as has been previously described [ 31 ]. The full set of protein coding genes from these genomes (including all the transcription factorsterm) was then clustered using OrthoFinderterm and OrthoMCLterm and the distribution of the transcription factorsterm within these orthogroupsterm was analysed. OrthoMCLterm was selected for comparison here as it is the method by which all orthogroupsterm of transcription factorsterm are currently defined [ 31 ]. An orthogroupterm of transcription factorsterm was defined as an orthogroupterm whose constituent genes comprised ≥50 % transcription factorsterm of the same domain classification.

To determine if OrthoFinderterm was correct in combining multiple separate OrthoMCLterm orthogroupsterm each orthogroupterm was subject to gene-tree—species-tree reconciliation. Using, gene-tree species-tree reconciliation it is possible to determine if OrthoFinderterm had incorrectly placed together genes that had diverged prior to the last common ancestor of the species being analysed. To do this, gene trees were inferred for each orthogroupterm by aligning the sequences using mafft-linsi [ 37 ] and inferring a maximum likelihood tree from this alignment using FastTree [ 38 ]. DLCpar [ 39 ] was used to reconcile these gene trees with the known species tree [ 28 ]. Using this method, each gene tree was assessed to determine if it contained bipartitions that occurred prior to the divergence of the last common ancestor of all the species being analysed. If such a bipartition was identified then the orthogroupterm was considered not to be a true orthogroupterm·def as it contained one or more genes that evolved by duplication prior to the last common ancestor of all species under consideration.

Discussion In this work we have presented OrthoFinderterm, a novel method for inference of orthogroupsterm. Our method is focused on a clear definition of an orthogroupterm, that is, that an orthogroupterm contains all genes descended from a single gene in the last common ancestor of the species whose genes are being analysed. This definition avoids conflating shared ancestryterm with other criteria that are not equivalent, such as functional conservationterm. Our method is designed to address the problem of orthogroupterm inference rather than categorise the disparate relationships that occur between individual genes within an orthogroupterm. These relationships are best resolved by first inferring orthogroupsterm using OrthoFinderterm and then using multiple sequence alignment and phylogenetic methods on these orthogroupsterm.

The two key novel features of our method are: (1) a method to automatically remove gene length biasterm and phylogenetic distanceterm from sequence similarity scores; and (2) a novel method to define the sequence similarity limits of an orthogroupterm. In the tests performed on the only publicly available orthogroupterm benchmark dataset (OrthoBenchterm) OrthoFinderterm outperformed all of the commonly used orthogroupterm assignment methods by between 8 % and 33 %. Moreover we have shown OrthoFinderterm to be scalable and robust to missing genesterm·var typical of incomplete genomes and de novo transcriptome assemblies. The software is freely available and can take pre-computed BLAST scoresterm as input making it easy to test on any newly developed benchmarks for which pre-computed BLAST scoresterm are available.

We further demonstrate the utility of OrthoFinderterm by providing a novel classification of all transcription factorsterm in the available, fully-sequenced plant genomes present in Phytozome V9. This analysis clusters 97.6 % of the 52,744 putative transcription factorsterm into orthogroupsterm. This novel analysis identifies millions of relationships that have not previously been reported providing new insight into the relationship and evolution of transcription factorterm gene families in plants.

Inferring orthologuesterm underpins much of modern biological research and is among the first steps in the annotation and analysis of genome and transcriptome sequencing projects. As sequencing technologies are now within the budgets of most research groups these data resources are increasing in number rapidly. Thus there is a requirement for an orthogroupterm inference method that is accurate, robust, scalable, and that can be run easily by independent research groups on conventional computing resources. Many orthogroupterm inference methods are not available for general use but are provided as static databases (for example, EggNog and TreeFam). Thus the most widely used methods are those that enable researchers to analyse their own data resources. With this in mind OrthoFinderterm has been developed with the aim of being easy to use. The method is executed as a single commandterm, has minimal dependencies and requires as input just the individual protein sequence FASTA filesterm for each species that is being clustered. The algorithmterm carries out all calculations (including BLASTterm searches and MCL clustering) and outputs the orthogroupsterm in both a plain tab delimited text file and in the OrthoXMLterm community format. The algorithmterm itself is small, fast and memoryterm·var efficient, making it suitable for use on linux desktop computers. Further information about the algorithmterm can be found at [ 19 ] and a standalone implementation of the algorithmterm is available under the GPLv3 licence at [ 20 ].