cRegions—a tool for detecting conserved cis-elements in multiple sequence alignment of diverged coding sequences

Implementation

The overall principle of the cRegions tool is to compare observed nucleotide frequencies to expected probability distribution and calculate appropriate metrics (described below) to detect regions where the coding sequence is more conserved than expected. Scripts used in the cRegions web tool are available in GitHub repository at https://github.com/bioinfo-ut/cRegions. The workflow of cRegions is as follows:

1. Two inputs are required: a protein multiple sequence alignment (MSA) and nucleic acid sequences containing coding sequences (CDS) of respective proteins. Both inputs have to be in FASTA format. mRNA or the full genome can be used instead of the exact CDS. However, the coding sequence must not contain introns.

2. Protein alignment is converted into a corresponding codon alignment using respective coding sequences with PAL2NAL (Suyama, Torrents & Bork, 2006).

3. Henikoff position-based sequence weights are calculated using the codon alignment (Henikoff & Henikoff, 1994).

4. Codon usage bias is calculated from the codon alignment. Calculated proportions are adjusted by Henikoff position-based sequence weights to account for non-uniform phylogenetic coverage (Fig. S1). Codon preference for serine in the TCN block and in the AG[A/G] are calculated separately.

5. Expected nucleotide proportions are calculated for each position in the codon alignment based on the amino acid sequence and the codon usage bias (calculated or provided by the user). Henikoff position-based sequence weights are used to adjust expected proportions (Fig. S2).

6. The observed nucleotide frequencies are compared to the expected probability distribution of nucleotides in each position. The comparison is made only for positions with amino acids having more than one codon. Three different metrics are used for this:

   1. The algorithm uses R (R Development Core Team, 2014) to calculate p-values of Chi-square goodness of fit test (chisq.test) for each column in the codon alignment. The test allows us to see whether the observed distribution of nucleotides is significantly different from expected distribution. We use the negative logarithm of the p-value of the Chi-square goodness of fit test as the metric. Bonferroni correction is used to show the threshold with significance level α = 0.05. $$p\text{-value}=\text{chisq}.\text{test}(\text{c}({\text{A}}_{\text{obs}},{\text{C}}_{\text{obs}},{\text{G}}_{\text{obs}},{\text{T}}_{\text{obs}}),\text{\hspace{1em}p}=\text{c}({\text{A}}_{\text{exp}},{\text{C}}_{\text{exp}},{\text{G}}_{\text{exp}},{\text{T}}_{\text{exp}}))$$ The subscript ‘obs’ indicates observed frequencies, the subscript ‘exp’ indicates expected proportions.

   2. The second metric is the root-mean-square deviation (RMSD). Only nucleotides which have a predicted probability over zero are included in the RMSD calculation. $$\text{RMSD}=\sqrt{\frac{1}{4}\left[{\left(A_{\text{obs}}-A_{\text{exp}}\right)}^2+{\left(C_{\text{obs}}-C_{\text{exp}}\right)}^2+{\left(G_{\text{obs}}-G_{\text{exp}}\right)}^2+{\left(T_{\text{obs}}-T_{\text{exp}}\right)}^2\right]}$$ The subscript ‘exp’ indicates expected frequencies.

   3. The third metric is the maximum difference (MAXDIF), which selects only a single nucleotide from each column having the highest absolute difference between predicted and observed values. $$\text{MAXDIF}=\max \left(|A_{\text{obs}}-A_{\text{exp}}|,|C_{\text{obs}}-C_{\text{exp}}|,|G_{\text{obs}}-G_{\text{exp}}|,|T_{\text{obs}}-T_{\text{exp}}|\right)$$

      The subscript ‘exp’ indicates expected frequencies.

In all cases, the larger numerical value of a metric indicates higher conservation at the nucleic acid level. Additionally, if a position in the codon alignment has more than 20% of gaps, the metric is not calculated for that position (Fig. S3).

Alphavirus dataset

