Sequence retrieval and quality checks

We retrieved 132 publicly available datasets (coding DNA sequence and protein-coding genes) for 46 (Supplementary Table S1) and assessed their gene annotation predictions using BUSCO Plants v3.0.2 (embryophita_odb10) (20). We checked that the number of CDS was consistent with the number of proteins and that they shared the same locus tag name. When protein-coding genes were missing, we generated the sequences from GFF files originating from data providers. Finally, alternate splices were filtered, and the longest sequence was conserved.

Pangene construction

Out of the full dataset, 105 genomes were considered to produce 19 pangenes computed with the get_homologues-est software v20092018 (21) (Table 1), based on NCBI Blast-v2.2 and using the following program options (-M -F -t 0 -m cluster) as previously applied to the Brachypodium distachyon dataset (13). We processed each cluster the following way:

For single-gene copy clusters (a single sequence per genome), protein sequences were aligned using MAFFT v7.313 (22) (parameters adjusted according to the number of sequences) and an automatic procedure to generate a consensus sequence was applied (Figure 1A). For each position of the alignment, we kept the most frequent amino acid. In case of a tie, the amino acid of the genome with the highest BUSCO scores (‘complete’ then ‘fragmented’ and finally lowest ‘missing’ scores) was selected. Finally, if more gaps than amino acids were present, this position was removed from the sequence.

For multi-copy clusters (multiple sequences per genome), we applied the same procedure as for single-gene copy clusters but added a preliminary step to select a representative sequence by cluster. Multiple sequence alignments were used to generate a distance matrix using distmat (Jukes-Cantor correction method) from EMBOSS v6.6 (Figure 1B). The matrix was required to define the distance for each sequence of all other genomes and the sequences with the smallest sum of distance was selected as representative of the considered genome (Figure 1C). Then, the consensus step was applied. It is worth mentioning that get_homologs-est generated sequence clusters of not too distantly related sequences. Large gene families can include several multi-copy clusters being grouped together at the sequence clustering step.

For genotype-specific clusters (paralogs in a single genome), we generated a distance matrix between all sequences and the sequence with the lowest average distance (min(d/sum(d))) of all sequences was putatively considered as the most representative sequence. Those sequences were added to the pangene.

Finally, singletons (cluster of one sequence) were searched for similarity using DIAMOND (23) with a default e-value on the protein-coding genes of all other species genomes to predict their putative prediction accuracy. Sequences with a minimum of one hit in at least two species were added to the pangene; otherwise sequences were excluded.

Table 1.

List of GreenPhylDB pangenes with associated statistics

Species Code	Genomes available by species name	# Genes	Busco complete (%)	# Pangenes	% singletons (not in pangene)
BRADI	Brachypodium distachyon (54 genomes)	44 858 (average)	95.5 (average)	61 622	0.8
BRANA	Brassica napus	101 040	98.3	77 456	5.0
		80 382	96.2		0.5
		70 162	91.2		0.6
BRAOL	Brassica oleracea	35 400	80.7	60 869	2.4
		61 279	96.9		7.7
		56 687	99.5		3.3
BRARR	Brassica rapa	46 250	98.1	49 916	4.4
		46 721	97		4.1
CAPAN	Capsicum annuum	35 336	90.6	41 828	7.4
		34 476	84.6		6.0
		35 884	88.9		8.2
CICAR	Cicer arietinum	28 269	94.8	25 013	7.5
		30 257	93		11.7
COCNU	Cocos nucifiera	52 931	87.8	38 584	1.4
		34 953	86.3		7.0
CUCSA	Cucumis sativus	22 324	89.5	23 446	11.2
		23 780	94.8		7.5
		22 935	94.8		5.7
IPOTF	Ipomoea trifida	32 301	96.6	21 417	8.0
		30 227	94.7		5.3
MAIZE	Zea mays	39 591	94.6	45 301	4.6
		40 003	92.6		6.7
		40 557	87.3		9.4
		36 509	87.8		4.9
MALDO	Malus domestica	45 116	98.3	54 987	7.5
		95 232	91		24.5
		44 677	95.2		15.8
MEDTR	Medicago truncatula	50 444	96.6	43 859	19.2
		44 623	98.8		4.7
MUSAC	Musa acuminata	35 276	98.5	45905	5.5
		44 702	60.3		17.5
		32 692	71.2		18.2
		45 069	71.9		21.0
ORYSA	Oryza sativa	55 986	95.1	56785	11.5
		36 140	87.8		11.3
		37 549	96.1		11.1
		60 897	93.7		9.9
		60 123	89.4		10.7
		35 495	89.9		5.1
		35 594	99.3		2.3
SORBI	Sorghum bicolor	34 129	99.2	45054	11.8
		36 110	97.4		15.4
SOYBN	Glycine max	54 175	99.5	34512	5.0
		52 130	99.3		4.1
THECC	Theobroma cacao	21 330	99.2	32917	1.9
		44 607	99.6		23.1
TRITU	Triticum turgidum	107 891	99.6	54687	8.6
		67 182	98.9		5.3
VITVI	Vitis vinifera	41 733	98.4	43766	22.8
		96 331	92.4		5.1
		73 109	95.9		5.6

Figure 1.

Pangene construction. (A) Schema illustrating the creation of consensus sequences. (B) Example of distance matrix network with five sequences from three genomes of Brassica napus (brana_pan_p029014). (C) Selection of representative sequences based on the minimum value of summed genetic distances. Sequences conserved are in bold.

As a unique identifier was required for each pangene, we defined a nomenclature with a prefix composed of the [5-letter UniProt taxonomy database code]_pan, followed by p (for protein) and an auto-increment of 6-digits (e.g. musac_pan_p029014 for Musa acuminata).

Sequence clustering and functional annotation

