ReGSP: a visualized application for homology-based gene searching and plotting using multiple reference sequences

Workflow

The overall workflow of ReGSP is shown in Fig. 1. The analysis begins by setting a query and reference sequences by the user. The program then evaluates every query contig against each reference sequence, based on two filtering parameters(see section Gene search parameters), following a BLAST analysis. BLAST is preferentially used for finding HSPs because of its speed and proven ability to find similarity regions (Mount, 2007). ReGSP uses BLASTx tool from the BLAST family, as it performs well for identifying similarities between the sequences of known genes in a reference genome and a provided nucleotide query sequence. Genes that satisfy the two filtering parameters are listed in the final output files, i.e., a .regsp text file and a MUMmer plot PDF file. In parallel, ReGSP performs nucleotide k-mer analysis using the cPlot visualization tool. A cPlot PDF file presents a sequence similarity dot plot of all queries plotted against the reference sequences. In addition, ReGSP provides a .sorted text file similar to the .regsp text file, which contains all the results before the second filtering step.

Complete and accurate reference genomes are key to obtaining information on novel genes in a homology search. In some cases, a single reference genome does not provide all the detailed and diverse information required for gene search (Ballouz, Dobin & Gillis, 2019; Yang et al., 2019). Accordingly, ReGSP accepts multiple reference sequences that are phylogenetically close to the query, to validate contigs and identify genes from BLAST fragments obtained after an alignment of each reference and the query. Query sequences can be contigs or scaffolds from newly sequenced or previously published genomes, whenever gene analysis is required.

Gene search parameters

ReGSP uses two parameters, the maximum number of matched BLAST hits per group (N) and the minimum number of sub-gene counts per contig (M), for gene searching. The maximum number of matched BLAST hits per group, N, controls the number of allowed overlapping BLAST hits per group, where a group is a non-overlapping section of the contig with a BLAST match against the reference. If a contig search yields overlapping BLAST hits, all matches are considered as one group. Hence, a group can have a maximum of N − 1 hits. If the number of hits exceeds this number, the hits will be sorted based on the BLAST score, and the top scoring N − 1 hits will be considered to be candidate genes. The second parameter for candidate gene identification is the minimum number of BLAST gene matches per contig, M. This parameter dictates the validity of a contig based on the total count of genes in a contig following BLAST analysis and satisfying the first filtering parameter N. Note that the genes (matched BLAST hits) in a group are counted individually. If a gene count for a contig does not match the value indicated by M, it will not be considered as a candidate gene. To put the two parameters in perspective, let us consider a scenario with two contigs as a query sequence and three reference genomes for which the coding sequences are known. Let us set the first parameter, i.e., the maximum number of matched BLAST hits per group (N), as 4 and the minimum number of sub-gene counts (M) also as 4. Figure 2 shows a visual representation of the two contigs, in which ReGSP attempts to identify candidate genes.

A detailed view of Reference 1 (top section of Fig. 2) reveals five BLAST hits with Contig 1. The first two matched genes (BLAST match a and BLAST match b) share an overlapping region; hence, they are considered a group. However, since the maximum number of allowed genes in a group is four, up to three overlapping hits are allowed (i.e., N − 1); therefore, both BLAST matches are considered hits. The two overlapping hits and the remaining three hits (c, d, and e) are five sub-gene hits for Contig 1. This qualifies Contig 1 to satisfy the second parameter (minimum sub-gene count), which was set to 4. Similarly, 8 BLAST hits against Reference 1 were detected for Contig 2. Here, BLAST hits 5, 6, 7, and 8 share an overlapping region, and belong to the same group. However, the maximum number of allowed genes per group exceeds the first parameter N (a total of 3 allowed, as N − 1). Hence, the hits are ordered in a descending manner based on the E-value, and the hit (fragment) with the lowest similarity (in this case, BLAST hit 8) is discarded, with seven BLAST hits remaining. Those are hits 1, 2, 3, and 4, and the overlapping areas 5, 6, and 7. Contig 2, with seven sub-gene counts, has a higher sub-gene count than that set in the second parameter, which makes it valid.

Consequently, no BLAST hits are identified for Contig 1 against Reference 2, whereas three are identified for Contig 2. For Contig 2, the number of BLAST hits (3) is below that set by M; hence, none of the two contigs contain a valid number of genes according to ReGSP. For Reference 3, however, only Contig 1 has a valid set of genes, with four fragments, whereas the single BLAST hit for Contig 2 fails to satisfy the M criterion.

