Assigning confidence scores to homoeologs using fuzzy logic

Polyploid genomes used

Upland cotton (G. hirsutum; A_tA_tD_tD_t; 2n = 4× = 52) is an important fiber crop, making up 90% of cotton production worldwide (Zhang et al., 2015). The diploid Gossypium progenitors diverged ∼5–10 MYA, hybridized ∼1–2 MYA, and was followed by a genome doubling, giving the allotetraploid cotton (Wendel & Cronn, 2003). The genome size of this allotetraploid is 2.25 Gb, assembled into 26 chromosomes. The “TM-1” G. hirsutum genome, NAU-NBI assembly and annotation v1.1 (Zhang et al., 2015), was obtained from ftp://ftp.bioinfo.wsu.edu/species/Gossypium_hirsutum/NAU-NBI_G.hirsutum_AD1genome/genes/. A total of 66,967 protein-coding genes were used for homoeolog inference.

Oilseed rape (B. napus; 2n = 4× = 38) is an allotetraploid member of the Brassicaceae (mustard/cabbage) family. The progenitors of B. napus diverged about four MYA, followed by a relatively recent hybridization, dating back to 7,500–12,500 years ago (Chalhoub et al., 2014). The 1.13 Gb B. napus genome was sequenced and assembled into 19 pseudomolecules, with 98,130 canonical protein coding genes (Chalhoub et al., 2014). The “Darmor-bzh v5” version of the B. napus genome was obtained from (http://www.genoscope.cns.fr/brassicanapus/data/).

Wheat (T. aestivum) is a staple food. Its genome is large (∼17 Gb) and allohexaploid, consisting of three subgenomes with seven sets of chromosomes each (6× = 2n = 42). The divergence of the progenitor species is estimated to be around seven MYA, and the two hybridization events <1 MYA (International Wheat Genome Sequencing Consortium (IWGSC), 2014). The wheat genome was obtained from Ensembl Plants 33 using the TGACv1 assembly (Clavijo et al., 2017) of the “Chinese Spring” cultivar; ftp://ftp.ensembl.org/pub/release-33/embl/triticum_aestivum/. The annotation consisted of 103,458 genes after removing alternative splice variants. Canonical genes among the splice variants were chosen based on the method described in (Altenhoff et al., 2011).

Homoeolog inference

In OMA, we inferred homoeologs in the polyploid species by treating each subgenome as a separate genome and inferring orthologs following the normal OMA pipeline (Train et al., 2017). Briefly, the steps for obtaining pairwise homoeologs are as follows:

• Subgenome delineation: Polyploid genomes were separated into two or more subgenomes based on their identifiers. From here on, a subgenome is treated as a standalone genome/proteome.

• Homology inference: Using proteomes from all the species in OMA, Smith–Waterman alignments were made with all possible pairs of sequences from all genomes (Smith & Waterman, 1981). Pairs of protein sequences from different genomes (or subgenomes) with sufficient alignment score and overlap were retained.

• Homoeolog and co-homoeolog inference: Pairs retained from the previous step that are the mutually evolutionary closest sequences between a pair of subgenomes are kept during this step. In order to include many-to-many homoeologous relationships, pairs found within a confidence interval of the mutually closest sequences are kept.

• Witness of “non-homoeology” verification: In order to avoid paralogs to be mistakenly identified as homoeologs due to differential gene loss (Dessimoz et al., 2006), a verification step is performed by searching for a third genome that retained both orthologous copies. The third genome can be any of the 2,103 eukaryotic, prokaryotic, or archaeal genomes currently in OMA, but are often very closely related species. This third genome acts as a witnesses of non-homoeology, and pairs that do not pass this test are filtered out.

Recently, several improvements to the OMA algorithm were introduced (Train et al., 2017). This included refinements to account for rapidly evolving duplicated genes, as well as fragmentary sequences. Additionally, those homoeologs which were considered outliers in the previous OMA algorithm (Altenhoff et al., 2015) are now retained.

All data, including the genome information and homoeologous/orthologous relationships, is stored in an HDF5 database (https://omabrowser.org/All.Dec2017/OmaServer.h5), and queried programmatically using python. Our code for the confidence score computations can be found in File S1. Skfuzzy (https://github.com/scikit-fuzzy/scikit-fuzzy) was used for all functions related to fuzzy logic. Only canonical genes were considered (no alternative splice variants) for the following analyses.

Variables used as input for fuzzy sets

Three inputs, that is, features of a given homoeolog pair, were used for the fuzzy logic variables: the evolutionary distance, synteny score, and total copy number.

Synteny is the overall conservation of chromosome order and location of genes when comparing two chromosomes. However, rearrangements may result in smaller regions of the chromosome being syntenic, rather than the whole chromosome. Thus, we computed a local synteny score for each pair of homoeologs (Fig. 1). This consisted of obtaining two windows of genes, one on each subgenome: the first one containing the 10 neighbor genes surrounding one homoeolog, and the second window containing the 10 neighbor genes surrounding the corresponding homoeolog. If the homoeolog was on a scaffold with less than 10 genes, all the genes on the scaffold were used. However, if there were less than two genes (including the homoeolog) on the scaffold, we assigned the synteny score to zero for these pairs. We then computed the proportion of genes (not including the homoeologs) in each window which were homoeologous to the corresponding window on the other subgenome and took the mean of these two as the synteny score.

The evolutionary distance is based on the number of nucleotide substitutions between two sequences. This is in PAM units, and calculated as part of the normal OMA algorithm (Roth, Gonnet & Dessimoz, 2008). A distance of one PAM unit describes an amount of evolution which will change, on average, 1% of the amino acids.

The “total copy number” is a metric to understand the degree of duplication for a pair of homoeologs. For a given pair, it is calculated as the number of homoeologs for the first gene + the number of homoeologs for the second gene.

Defining the universe of discourse and membership functions for all variables

For each genome, input variable (synteny, evolutionary distance, total copy nr), and output variable (confidence score), the universe of discourse is the range of possible values. The universe for each variable was are defined in Table 3, using the “Antecedent” and “Consequent” functions from the skfuzzy control module. Gaussian curve membership functions were defined for each genome using the skfuzzy “gaussmf” function. The membership classes for each fuzzy set are defined in Table 4.

Fuzzy rules and control system simulation

We created five rules based on the universes defined above and stored them in the lookup table (see Results). The control system and simulation were made using the skfuzzy control module and the rules. This simulation was then used to defuzzify, that is, return a crisp output, using the centroid defuzzify method. It takes the inputs and returns a confidence score between 0 and 100. We then kept the smallest confidence score returned as the minimum and scaled the maximum confidence score to be 100. A set of 30 homoeolog pairs were manually evaluated in G. hirsutum. Ten pairs were randomly chosen from the 0–60, 60–90 to 90–100 confidence score ranges. Putative functions were found by searching the protein sequence in the NCBI Conserved Domain Database (CDD) (Marchler-Bauer et al., 2017).

Polyploid species and homoeolog inference

We inferred homoeologs for three polyploid species using OMA: G. hirsutum (GOSHI), B. napus (BRANA), and T. aestivum (WHEAT). In the OMA pipeline, homoeologs are predicted as orthologs between subgenomes, thus, the genes have to be annotated as belonging to a particular subgenome in order to infer homoeology. However, as assemblies can be in different states of draft, not all genes are always mapped to a specific chromosome or subgenome. We therefore discarded 2,002, 515, and 4,151 genes from GOSHI, BRANA, and WHEAT, respectively, as they were not previously mapped to any particular subgenome. This left between 96% and 99.5% of the total genes in each genome to infer homoeology. Using these genes, we predicted pairs of homoeologs between subgenomes for each of the three species. This resulted in 30,230, 35,661, and 88,513 pairs of homoeologs for GOSHI, BRANA, and WHEAT, respectively.

Investigation into global versus local synteny

Rearrangements may result in smaller regions of the chromosome being syntenic, rather than the whole chromosome. In order to justify using a local synteny score rather than a global synteny based on chromosome matching between subgenomes, we computed the number of homoeologs across pairs of chromosomes between two subgenomes in each species. With OMA we predicted many homoeologs across different chromosome groups (Fig. 2).

Non-homoeologous chromosomes, that is, different chromosome groups, with increased frequency of homoeologs are consistent with known reciprocal translocations. For example, in GOSHI, there are two known large reciprocal translocations: parts of the chromosomes were exchanged between A02 and A03, as well as between A04 and A05. We would be able to see this by an increased frequency of homoeolog pairs inferred between chromosomes not belonging to the same chromosome group (Fig. 3). Indeed, with OMA we inferred 692, 1,208, 598, and 696 pairs of homoeologs between chromosomes A02/D03, A03/D02, A04/D05, and A05/D04, respectively (Fig. 2). This is significantly higher than the mean number of homoeolog pairs for non-homoeologous chromosomes. We observe even more levels of chromosomal translocation in WHEAT and BRANA as well, for example, detecting the known reciprocal translocations in WHEAT between 4A, 5A, and 7B (Ma et al., 2013).

Additionally, chromosome pairs with a few number of homoeologs may represent single-gene translocations between non-homoeologous chromosomes and should not be discarded from homoeology prediction via synteny-based methods. Taken together, these results suggest that a local synteny score is more robust than global synteny in order to account for large and small translocations.

Description of fuzzy logic inputs

For each homoeolog pair, we used the synteny score, the evolutionary distance, and the total copy number as input.

The synteny score is the degree of local gene neighborhood conservation. Although synteny is not a hard requirement for homeologs, a conservation of synteny is a good indicator of correct homoeolog predictions. In order to account for chromosomal rearrangements, as well as genome assemblies which are not yet fully assembled into pseudomolecules, we computed a local synteny score. This technique, however, only works when both genes of a homoeolog pair have at least one neighbor gene. Therefore, we could not compute the synteny score for homoeologs that were on small scaffolds with only one gene annotated. For those pairs we set the synteny score to zero. This was 490, 0, and 36,250 pairs for GOSHI, BRANA, and WHEAT, respectively.

The evolutionary distance is based on the number of nucleotide substitutions between two sequences (in PAM units). Because of the relatively short divergence between subgenomes, we expect there to generally be a low distance between homoeologs. Additionally, genes which have a high number of predicted homoeologs could indicate something suspect, such as a transposable element (TE) misannoted as a gene. All distributions of the inputs are shown in Fig. 4.

Fuzzification of data into membership functions

The membership functions allow us to translate an input value into a degree of membership between 0 and 1. The membership curves are overlapping to account for the fuzziness between categories, and for all genomes we defined what we consider to be low, medium, or high input values (Fig. 5). For example, a homoeolog pair with a synteny score of 0.6 could then be translated linguistically to not at all low synteny, somewhat medium synteny, and mostly high synteny. For synteny, the curves are statically defined based on the synteny input universe (possible scores only from 0 to 1), while the distance and total copy number are based on the minimum, maximum, and median values for that particular genome (see Materials and Methods).

The fuzzy rules and the control output

The output from our fuzzy inference process, Confidence, also has a membership function. It is used to map the fuzzy confidence to a crisp confidence score, between 0 and 100 (Fig. 6). Fuzzy logic uses a set of IF THEN statements (rules) for mapping the input space to an output space. We set the rules as to what we think should determine the five categories of confidence (Fig. 6). The very high confidence homoeologs are those which have a high synteny, a low distance, and a low total number of homoeologs. The very low confidence homoeologs have a low synteny score, a high distance, and a high copy number of homoeologs. High, medium, and low confidence is everything in between, based on what we would consider low or high inputs. We put a lot of emphasis on synteny score because the chances of a non-homoeologous pair being syntenic but not originating from a common ancestral gene is very low.

Defuzzification

After defining the rules, we created a control system and simulation for each of the genomes. The inputs for each homoeolog pair were then fed into the simulation which contains the rules, and defuzzified. The defuzzification process converts the fuzzy linguistic confidence to a crisp confidence score, which we then scaled between the minimum value and 100. The reason for scaling to 100 is so that there would not be a sharp cutoff and a maximum confidence score around 80. This facilitates comparison of homoeolog confidence scores within genomes, as people naturally tend to associate the best score with 100. The resulting distribution of confidence scores is shown in Fig. 7. Unscaled confidence scores are shown in Fig. S1.

Evaluation of homoeologs confidence scores

We assessed the homoeolog predictions by looking at the correlation between the total number of orthologs per homoeolog pair and the confidence score. We also manually evaluated a set of 10 homoeolog pairs from 0–60, 60–90, to 90–100 confidence score ranges (Table S1).

The total number of orthologs takes into account the ortholog predictions for all of the species in OMA. Homoeolog pairs with few orthologs are either lineage-specific or dubious, whereas pairs with many orthologs represent those more likely to be true. Although the correlation is low between the total number of orthologs and the confidence score (R = 0.005), it holds as a general trend when looking at confidence scores between 0 and 70 vs. confidence scores between 70 and 100 (Fig. 8). In the set of GOSHI homoeologs manually evaluated, the total number of orthologs was significantly lower for the low confidence homoeologs (0–60), than the 60–90 and 90–100 samples (Table S1).

Interestingly, in GOSHI, for the set of 10 manually evaluated pairs from 0 to 60 confidence, half had RVT-3 (reverse transcriptase-like) domains. According to the CDD description, “This domain is found in plants and appears to be part of a retrotransposon”. This could explain the few number of orthologs for those with low confidence scores, because TEs rapidly evolve and may have lost homology in the other species. By contrast, none of the sampled homoeologs with a confidence score above 60 had RVT-3 annotation, or any functional description associated with TEs.

Finally, in order to compare our new confidence scores to the previous confidence classes in OMA, we looked at the proportion of pairs in each confidence score bin which were previously marked as either high or low confidence (Fig. 9). There are a significant amount of homoeolog pairs which were previously classified as low confidence, but are now between 90 and 100 confidence score (37%, 31%, 5% of the formerly low confidence pairs in GOSHI, BRANA, and WHEAT, respectively). On the other hand, there was a much smaller proportion of pairs that were considered high confidence before, but now have less than 70 as a confidence score (0.5%, 1%, and 4% for GOSHI, BRANA, and WHEAT, respectively). In wheat, the large peak of scores between 70 and 80 that were formerly high confidence is because 68% of these pairs are on small scaffolds for which we were unable to compute synteny; however, they belong to the same chromosome group.