In the present work, we used the non-structural and structural polyprotein of alphaviruses as an example. Sequences were downloaded from the NCBI viral genome database (non-redundant dataset, 16 April 2018). The first dataset consists of 24 known alphaviruses (Table S1). The non-structural polyprotein dataset was further divided into two sub-datasets. The first subset (SFV dataset) consists of seven viruses from genus Alphavirus, all belonging to the ‘SFV Complex’ (Table S1) (Forrester et al., 2012). The second subset was formed by nine ‘New World’ alphaviruses (Table S1).

cRegions and synplot2

Protein sequences were aligned with MAFFT (Katoh et al., 2002) using the default settings at http://www.ebi.ac.uk/Tools/msa/mafft/. Graphs were created with default settings using a sliding window with size one unless stated otherwise. Codon alignments for the synplot2 (Firth, 2014) were created with PAL2NAL (Suyama, Torrents & Bork, 2006), thus the input alignments for synplot2 and cRegions are identical. In this study the smallest possible window size was used for synplot2, giving the resolution of three codons (2n + 1 codons). In case of synplot2, significant hits were selected at threshold p < 10⁻⁵. It is less conservative compared to the threshold used in the synplot2 paper (Firth, 2014).

Sequence weighting

Henikoff position-based sequence weights are used to compensate for the over-representation of well-sequenced taxa in the MSA (Henikoff & Henikoff, 1994). Contrary to the original work of Henikoffs, we applied position-based weights to nucleotide sequences in the codon alignment, not to protein sequences. Thus, including variance at the codon level. Predicted nucleotide proportions for each position in the codon alignment are adjusted with sequence weights.

Sliding window mode

Cis-acting elements may be longer than a single codon, for example, dual-coding regions, thus the possibility to calculate a single metric over consecutive codons may be preferred. The cRegions web tool allows the user to set the window size from one to 1/9 of the length of the codon alignment. By default, the third codon position is used in the sliding window mode, as it is most informative. An additional threshold exists in the sliding window mode. The threshold is for skipping columns instead of terminating the metric calculation for these consecutive positions. For example, if there is an insertion in a single sequence, the position should be skipped and the next codon included in the current window instead. It should be noted that skipping can happen several times in a row. By default, if a position in an MSA has more than 90% of gaps, it is skipped in the sliding window mode. It should be noted that the threshold for skipping gaps and the threshold for metric calculation are different parameters (Fig. S3). In case of RMSD and MAXDIF, an arithmetic mean is calculated over consecutive codons. However, a p-value of Chi-square test is calculated based on observed values that are added over all consecutive codons. Again, Bonferroni correction is used to show the threshold with significance level α = 0.05. It should be kept in mind that adjacent positions with low metric values will decrease the value of a single conserved position if the window size is larger than one.

Visualisation