Once the genes have been identified for each contig using the two search parameters, ReGSP records the analysis result in a .regsp text file for the user to download. Choosing the appropriate values of the gene search parameters is an adaptive process, and depends on the length or gene content of the references, and the evolutionary distance between each reference and the query.

Plotting tool

ReGSP uses PROmer to identify regions of conserved protein sequences in a given contig and MUMmer plot to display the results. MUMmer enables rapid alignment of DNA and amino acid sequences, and its graphical module allows assembly viewing, and comparison of alignments of sequences from closely or distantly related species. Hence, in ReGSP, following contig gene validation using the two filtering parameters, a MUMmer plot will reveal conserved regions of contigs against each reference sequence.

Further, ReGSP incorporates a new dot plot tool called cPlot (Ji, Jung & Yi, 2021). This novel tool identifies nucleotide sequence similarity between the provided query sequences and the selected or provided references, to generate a dot plot. Initially, cPlot extracts common reads from the query and reference sequences based on k-mer size, and then merges the query sequences that are consecutive when aligned. The aligned reads are k-mer–sized and indexed for use at subsequent analysis steps. Next, cPlot identifies an optimal read alignment among the overlapping read candidates (see Fig. 1). To do this, reads that match a provided reference are selected. During this process, the longest base alignment read that is most similar to the reference is selected based on the sorted read index generated at the initial analysis stage. cPlot processes all the reads, starting with the read closest to the region at the beginning of the reference sequence, and collates overlapping reads into a subset. After optimal alignment, the query sequence is rearranged and sorted according to the reference sequence.

cPlot operates in the nucleotide domain, provides additional level of sequence similarity information for ReGSP, and helps to visualize short-read assembly. The original cPlot script was developed to handle a single reference sequence for contig plotting. However, we have modified it to handle multiple references to suit the needs of ReGSP. The final dot plot result of cPlot is provided as a PDF file for downloading.

Web interface

ReGSP can be accessed at https://ds.mju.ac.kr/regsp. The web interface is freely accessible, and allows anyone who is interested to perform reference-based gene search in a simple and convenient way. The first requirement for using the interface is to provide a reference sequence, and the user has multiple options for reference inputting. The accession number of the reference sequence can be simply typed in a text box. Since ReGSP can handle multiple reference sequences for gene searching, multiple references can be entered separated by a comma. The other option for selecting a reference is using the search tree button, which brings up a list of all organelle genomes from the NCBI Organelle Genome Resources. The organelle search tree contains more than 19,000 organelle genomes, and the user can search and select an organism based on its name. Furthermore, related species from the same genus, family, and kingdom can be selected. On the other hand, if the users want to use their own custom reference sequences, they are required to provide a nucleotide reference sequence or multi-FASTA nucleotide sequences, and an amino acid sequence of the references. The nucleotide and amino acid sequences are uploaded separately. The two input parameters of ReGSP, N and M, must be set using natural numbers in the two side-by-side text boxes on the interface. The web interface additionally offers an option to select an appropriate k-mer size for use in cPlot visualization. A default k-mer value of 13 will be used unless the user provides a specific number between 10 and 30. Figure 3 shows the ReGSP interface. Current tasks in progress, completed analysis results, and result history can be accessed on the side menu of the application page. The ReGSP web interface was built using an angular web framework (TypeScript programming language), with a Flask API to communicate with its Python back-end.

Data set

To evaluate the ReGSP performance, chloroplast genomes of six algal species were analyzed, with four species from the Chromista group and two species from the Rhodophytes. The genomes were retrieved from the public NCBI Nucleotide database. The complete genome of Heterosigma akashiwo from the Raphidophyceae clade was used as a query sequence to showcase the analysis. Of note, ReGSP is suitable for analyzing contigs and scaffolds against multiple references. For reference purposes, sequences of two species from the Bacillariophyceae (diatom) class were used: the pennate diatom Haslea silbo and Seminavis robusta. Another reference was Nannochloropsis limnetica from Eustigmatophyceae class. The remaining reference sequences represent two classes of red algae (Rhodophyta) which are Bangiopsis subsimplex and Cyanidium caldarium.

