Genomic phylostratography based gene age inference

Discussed intensely in the past years (Capra et al., 2013; Altenhoff et al., 2016; Liebeskind et al., 2016), gene age inference is not a trivial task and might be biased in some currently existing approaches (Liebeskind et al., 2016; Yin et al., 2018; Casola 2018).

In particular, Moyers & Zhang argue that genomic phylostratigraphy (a prominent BLAST based gene age inference method)

  1. underestimates gene age for a considerable fraction of genes,

  2. is biased for rapidly evolving proteins which are short, and/or their most conserved block of sites is small, and

  3. these biases create spurious nonuniform distributions of various gene properties among age groups, many of which cannot be predicted a priori (Moyers & Zhang, 2015; Moyers & Zhang, 2016; Liebeskind et al., 2016).

However, these arguments were based on simulated data and were inconclusive due to errors in their analyses. Furthermore, Domazet-Loso et al., 2017 provide convincing evidence that there is no phylostratigraphic bias. As a response, Moyers & Zhang, 2017 recently published a counter-study stating that a phylostratigraphic trend claimed by Domazet-Loso et al., 2017 to be robust to error disappears when genes likely to be error-resistant are analyzed. Moyers & Zhang, 2017 further suggest a more robust methodology for controlling for the effects of error by first restricting to those genes which can be simulated and then removing those genes which, through simulation, have been shown to be error-prone (see also Moyers & Zhang, 2018).

In general, an objective benchmarking set representing the tree of life is still missing and therefore any procedure aiming to quantify gene ages will be biased to some degree. Based on this debate a recent study suggested to perform gene age inference by combining several common orthology inference algorithms to create gene age datasets and then characterize the error around each age-call on a per-gene and per-algorithm basis. Using this approach systematic error was found to be a large factor in estimating gene age, suggesting that simple consensus algorithms are not enough to give a reliable point estimate (Liebeskind et al., 2016). This was also observed by Moyers & Zhang, 2018) when running alternative tools such as PSIBLAST, HMMER, OMA, etc. However, by generating a consensus gene age and quantifying the possible error in each workflow step, Liebeskind et al., 2016 provide a very useful database of consensus gene ages for a variety of genomes.

Alternatively, Stephen Smith, 2016 argues that de novo gene birth/death and gene family expansion/contraction studies should avoid drawing direct inferences of evolutionary relatedness from measures of sequence similarity alone, and should instead, where possible, use more rigorous phylogeny-based methods. For this purpose, I recommend researchers to consult the phylomedb database to retrieve phylogeny-based gene orthology relationships and use these age estimates in combination with myTAI. Alternatively, users might find the simulation based removal approach proposed by Moyers & Zhang, 2018 more suitable.

Evidently, these advancements in gene age research are very recent and gene age inference is a very young and active field of genomic research. Therefore, many more studies need to address the robust and realistic inference of gene age and a community standard is still missing.

Despite the ongoing debate about how to correctly infer gene age, users of myTAI can perform any gene age inference method they find most appropriate for their biological question and pass this gene age inference table as input to myTAI. To do so, users need to follow the following data format specifications to use their gene age inference table with myTAI. However, even when users rely on established procedures such as phylostratigraphy the gene age inference bias will be present as ‘systematic error’ in all developmental stages for which TAI or TDI computations are performed. Thus, stages of constraint will be detectable in any case. Since TAI or TDI computations are indended to enable screening for conserved or constrained stages in developmental or biological processes for further downstream experimental studies, even simple approaches such as phylostratigraphy can give first evidence for the existence of transcriptomic contraints within a biological process. If researchers then wish to extract the exact candidate genes that might potentially cause such transcriptome constraints then I would advise to rely on more superior approaches of gene age inference as discussed above.

In my opinion, what is completely missing in this entire debate is the bioinformatics/technical aspect of using BLAST or any other BLAST-like tool for gene age inference. Recently, it was intesely discussed how BLAST hits are biased by the use of the default argument max_target_seqs (Shah et al., 2018). The main issue of how this max_target_seqs is set is that:

According to the BLAST documentation itself (2008), this parameter represents the ‘number of aligned sequences to keep’. This statement is commonly interpreted as meaning that BLAST will return the top N database hits for a sequence query if the value of max_target_seqs is set to N. For example, in a recent article (Wang et al., 2016) the authors explicitly state ‘Setting ’max target seqs’ as ‘1’ only the best match result was considered’. To our surprise, we have recently discovered that this intuition is incorrect. Instead, BLAST returns the first N hits that exceed the specified E-value threshold, which may or may not be the highest scoring N hits. The invocation using the parameter ‘-max_target_seqs 1’ simply returns the first good hit found in the database, not the best hit as one would assume. Worse yet, the output produced depends on the order in which the sequences occur in the database. For the same query, different results will be returned by BLAST when using different versions of the database even if all versions contain the same best hit for this database sequence. Even ordering the database in a different way would cause BLAST to return a different ‘top hit’ when setting the max_target_seqs parameter to 1. - Shah et al., 2018

The solution to this issue seems to be that any BLAST search must be performed with a significantly high -max_target_seqs, e.g. -max_target_seqs 10000 (see https://gist.github.com/sujaikumar/504b3b7024eaf3a04ef5 for details) and best hits must be filtered subsequently. It is not clear from any of the studies referenced above how the best BLAST hit was retrieved and which -max_target_seqs values were used to perform BLAST searches in the respective study. Thus, the comparability of the results between studies is impossible and any individual claim made in these studies might be biased.

In addition, the -max_target_seqs argument issue seems not to be the only issue that might influence technical differences in BLAST hit results. Gonzalez-Pech et al., 2018 discuss another problem of retrieving the best BLAST hits based on E-value thresholds.

Many users assume that BLAST alignment hits with E-values less than or equal to the predefined threshold (e.g. 105 via the specification of evalue 1e-5) are identified after the search is completed, in a final step to rank all alignments by E-value, from the smallest (on the top of the list of results) to the largest E-value (at the bottom of the list). However, the E-value filtering step does not occur at the final stage of BLAST; it occurs earlier during the scanning phase (Altschul et al., 1997; Camacho et al., 2009). During this phase, a gapped alignment is generated using less-sensitive heuristic parameters (Camacho et al., 2009); alignments with an E-value that satisfies the defined cut-off are included in the subsequent phase of the BLAST algorithm (and eventually reported). During the final (trace-back) phase, these gapped alignments are further adjusted using moresensitive heuristic parameters (Camacho et al., 2009), and the E-value for each of these refined alignments is then recalculated. - Gonzalez-Pech et al., 2018

This means that if one study mentioned above ran a BLAST search with a BLAST parameter configuration of lets say -max_target_seqs 250 (default value in BLAST) and evalue 10 (default value in BLAST) and then subsequently selected the best hit which returned the smallest E-value and another study used the parameter configuration -max_target_seqs 1 and evalue 0.0001 then the results of both studies would not be comparable and the proposed gene age inference bias might simply result from a technical difference in running BLAST searches.

In more detail, even if one study for example ran BLAST with evalue 10 (default value in BLAST) and then subsequently filtered for hits that resulted in evalue < 0.0001 whereas another study ran BLAST directly with evalue 0.0001, according to Gonzalez-Pech et al., 2018 these studies although referring to the same E-value threshold for filtering hits will result in different sets of filtered BLAST hits.

A recently introduced approach is called synteny-based phylostratigraphy (Arendsee et al., 2019). Here, the authors provide a comparative analysis of genes across evolutionary clades, augmenting standard phylostratigraphy with a detailed, synteny-based analysis. Whereas standard phylostratigraphy searches the proteomes of related species for similarities to focal genes, their fagin pipeline first finds syntenic genomic intervals and then searches within these intervals for any trace of similarity. It searches the (in silico translated) amino acid sequence of all unannotated ORFs as well as all known CDS within the syntenic search space of the target genomes. If no amino acid similarity is found within the syntenic search space, their fagin pipeline will search for nucleotide similarity. Finding nucleotide sequence similarity, but not amino acid similarity, is consistent with a de novo origin of the focal gene. If no similarity of any sort is found, their fagin pipeline will use the syntenic data to infer a possible reason. Thus, they detect indels, scrambled synteny, assembly issues, and regions of uncertain synteny (Arendsee et al., 2019).

Hence, all of the above mentioned approaches are far from being perfect and much more research is needed to systematically compare different approaches for gene age inference.

The overall rational behind gene age inference is to assign each protein coding gene of an organism of interest with an evolutionary age estimate which aims to quantify its potential origin within the tree of life (detectable sequence homolog; orphan gene (see Tautz & Domazet-Loso, 2011)). Hence, gene age inference generates a table storing the gene age in the first column and the corresponding gene id of the organism of iterest in the second column. This table is named phylostratigraphic map.