Pangenes and protein-coding genes of reference genomes (without pangenes) were searched all against all using DIAMOND. We then performed a clustering using TribeMCL (24) (M = 1.2, 2, 3 and 5), defining 4 levels of stringency (from 1 to 4) to take into account potential sub-classification and we obtained 9419, 18 805, 23 409 and 29 345 clusters, respectively.

We then scanned all sequences for protein domain signatures using InterProscan (25,26) and also crossed linked matches with UniProtKB-SwissProt entries (27). Cluster names resulting from curation from previous GreenPhylDB versions (2,28) were transferred when at least 51% of sequences were found clustered together as before (based on species in common between releases). In addition, for this release, we implemented an automatic method to name clusters based on the name of InterPro domains (family type only) that were found specific to clusters. In other words, when detected in at least 51% of the sequences composing an unannotated cluster, the name of the InterPro signature was assigned to it. In total, GreenPhylDB comprises 3538 clusters functionally characterized across the four levels.

Homology inference

The previous phylogenetic-based methodology that we applied in the previous version has been conserved but uses a larger set of genomes to update our automated pipeline. The pipeline uses MAFFT for the multiple alignment step. FastTree 2 (v2.1.11) (29) was preferred over PhyML (30) due to the size of the clusters. Gene rooting and orthologous scoring was computed with Rap-Green (31) using the viridiplantae species tree extracted from NCBI taxonomy and converted into PhyloXML (2). We successfully produced gene trees at level 1 for more than 99.8% of the clusters (n = 9413) which enabled us to predict ∼17.8 million of orthologs and ∼1.8 million of in-paralogs (or ultra-paralogs) relationships. The pipeline was complemented by a Reciprocal Best Hits (RBH) method—computed between all pairs of genomes—that resulted in more than ∼12.1 million orthologous relationships.

Gene family pages

All cluster (or gene family) pages present the same type of information divided into several tabs :

Gene family composition: a bar chart allows users to visualize at a glance the composition of the gene family by species (Figure 2A). Species are ordered taxonomically to easily detect possible variations between phyla. Each bar is clickable and produces a table with the list of sequences and associated cross-references (i.e. InterPro, UniProt). Sequences can be exported in multiple formats and/or stored in a user list.

Gene family structure: sequences are clustered at four levels of clustering, from less stringent to more stringent, in most cases narrowing the number of sequences (Figure 2B).

Protein domains: here, InterProscan was used to assess the domain conservation consistency and the specificity of the sequence clusters (Figure 2C). For each cluster, we performed statistical analyses to determine whether InterPro signatures were specific and therefore not found in any other sequences of clusters of the same level.

Phylogenomic analyses: this section includes multiple sequence alignments that can be downloaded and visualised using MSAviewer (33) and gene trees can be explored with InTreeGreat (https://www.southgreen.fr/content/intreegreat-tool) and PhyD3 (34) (Figure 2D). Some gene trees can be very large, and the interface proposes an option to prune automatically the tree based on user choice for a range of species.

Homologous predictions: the interface enables users to display and refine all the homologies detected by the phylogenetic-based approach and Reciprocal Best Hits (RBH). It is possible to filter and select only a subset of species of interest.

Figure 2.

Overview of gene family interfaces for GP001047: U2 auxiliary factor small subunit family. (A) diagram of the sequence count (327) by species (46) (B) Sequence flow in cluster structure (from level 1 to 4). (C) InterPro domain specificity statistics. About 92% of sequences have the U2 auxiliary factor small subunit which is uniquely found in this cluster. (D) Viewers for multiple sequence alignment (MSAViewer) and phylogenetic tree (InTreeGreat).

Pangene pages

The new page type for pangene sequences is a central and unique concept in this version as they were used for the clustering and homology predictions instead of all individual sequences that compose it. When browsing these pages, users can quickly see which genes are present or missing by looking at the status: core or dispensable compartments. Then, information related to the sequence composition is reported. Users can access information about pangene classification, the consensus sequence (except for singletons) with the multiple sequence alignments used to create it as well as related homology predictions (Figure 3). In the case of multi-copy clusters, it is possible to see which sequence was selected as representative (.rep) or participant (.p) and also why they were selected by browsing the distance matrix.

Figure 3.

Example of a pangene page (i.e. maize_pan_p014093) member of the U2 auxiliary factor small subunit gene family (GP001047). (A) Gene composition tab: all genes are part of core compartment for those four genomes since a representative sequence exists for each of the reference genomes. (B) Consensus sequences and associated multiple alignment. (C) List of homologs: the Zea mays pangene is predicted orthologous to pangenes in Triticum turgidum and in Oryza sativa. The Popups (green rectangles) display the sequence compositions of the respective pangenes. (.rep) refers to the representative sequence kept to create the consensus and (.p) to paralog sequences not used in the consensus.

New tools

The database can still be searched via keyword searches or by entering a query sequence for similarity search using DIAMOND (which replaces BLAST for faster processing), enhanced by new tools to further explore those datasets.

Manual curation and sharing of gene families

While automatic clustering is a relevant and efficient starting point, sometimes limitations (e.g. missing sequences, errors in gene annotations) are present and prevent access to ready-to-use datasets, justifying a deeper characterization that will eventually lead to a refinement of the automatic clustering. As a result, knowledge generated on individual gene families is often available only in publications and their supplementary information as PDFs. To encourage knowledge capture, we developed a section for advanced users to create and share their own gene families. Two methods are possible: users can either start from scratch and upload their data or use existing clusters and take advantage of multiples operations implemented in the ‘MyList’ features: such feature was indeed developed to facilitate intersecting, combining clusters. This new tool can be valuable during the review process by providing a unique identifier to referees—and eventually to users—to explore the structure and composition of the submitted gene family.