The dataset can be found in the ReGSP sample page (https://ds.mju.ac.kr/regsp/#/sample), and the details of all the sequences are shown in Table 1. The amino acid sequences for each reference that are used for BLASTx alignments and all the nucleotide sequences including the query were downloaded from NCBI Nucleotide database.

To identify the evolutionary relationship between the query and the reference sequences, 12 common genes harbored by the analyzed algae were selected, namely: psaA, psaB, psbA, psbB, psbC, atpA, atpB, atpI, petA, petB, rbcL, and rpl2. According to UniProtKB (Consortium, 2021), the psa and psb gene families encode protein subunits of photosystem I (PSI) and photosystem II (PSII), respectively, which drive oxygenic photosynthesis, whereas petA and petB encode components of the cytochrome b6-f complex, which mediates electron transfer between PSI and PSII. The genes atpA, atpB, and atpI are involved in the production of ATP, while the rbcL gene catalyzes carbon dioxide fixation and rpl2 is ribosomal protein large subunit involved in translation. For each species, the sequences of the 12 genes were combined and subsequently aligned using the MUSCLE program in MEGA X (Kumar et al., 2018). Figure 4 shows the phylogenetic tree constructed based on these genes, with their evolutionary history inferred using maximum-likelihood method and Tamura–Nei model (Tamura & Nei, 1993).

Text file output

ReGSP presents homology-based gene search results for the submitted query against the provided references in two text files (.regsp and .sorted). In the current study, all counts of gene matches for the Heterosigma akashiwo query against the references extracted from the .sorted text file are shown in Table 2. Customized and sorted BLASTx output is shown in the .sorted text file, before proceeding with the second parameter analysis, i.e., the minimum gene count per contig (M). The .regsp text file is similar to the .sorted text file, except that it includes the output after M parameter analysis. As shown in Table 2, the minimum number of genes per contig was greater than 50 for all references. Hence, in this case, the .regsp and .sorted text files showed identical results. However, for cases where the number of hits for contigs and/or scaffolds is smaller than M, the two text files show different gene counts. Further, reference sequences that are phylogenetically close to the query sequence have a higher number of BLAST hits than those for distant relatives with similar base length.

In addition to the number of gene hits, the two text files contain detailed gene information, such as sequence hit region, alignment length, bit score, translation frame, and the resulting protein alignment. These details provide valuable information for the analysis of the query sequence against the reference(s), and are sorted based on similarity to each reference. Figure 5 shows a small section of the .regsp text output. The various pieces of information were retrieved for all references, and gene search results were sorted based on similarity score.

Plot output

ReGSP generates two plot files. The first plot file is generated using the PROmer alignment program of MUMmer. For the analysis, all query contigs are aligned against each provided reference sequence in all (six) amino acid frames to reveal conserved regions in the sequences. The high sensitivity of PROmer leads to an improved alignment of sequences that are positioned close in the phylogenetic tree. Using the Heterosigma akashiwo query example, the analysis detected multiple conserved genes in all reference sequences, with a slight preference for species more closely related to the query, as shown in Fig. 6 and the phylogenetic tree in Fig. 4. Figure 6 shows the alignment of Heterosigma akashiwo query genome sequence, with every reference in a separate grid.

The second plot output of ReGSP is a nucleotide domain dot plot generated by cPlot. Figure 7 presents different k-mer–size plots of the Heterosigma akashiwo chloroplast genome with the provided reference sequences: a k-mer value of 10 (Fig. 7A) and a k-mer value of 13, the default k-mer value for ReGSP (Fig. 7B). Typically, multiple instances of sequence similarity are identified at small k-mer values for species that are evolutionarily close to the query, as well as those that are distant from the query. Small k-mer values are helpful when most of the references are distantly related to the query. However, as the k-mer value increases, only highly similar sequence alignments are detected, and the dot plot fades out for less similar sequences. Figures 7C and 7D illustrate this property, at k-mer values of 20 and 30, respectively. Therefore, high k-mer values reveal long segments of sequence similarity, with dot plots available for sequences that are closely related. On the other hand, low k-mer values, which require a relatively longer ReGSP execution time, will result in multiple short matched segments. For the Heterosigma akashiwo example, the query shares high sequence similarity with Haslea silbo and S. robusta sequences, which belong to the same clade it occupies. Further, the Heterosigma akashiwo chloroplast sequence shared the lowest similarity with algae from the Plantae kingdom. The ReGSP analysis findings concur with the phylogenetic relationship of the species, as shown in Fig. 4.

Finally, based on testing of different lengths of query and reference sequences, for most optimized visual plot results, we recommend providing no more than 30 multi-FASTA files as aquery and less than 10 reference sequences.