The cRegions web tool uses ‘highcharts’ libraries to visualise results (http://www.highcharts.com/). Alignment visualisation is provided by MSAViewer (Yachdav et al., 2016). The combination of highcharts and MSAViewer allows the user to pinpoint (by clicking on the point) and navigate directly to a conserved region or nucleotide. In addition to scatter plots of different metrics, cRegions web tool displays an interactive graph of codon usage. Codon usage is calculated over all analysed sequences. A table of codon frequencies and a file with tab-separated values are included in the downloadable zip container. The same codon table can be used as an input for the cRegions algorithm.

VEEV sequences

The Venezuelan equine encephalitis virus (VEEV) neighbour sequences were downloaded from the NCBI viral genomes database (https://www.ncbi.nlm.nih.gov/genomes/GenomesGroup.cgi?taxid=11018, 18 May 2018). A total of 11 entries were removed as they did not have annotated full-length non-structural polyprotein. Identical protein sequences were removed using jalview ‘remove redundancy 100’ (Waterhouse et al., 2009). The final dataset contained a total of 94 isolates, including 93 VEEV neighbour sequences and a reference VEEV sequence (NC_001449).

Randomly mutated protein-coding sequences

A random 3,000 nt long protein-coding sequence was created with SMS v2 tool (http://www.bioinformatics.org/sms2/random_coding_dna.html). A different number of random mutations (25–1,600) were introduced into that sequence with the SMS2 mutate tool (http://www.bioinformatics.org/sms2/mutate_dna.html) (Stothard, 2000). In total, we generated 3 × 8 datasets. Each set consisted of one original randomly generated protein-coding sequence and six, nine or 14 randomly mutated sequences with a different number of mutations per bp. For the set with 15 sequences we generated different initial protein-coding sequence to remove a bias which could occur if we only use one protein-coding sequence as a seed. These sequences do not need aligning, because homologous nucleotides are already aligned.

Threshold correction for nearly identical sequences

cRegion algorithm assumes that, in general, the distribution of nucleotides at each position in the MSA correspond to an average codon usage of the same sequences under analysis. The assumption is reasonable if the sequences under analysis have diverged. However, in the case of nearly identical sequences, the distribution of nucleotides in each position in the MSA is more similar to observed nucleotide proportions rather than to average codon usage. Therefore, the expected proportions are more similar to observed proportions. Thus, even if using Bonferroni correction, the Chi-square test may still give many potentially false positive signals. Therefore, we need to adjust the threshold such a way that in the case of nearly identical sequences the threshold is stricter. For that, we also include the average pairwise identity of nucleotide sequences into calculations of expected nucleotide proportions at each position. The expected nucleotide proportions are adjusted by the observed proportions which depend on the average pairwise identity of nucleotide sequences. The exponent value was found empirically. $${\left(A_{\text{exp}},C_{\text{exp}},G_{\text{exp}},T_{\text{exp}}\right)}_{\text{adj}}=\left(A_{\text{exp}},C_{\text{exp}},G_{\text{exp}},T_{\text{exp}}\right)+i_n{}^7*\left(A_{\text{obs}},C_{\text{obs}},G_{\text{obs}},T_{\text{obs}}\right)$$ i_n = the average pairwise identity of nucleotide sequences

In addition to the previous adjustment of expected values, also the threshold itself is adjusted. The threshold correction depends on the ratio between the average pairwise identity of nucleotide and protein sequences and the number of sequences in the MSA.

• t = threshold − log10(p–value)

• i_n = the average pairwise identity of nucleotide sequences

• i_p = the average pairwise identity of protein sequences

• d = the number of sequences in the multiple sequence alignment

Alphaviruses

First, we applied cRegions and synplot2 on all 24 non-structural polyproteins of alphaviruses (see also ‘alphavirus dataset’ example on the cRegions homepage). We detected a total of six significant signals with cRegions (Fig. 1A) and three significant signals with the synplot2 (Figs. S4A and S5A). The first signal from the 5′ end was recognised by both programs (Fig. 1A) and spanned from positions 138 to 174 in the codon alignment (Table 1). It is a conserved sequence element (CSE) called ‘51 nt CSE’, which acts as an enhancer for the RNA synthesis, affecting viral replication. This CSE forms two stem-loops and is located at positions 155–205 in the Sindbis virus (SINV) genome (Niesters & Strauss, 1990). Thus, the detected signal lies exactly in the region (Table 1).

The second significant hit, a single nucleotide at position 1,149, was detected only by the cRegions algorithm. However, two adjacent positions 1,143 and 1,146 were just below the threshold. In the New World alphavirus dataset, in addition to position 1,149, 1,143 was also significant. The signal is just adjacent to the packaging signal of SINV and New World alphaviruses (see also ‘New World alphavirus dataset’ example on the cRegions the homepage). It has been shown that a 570 nt fragment positions 684–1,253 from the SINV binds to the viral capsid protein and is required for packaging of SINV. The detected signal lies in this region (Weiss et al., 1989). However, when we analysed VEEVs separately we were able to detect the positions of phylogenetically conserved predicted stem-loops (Fig. S7). The results are similar to the work done by Kim et al. (2011).

The third and the fourth signal are also single nucleotides at positions 2,835 and 2,967, respectively (Fig. 1A). Both signals are located inside nsp2 conserved region called region b. This 266-nucleotide region is located from 2,726 to 2,991 in the SFV genome (White, Thomson & Dimmock, 1998). Previous deletion mutation analysis has shown that nucleotides from 2,767 to 2,824 in the b region are required for efficient packaging of SFV genome. (White, Thomson & Dimmock, 1998). The first signal is located in that region. Additionally, analysis of ‘SFV Complex’ viruses separately led to increased significance of the first signal (Fig. 1C) and the same signal became visible with synplot2 (Figs. S4B and S5B). Expectedly, both signals disappeared in the New World dataset, as the packaging signal is in a different location in these viruses (Fig. 1B). Therefore, dividing datasets to different subsets may help to detect signals that are only characteristic to smaller subgroups.

The fifth significant hit is a single nucleotide at position 6,834. It is downstream of the ‘leaky’ stop codon (stop codon is at 6,814–6,816 on the codon alignment and in the SINV genome at nt 5,748–5,750). Synplot2 was able to detect a much larger region compared to cRegions in the same area (Figs. S4A and S5A). The detected signal is a 3′ stem-loop RNA secondary structure immediately adjacent to the stop codon (+13 nt downstream of the stop codon in SINV). For many alphaviruses, including VEEV and SINV, it has been reported to influence read-through. In the SINV genome, the double helix part (the stem) of the stem-loop is predicted to form between the two regions: 5,763–5,772 and 5,928–5,939 (Firth et al., 2011). Therefore, the detected signal at position 6,834 (5,768 in SINV) is inside the first region. However, when we analysed VEEVs separately, we were able to detect multiple significant signals inside this stem-loop region (Fig. S7).

The sixth signal consists of two positions 8,658 and 8,664 on the codon alignment (Fig. 1A; Figs. S4A and S5A). The signal is located within the subgenomic promoter of alphaviruses (Raju & Huang, 1991; Rupp et al., 2015).

The cRegions and the synplot2 were also applied to the structural polyproteins of alphaviruses (see also ‘alphavirus structural dataset’ example on the cRegions homepage). Sliding window size 2 was used with cRegions. A strong signal was detected in positions 2,643–2,649 on the codon alignment, which corresponds to a UUUUUUA motif (Fig. 2). The motif is responsible for a frameshift in a structural protein (Firth et al., 2008; Chung, Firth & Atkins, 2010). The same signal was detected with the synplot2 (Fig. S6).

Requirements on sequences

The method used in cRegions has some limitations and prerequisites (Puustusmaa & Abroi, 2016). First, the sequences under study must have diverged. Second, the embedded functional element must have been under selection. To help users to evaluate their sequences in these aspects, we added an interactive version of Fig. 3 to the web tool. The plot visualises the sequences under study in comparison to randomly mutated sequences and sequences thoroughly analysed in the previous or current study with respect to divergence and selection. To evaluate divergence and selection we used the relationship between average pairwise nucleotide identity and average pairwise amino acid identity (Fig. 3). As shown in Fig. 3, randomly mutated simulated sequences form a clear and narrow assembly on the plot. Randomly mutated sequences with a defined extent (N mutations per bp) were used to model neutral evolution and/or non-diverged sequences (more details in ‘Materials and Methods’). The naturally occurring sequences used in the previous and in the current study locate clearly away of the simulated sequences.

Sequences having low divergence or/and having close to neutral evolution

As cRegions was designed to work on diverged sequences, the method may give potential false positive signals in low divergence sequences or in sequences locating close to neutrally evolving sequences (Fig. 4; Fig. S8). To avoid this, we recommend enabling threshold correction on the web tool. By enabling this option expected values are corrected with observed values and the adjusted threshold is calculated (see ‘Materials and Methods’). This removes most of the false positive signal from sequences that are close to randomly mutating sequences (Fig. 4; Fig. S8). We would like to note that the correction is needed only in the case of sequences which are close to the neutrally/randomly evolving sequences (Fig. 3). Another option is to use synplot2 which uses neutral evolution as its null hypothesis (Firth, 2014)