RGAugury: a pipeline for genome-wide prediction of resistance gene analogs (RGAs) in plants

Pipeline design

NBS and TM-CC containing proteins and membrane associated RLPs and RLKs belong to the main four known RGA families [5], and were thus included in the current version of RGAugury. RGA identification generally includes two basic steps: the identification of all conserved domains and/or motifs from protein sequences translated from gene sequences, followed by the analysis of the domain and/or motif composition including the classification of the genes into one of the four known types of RGAs based on such domain and/or motif structures. To improve the computing performance for the genome-wide identification of RGAs, three additional solutions were implemented in the pipeline: (1) an initial filtering of RGA candidates to remove a large portion of non-RGA genes, (2) parallel computing for time-consuming calculations, and (3) selection of the most efficient protein databases for domain detection. The design and the workflow of the pipeline are summarized in Fig. 1 and described in the following sections.

Initial filtering of RGA candidates

Because RGAs occupies a small percentage of the total genes in a genome, initial filtering to remove non-RGA genes can dramatically reduce the number of genes for the downstream domain/motif detection which saves considerable computing time. To do so, we used BLASTP to identify potential RGA candidates against an RGA database called RGAdb (see description below). The BLAST+ package was selected for sequence alignment as it outperforms BLAST in calculation power under the same conditions [36]. In the initial step, the input protein sequences are aligned against RGAdb by BLASTP using an E-value cut-off of 1e-5 (this E-value cut-off may be adjusted according to different species). Non-RGA proteins are filtered out (Fig. 1). Based on the analyses of 50 sequenced plant genomes, an average of 76.4 % of the annotated genes were removed, significantly abridging the downstream analysis time.

RGAdb was constructed using protein sequences and their annotations from the NCBI non-redundant protein database (nr) (http://www.ncbi.nlm.nih.gov/), plant resistance gene database (PRGdb, http://www.prgdb.org) [19]) and other cloned R-genes. Sequence entries were retrieved from nr with the key words ‘resistance’ and ‘disease’. Disease-irrelevant entries such as resistance to aluminum, drought, cold, DNA-damage, herbicide and UVB were removed. A total of 14,906 disease resistance related sequences were retained. In addition, ‘contributed’, ‘putative’ and ‘references’ entries from PRGdb were merged and the resulting dataset was filtered to remove redundant and non-annotated entries. As a result, 44,109 entries from PRGdb were appended into RGAdb. A few NBS encoding proteins from wheat and Brachypodium were also added into RGAdb [37–39]. Finally RGAdb contains 59,597 entries which were derived from more than 300 plant species (Additional file 1) and 80 % of the database entries are NBS coding proteins, RLKs, RLPs, MLO- and RPW-like proteins. In addition, this database will be regularly maintained to remove irrelevant entries and add new RGA entries.

Domain and motif detection

The initial filtering results in a reduced set of potential RGA candidates for the downstream domain and motif detection. There are seven RGAs-related domains and motifs, including NB-ARC or NBS, LRR, TM, STTK, LysM, CC and TIR in proteins. To detect them, four third-party tools were chosen and integrated into the pipeline. nCoils program was used to identify the CC domain present in CNL, CN and TM-CC types [35]. Minor modifications on the original source code of nCoils were made under its redistribution permission to facilitate the RGAugury pipeline to call nCoils through command-line. For detection of the NB-ARC domain, the pfam_scan toolkit was preferred because it outperformed InterProScan according to the test results obtained using the Arabidopsis genome data (data not shown). In addition, unlike InterProScan, the P-value parameter in pfam_scan is adjustable offering flexibility for the adjustment of the P-value cut-off across plant species. However, InterProScan was chosen to identify LRR and LysM, which are two components of RLPs and RLKs that play a role in pathogen signal recognition [16]. InterProScan is also suitable for the detection of the STTK domains of RLKs [25, 40–42]. Two tools, Phobius [32] and TMHMM [33], are available for TM domain identification but Phobius was elected because it performed better than TMHMM [32]. RGAugury’s configuration file can be easily modified to add any additional domains or motifs for detection, thus offering flexibility and extendibility.

InterProScan is a protein domain identifier that uses up to 14 databases for detection [43]. The accuracy and computing performance depend on the number and nature of the databases used for analysis. Pfam, Gene3D, SMART, Superfamily and the external database Panther were chosen for RGA domain detection based on accuracy and computing performance for the detection of LRR, LysM and STTK domains/motifs (data not shown). Three modes for database selection are provided in the pipeline. The default ‘Quick’ mode uses only Pfam and Gene3D which can identify most domains and motifs from input protein sequences. In the ‘Deep’ mode, Superfamily and SMART are also included for domain detection. As a consequence of the additional databases, the ‘Deep’ mode may perform slightly more slowly but this reduction in performance must be weighed against the accuracy and completeness of the results coveted. Differences between the two modes were mainly in the numbers of identified proteins with LRR motifs. The third mode is called ‘Free’ mode. This mode allows users to select one or more databases from a list. This mode is usually intended for result confirmation but can also be of use in custom applications. Here, the Panther (Protein Analysis Though Evolutionary Relationships) [44] database is introduced to the pipeline. Panther is a large protein database for comprehensive protein evolutionary and functional classification that contains more than 12,000 protein families from 104 sequenced genomes [45]. Our test results indicated that Panther would require huge calculation resource upon CPU threads for genome-scale protein analysis, which may take days to weeks depending on the input sequences. Thus, it is not suitable for genome-scale RGA identification but it can be accessed in the ‘Free’ mode for specific purposes such as the confirmation of small numbers of identified RGA candidates for example.

RGA identification

Once all domains and motifs are identified from the input protein sequences, RGAugury creates a table that lists the RGA candidates initially identified by BLASTP, along with their identified domains and coordinates on genes. The RGA identification module was developed to classify genes as potential RGAs and to classify them into specific RGA families (Fig. 2). First, the program classifies genes containing an NB-ARC or a TIR domain into the NBS encoding family. Genes that contain a TM domain may belong to the RLP, RLK or TM-CC families depending on the presence of other domains in the gene structures. If no NB-ARC, TIR or TM is observed, the gene is discarded as a non-RGA. The NBS-encoding gene family members are further divided into several subgroups according to their domain architecture, namely NBS, CNL, TNL, TN, CN, NL, TX and OTHER that may have chimeric domain/motif architecture. Here C, N, L, T and X represent CC, NBS or NB-ARC, LRR, TIR and unknown domains, respectively. For example, a gene with CC, NB-ARC and LRR domains is classified in the CNL subgroup. However, a gene with both TIR and CC domains ends up in the OTHER subgroup. In this case, a user needs to manually check this unexpected domain combination using the exported domain information. A gene with only an NB-ARC domain would require confirmation by InterProScan.

Data input and result output

The command-line pipeline requires three parameters. The first parameter is a protein sequence file in FASTA format from either a whole genome annotation project or manually annotated protein sequence data. The second one is the P-value cut-off for the initial RGA filtering using BLASTP. The default value is set to 1e-5. The last parameter informs on the database(s) to be queried for domain detection using InterProScan. RGAugury provides an optional parameter for specifying a companion Generic Feature Format (GFF3) or Gene Transfer Format (GTF) file when they are available. GFF3 or GTF files are helpful to draw RGA distribution plots using CViT [46] and draw gene and RGA domain structure plots. For most sequenced plant genomes, GFF3 annotation files can be downloaded from Phytozome [47].

Processing files derived from domain detection step and final RGA identification results are exported to plain text files. When a GFF file is provided, distribution plots of RGA families on chromosomes and the gene and domain structures for each RGA candidate are also generated.

Implementation of the command-line version

A single command-line pipeline program was written in Perl to seamlessly process all analyses. As domain detection for genome-scale protein sequences is computationally intensive, performance was a critical factor considered during the design of RGAugury. Benefitting from the development of programming techniques and hardware, parallel computing can significantly boost large-scale data analysis. To parallel heavy analyses, a Perl module, tool.pm, which invokes the fork function to parallel the execution of BLAST searches in the initial RGA filtering and domain detection, was implemented.

Web interface

To ease the command-line operations, a user-friendly web interface was designed and implemented using PHP and Java script (Fig. 3a). The necessary FASTA formatted protein sequences can be copied and pasted into the sequence text box or uploaded from a sequence file. A GFF3 file corresponding to the input protein sequences is optional but recommended. Databases for InterProScan can be selected through one of three modes: ‘Quick’, ‘Deep’ or ‘Free’ (flexible selection of up to five available databases). The default P-value cut-off for initial RGA filtering was set to 1e-5 but can be altered at user’s wish.

A project name is used as an intuitive ID for project management purposes. A ‘Job List Status’ page is dynamically generated upon project submission or following a mouse click on the ‘Status of Submitted Projects’ link (Additional file 2: Figure S1A). Finished jobs and current status of jobs being processed are displayed. Finished job are marked with the status “complete” in green that can be updated with the “Check” link to the ‘Results and Summaries’ page that shows the summary information, such as the numbers of RGAs and the detailed domain structure thumbnails for individual RGAs (Fig. 3b). The ‘Results and Summaries’ page also provides a link to download intermediate and final result files. All RGA candidates listed can be browsed. A search box is provided to query the candidates using full or partial letters or numbers from the table. For each RGA, an individual page will be generated by clicking on the Gene Structure Viewer (GSV) icon to show the details of the gene structure (Additional file 2: Figure S1B). GSV also provides links to the EnsemblePlants and GenBank websites for the genes that are available in these two databases, a useful feature for accessing additional gene information. When a GFF3 file is submitted along with the protein sequences, figures of RGA distribution on chromosomes are shown by clicking on the ‘Whole genome RGA distribution can be viewed’ link on Results and Summaries page (Additional file 2: Figure S1C).

The web version of RGAugury provides a new feature for InterProScan version check. When InterProScan runs on a web server, the web program will automatically check over the Internet whether the local InterProScan and its domain databases are up-to-date. If not, the web program pop-up dialog box reminding the user to update the software and its databases. In addition, detailed help pages are provided and synced with the RGAugury Bitbucket Wiki page (https://bitbucket.org/yaanlpc/rgaugury/wiki). Any updates on the Wiki page will also be reflected on the web help pages.

The RGAugury web program provides a project management function. Users can cancel executing jobs or delete finished projects from the server. Overall, the web version of RGAugury provides flexibility, convenience and interactive functionality over the command-line pipeline without compromising the RGA identification capabilities. The RGAugury web program can run on an HTTP server such as Apache. As hours of computing time are needed to finish a genome-scale RGA prediction, RGAugury is not suited for installation on a public web server. However, users are advised to download the pipeline and its web program from Bitbucket (https://bitbucket.org/yaanlpc/rgaugury) and install them in users’ local servers which can be used through intranet.

RGA identification accuracy

A similar evaluation was performed for RLP- and RLK-type RGAs. The previous report [50] indicated that the Arabidopsis genome encodes 57 RLPs, of which three have no sequence in TAIR 10, precluding their further analysis. Out of the remaining 54 RLPs, 46 (85.2 %) were predicted and 8 remained unclassified as a consequence of the absence of a TM domain (Table 1, Additional file 3). For RLKs, a total of 608 members were reported in the Arabidopsis genome, but only 75 % (456) of them were claimed to be associated with membranes [25]. RGAugury validated all 456 membrane-associated RLKs and identified additional four RLKs (Table 1, Additional file 4). The recently cloned wheat stem rust disease resistance genes sr33 and sr35 were also validated by RGAugury as NBS-LRR resistance genes [37, 38]. These case studies demonstrate that RGAugury accurately predicts RGAs on a genome-wide scale, and can be used for R-gene annotation.

Computing performance and multiple thread optimization

RGAugury can run in either ‘Quick’ or ‘Deep’ mode which rely on the number of databases used for domain detection. Differences between the two modes resided mostly in the numbers of identified LRR containing proteins (Table 2). For example, LOC_Os11g42590.1, an NL type in rice, was identified as NBS type using the ‘Quick’ mode but as an NL type using the ‘Deep’ mode. The LRR domain in this gene can be identified by either Superfamily or SMART database through InterProScan. However, these differences were minor in most cases (Table 2). Thus, the ‘Quick’ mode is recommended for fast RGA annotation.

Note: Two modes for database selection were used: the ‘Quick’ mode (Pfam + Gene3D) and the ‘Deep’ mode (Pfam + Gene3D + SMART + Superfamily). Results were separated by a slash if differences existed between the two modes. Plants were sorted by taxonomic groups which are labelled on the left side of the table. A. coerulea: Aquilegia coerulea; A. halleri: Anemone halleri; A. lyrata: Arabidopsis lyrata; A. thaliana: Arabidopsis thaliana; A. trichopoda: Amborella trichopoda; B. distachyon: Brachypodium distachyon; B. rapa: Brassica rapa; B. stricta: Boechera stricta; C. clementina: Citrus clementina; C. grandiflora: Capsella grandiflora; C. papaya: Carica papaya; C. reinhardtii: Chlamydomonas reinhardtii; C. rubella: Capsella rubella; C. sativus: Cucumis sativus; C. sinensis: Citrus sinensis; C. subellipsoidea: Coccomyxa subellipsoidea; E. grandis: Eucalyptus grandis; E. salsugineum: Eutrema salsugineum; F. vesca: Fragaria vesca; G. max: Glycine max; G. raimondii: Gossypium raimondii; K. laxiflora: Kalanchoe laxiflora ; L. usitatissimum: Linum usitatissimum; M. acuminata: Musa acuminata; M. domestica: Malus domestica; M. esculenta: Manihot esculenta; M. guttatus: Mimulus guttatus; M. pusilla: Micromonas pusilla; M. truncatula: Medicago truncatula; O. lucimarinus: Ostreococcus lucimarinus; O. sativa: Oryza sativa; P. hallii: Panicum hallii; P. patens: Physcomitrella patens; P. persica: Prunus persica; P. trichocarpa: Populus trichocarpa; P. virgatum: Panicum virgatum; P. vulgaris: Phaseolus vulgaris; R. communis: Ricinus communis; S. bicolor: Sorghum bicolor; S. italica: Setaria italica; S. lycopersicum: Solanum lycopersicum; S. moellendorffii: Selaginella moellendorffii; S. polyrhiza: Spirodela polyrhiza; S. purpurea: Salix purpurea; S. tuberosum: Solanum tuberosum; T. cacao: Theobroma cacao; V. carteri: Volvox carteri; V. vinifera: Vitis vinifera; Z. mays: Zea mays

Though BLASTP has its own multiple thread parameter (-n), the paralleled tool, tool.pm, developed herein, outperformed it. Two threading methods were compared using 40 CPUs. BLASTP searches were performed for 5000 and 10,000 protein sequences randomly selected from the Arabidopsis genome against RGAdb using two threading methods (Additional file 5). An improvement of 3.8- and 3.4-fold in execution speed for the two datasets, respectively, was observed for the tool.pm module when compared to the thread parameter (-n) of BLASTP. Speed improvement of the domain detection programs nCoil, Phobius and pfam_scan theoretically depends on the called CPU numbers because they lack embedded multiple threads function. Thus, application of tool.pm to these programs significantly reduces the time needed for domain detection, making it possible to shorten computing time for genome-scale RGA identification to mere hours.

To further evaluate the performance of RGAugury, we predicted RGAs for 50 sequenced plant genomes, including 44 angiosperm species and 6 green algae. All protein sequences were downloaded from the Phytozome (ver11.0) database website [47]. With both ‘Quick’ and ‘Deep’ modes, the processing time was highly correlated with the number of coding genes in genomes (R ²  = 0.76 for the ‘Quick’ mode and 0.75 for the ‘Deep’ mode) (Fig. 4). The total processing time averaged 2.4 h in ‘Quick’ mode. The ‘Deep’ mode averaged an extra 0.5 h but resulted in the enhanced prediction of LRR domains in some cases. Maize (Zea mays) data was not included in the correlation analysis because only 15.9 % of the input protein sequences were retained after the initial BLATP filtering compared with an average of 23.6 % among these 50 plant genomes, resulting in significantly fewer hours needed for the downstream domain identification. In general, RGAugury saves an average of 76 % time in the downstream domain/motif analysis if the initial BLASTP filtering step is performed before domain detection, demonstrating that the RGAugury pipeline is efficient for large-scale genome-wide RGA identification.

Large-scale genome-wide RGA identification for comparative analysis

A summary of the RGA prediction results for the 50 sequenced plant genomes is shown in Table 2. NBS encoding genes are divided into two major groups based on their variable N-terminals: TIR and non-TIR [51]. It has been hypothesized that TNL have never evolved in monocots [52] or have been discarded during evolution [51, 53–55]. Tarr et al. [56] developed an indirect method to search for TNL in monocots and magnoliids based on degenerate PCR and, they did not observed any TNL in the studied monocots. Among the 50 genomes studied here, nine belong to three different monocot orders: Poales (Brachypodium distachyon, Oryza sativa, Panicum hallii, Panicum virgatum, Setaria italic, Sorghum bicolor, Z. mays), Zingiberales (Musa acuminata) and Alismatales (Spirodela polyrhiza). Our RGAugury prediction of these monocot species did not identify any TNL type confirming the previous findings (Table 2). Genes encoding TN or TX proteins were detected in these monocots but the numbers were smaller compared to genes predicted to encode other types of NBS motifs and compared to those from species other than monocots, implying that TNL may have been discarded during the evolution in the monocot lineage.

No genes coding for NB-ARC domains were predicted from the model green algae Chlamydomonas reinhardtii, Volvox carteri, and Cyanidioschyzon merolae because TIR and non-TIR NBS encoding genes arose only in the Plantae [57]. Our prediction results for all six green algae organisms confirmed this observation (Table 2), providing additional validation for the RGAugury prediction ability.