Advanced Options and Details

PDB 3D structure

Three-dimensional structures of biological macromolecules are available in the Protein Data Bank (PDB)(Rose et al., 2015), a copy of which is accessible from the ConSurf server. AlphaFold model structures of most other proteins are available in UniProt. The HMMER, CSI-BLAST or PSI-BLAST (Finn et al., 2011; Biegert and Söding, 2009; Altschul et al., 1997, 1990) search for homologues is done using the target chain sequence extracted from the SEQRES record of the PDB file as the input query. If the SEQRES record does not exist, ConSurf extracts the sequence from the ATOM record. Thus, a user-provided PDB file must include the ATOM record.

Searching for homologous sequences

For proteins, the server uses the MMseqs2, HMMER, CSI-BLAST or PSI-BLAST (Finn et al., 2011; Altschul et al., 1997; Biegert and Söding, 2009; Steinegger and Söding, 2017) heuristic algorithm with default parameters to collect homologous sequences of a single polypeptide chain of known 3D-structure. The search can be carried out using the following databases:

1. UniRef90 - consists of cluster sequences and sub-fragments with 11 or more residues that have at least 90% sequence identity with each other (from any organism) into a single UniRef entry, displaying the sequence of a representative.
2. SWISS-PROT - a curated protein sequence database which strives to provide a high level of annotation.
3. UniProt - the universal protein resource, a central repository of protein data created by combining SWISS-PROT, TrEMBL and PIR.
4. Clean UniProt - a modified version of the UniProt database aimed to screen the more reliable sequences using two criteria:
   
   • a. If the "Decription" (DE) field contain "Disease", "RIKEN", "variant", "mutation", "mutant" or "whole genome shotgun sequence" the sequence is removed;
   • b. If the database is "TrEMBL" and the "Comments" (CC) lines contain the word "CAUTION" the sequence is removed.
5. NR database - 'non-redundant' database (i.e. with duplicated sequences removed). NR contains non-redundant sequences from GenBank together with sequences from other databanks (Refseq, PDB, SWISS-PROT, PIR and PRF).

The user can also set the maximum number of homologs to collect and the number of iterations while the default runs one HMMER iteration, with maximum of 150 homologs and E-value cutoff of 0.0001 (E-value describes the number of hits one can "expect" to see just by chance when searching a database of a particular size. The higher the E-value, the more hits will be expected, but the pairwise distance between them and the query sequence will increase). The user can also control the maximal percentage identity between sequences, removing redundant sequences and specify the level of redundant sequences for removal. The sequences found are clustered by their level of identity using CD-HIT (Li and Godzik, 2006) and the cutoff specified by the user (default level is 95% identity). The minimal percentage for homologs is set by default to 35% which is the level of the upper bound of the 'twilight zone' for protein structures (Rost, 1999). See figure 2 for user's MSA parameters screen on the web server.

Figure 2: Screenshot of the parameters window to build the Multiple Sequence Alignment (MSA) on the Consurf web server.

Note: The best performance is obtained when single domains are used. Using a polypeptide chain that contains several domains often results in finding proteins that are homologous to one of the domains but not to the others. Therefore it is advisable to analyze multi-domain proteins one domain at a time.

For nucleic acid sequence the search is performed using the HMMER/BLAST algorithm (Biegert and Söding, 2009; Altschul et al., 1990) on the NCBI NT database. The default runs one BLAST iteration with E-value cutoff of 0.001.

Generating the phylogenetic tree

The server constructs a phylogenetic tree that is consistent with the available MSA, using the neighbor joining (NJ) algorithm (Saitou and Nei, 1987), as implemented in the Rate4Site program (Pupko et al., 2002). The server can also accept a user-provided phylogenetic tree (in Newick format). This option can be used only if a corresponding MSA is provided.

Calculating the conservation scores

The conservation score at a site corresponds to the site's evolutionary rate. The rate of evolution is not constant among amino (nucleic) acid sites: some positions evolve slowly and are commonly referred to as "conserved", while others evolve rapidly and are referred to as "variable". The rate variations correspond to different levels of purifying selection acting on these sites. For example, in proteins, the purifying selection can be the result of geometrical constraints on the folding of the protein into its 3D structure, constraints at amino acid sites involved in enzymatic activity or in ligand binding or, alternatively, at amino acid sites that take part in protein-protein interactions.

In ConSurf, the rate of evolution at each site is calculated using either the empirical Bayesian (Mayrose et al., 2004) or the Maximum Likelihood (Pupko et al., 2002) paradigm. In both of these methods, the stochastic process underlying the sequence evolution and the phylogenetic tree are explicitly taken into account. The Bayesian method was shown to significantly improve the accuracy of conservation scores estimations over the Maximum Likelihood method, in particular when a small number of sequences are used for the calculations (Mayrose et al., 2004). An additional advantage of the Bayesian method is that a confidence interval is assigned to each of the inferred evolutionary conservation score.

Analyzing nucleic acid sequences

Despite increasing interest in the non-coding fraction of transcriptomes, the number, the level of conservation, and functions, if any, of many non-protein-coding transcripts remain to be discovered. However, it has already been shown that many of the non-coding sequences are connected to regulatory processes. ConSurf offers estimations of the evolutionary rate for each position of nucleic acid sequences in the same manner used for amino acid residues. For that purpose, four evolutionary models were implemented in the Rate4Site program: (i) the Juke and Cantor 69 model (JC69), which assumes equal base frequencies and equal substitution rates (Jukes and Cantor, 1969). (ii) The Tamura 92 model that uses only one parameter, which captures variation in G-C content (Tamura, 1992). (iii) The HKY85 model, which distinguishes between transitions and transversions and allows unequal base frequencies (Hasegawa et al., 1985). (iv) The General Time Reversible (GTR) model, which is the most general time-reversible model. The GTR parameters consist of an equilibrium base frequency vector, giving the frequency at which each base occurs at each site, and the rate matrix (Tavare, 1986). When enough data (i.e. sequences) are available, the GTR model is superior over the more simplified Tamura 92 model. However, the Tamura 92 model is recommended in cases in which the data are not sufficient for reliable estimation of the model parameters and thus it is the default option for analyzing nucleic acid sequences in ConSurf. ConSurf also allows automatic selection of the model that best fits the analyzed sequences, as determined by the Akaike information criterion (AIC) (Akaike, 1974; Darriba et al., 2012, 2011).

Model of substitution for proteins

The inference of evolutionary conservation relies on a specified probabilistic model of amino-acid replacements (Pupko et al., 2002; Mayrose et al., 2004). The server supports a few models of substitution for nuclear DNA-encoded proteins as well as models of non-nuclear DNA-encoded proteins. The model of substitution can be chosen from the "Evolutionary Substitution Model" drop-down list. The LG (Le and Gascuel, 2008), JTT (Jones et al., 1992), Dayhoff (M.O. Dayhoff, R.M. Schwartz, and B. C, 1978) and WAG (Whelan and Goldman, 2001) matrices are suited for nuclear DNA-encoded proteins. The WAG matrix has been inferred from a large database of sequences comprising a broad range of protein families and is thus suited for distantly related amino acid sequences (Whelan and Goldman, 2001). The LG matrix (Le and Gascuel, 2008), incorporates variability of evolutionary rates across sites in the matrix, and was shown to outperform other substitutions matrices for proteins. The mtREV (Adachi and Hasegawa, 1996) and cpREV (Adachi et al., 2000) matrices are suitable for mitochondrial, and chloroplast DNA-encoded proteins, respectively. As in the case of analyzing nucleic acid sequences, the model that best fits the analyzed sequences can be automatically selected by ConSurf.

Conservation Scores

The conservation scores calculated by ConSurf appear in the SCORE column in the "Per-Residue Details" link in the outout page. The scores are normalized, so that the average score for all residues is zero, and the standard deviation is one. The conservation scores calculated by ConSurf are a relative measure of evolutionary conservation at each sequence site of the target chain. The lowest score represents the most conserved position in a protein (DNA/RNA). It does not necessarily indicate 100% conservation (e.g. no mutations at all), but rather indicates that this position is the most conserved in this specific protein (DNA/RNA) calculated using a specific MSA.

Confidence interval of the inferred conservation scores

Positions in the MSA exhibiting too little variation caused by too few sequences or too little diversity among sequence homologs can render evolutionary analysis meaningless (Thornton et al., 2000). Using the Bayesian method to calculate evolutionary conservation, confidence intervals for the conservation scores estimations are obtained (Susko et al., 2002). When the number of sequences is small, the confidence interval tends to be large, meaning a low level of support to the inferred conservation score. When the number of sequences increases, confidence intervals become smaller, and the point score estimates are more assured. In ConSurf, a confidence interval is assigned to each of the inferred evolutionary conservation scores. The confidence interval is defined by the lower and upper quartiles (the 25th and 75th percentiles of the inferred evolutionary rate distribution, respectively). This measure gives the 50% confidence interval and also indicates on the dispersion of each of the estimated scores. Amino acid positions that are assigned confidence intervals that are too large to be trustworthy are marked as having insufficient data in the output files of the server.

Coloring Scheme

The continuous conservation scores are partitioned into a discrete scale of 9 bins for visualization, such that bin 9 contains the most conserved positions and bin 1 contains the most variable positions. The color grades (1-9) are assigned as follows: The conservation scores below the average (negative values, which are indicative of slowly evolving, conserved sites) are divided into 4.5 equal intervals. The same 4.5 intervals are used for the scores above the average (positive values, which are indicative of rapidly evolving, variable sites). Thus, 9 equally sized categories of conservation are obtained. However, because the conservation distribution is asymmetrical around the average, the range of grade 1 is extended to include the most variable grades. Namely, the range of grade 1 is larger than the rest. Colors are then assigned to the 9 grades for graphic visualization. The width of each color grade varies for different polypeptide chains using this procedure. That is, the coloring results of a ConSurf run do not indicate the absolute magnitudes of evolutionary distances, but rather the relative degree of conservation of each amino acid position. ConSurf scaling procedure does not guarantee that grades 1-8 will always be occupied, although grade 9 is always occupied by at least one residue.

Dealing with unreliable positions

Conservation scores that are obtained for positions in the alignment that have less than 6 un-gapped amino acids are considered to be unreliable. When using the Bayesian method (Mayrose et al., 2004) for the conservation scores calculations, confidence intervals around the estimated rates are computed. The high and low values of each interval are assigned color grades according to the 1-9 coloring scheme. If the interval in a specific position spans 4 or more color grades the score is considered as unreliable. Such positions are colored light yellow in the graphic visualization output.

Output

If a 3D structure of the protein (or of the DNA/RNA) is provided:

1. The nine-color conservation scores are projected onto the 3D structure of the query protein and the colored protein structure is shown by FirstGlance in Jmol or NGL.
2. Scripts for visualizing the protein colored with ConSurf scores are generated for PyMol; (DeLano, 2008), UCSF-ChimeraX (Pettersen et al., 2021), Jmol; (Herráez, 2006).

For proteins, ConSurf also automatically suggests the possibility to map the calculated evolutionary conservation grades of the amino acids not only onto the single asymmetric unit chain, which is often deposited in the PDB, but also on the protein assembly as predicted by PISA

For proteins in which the 3D structure was not provided by the user, an up-to-date version of the Protein Data Bank (Rose et al., 2015) is searched for relevant homologues. If a structure of at least one homologous protein is available, the user may map the conservation scores on the structure. This option should ease the procedure for the non-expert users, who may be unfamiliar with the 3D structure. In addition, ConSurf goes one step further and uses HHPred (Meier and Söding, 2015) and MODELLER (Sali and Blundell, 1993) to produce a homology-model of the query (using default parameters). Briefly, HHPred uses a hidden Markov model (HMM) to search for potential templates of known 3D structure in the PDB (Remmert et al., 2012), and the MODELLER algorithm (Sali and Blundell, 1993) uses them to predict a 3D model for the query sequence. The ConSurf conservation grades are then mapped onto the predicted model.

For RNA sequence queries, ConSurf offers the possibility to predict the secondary structure. Structures are predicted using the RNAfold program of the Vienna package (Lorenz et al., 2011; Mathews et al., 2004), and the structure with the lowest free energy is selected. The ConSurf conservation grades are mapped onto the predicted secondary structure. Correlating the evolutionary data with the structural model offers the means to quickly detect functional regions within the RNA query.

For all cases, ConSurf creates the following outputs:

1. The sequence and MSA colored by ConSurf conservation scores.
2. A text file that summarizes for each position the normalized score calculated, the assigned grade (1-through-9), the reliability estimation (for the Bayesian method) and the amino acids/nucleotides observed in the respective MSA column.
3. The sequences selected for the MSA and the MSA constructed (unless those files were uploaded by the user).
4. A file with the frequency of each amino acid/nucleotide observed in each column of the MSA.
5. The evolutionary tree, which was calculated by the server or uploaded by the user, together with the MSA are shown using the WASABI platform (Veidenberg et al., 2015). Moreover, using WASABI, a user can select a subtree containing a fraction of the homologous sequences and conduct a follow-up ConSurf analysis with these selected sequences. To refine a ConSurf analysis to a selected subtree, a user can choose any internal node on the phylogenetic tree and open a WASABI menu by using a right mouse click. Selecting the option 'run ConSurf on subtree' will issue a new window with the new ConSurf run for the selected subtree sequences (see example in Figure 3).

Figure 3: Screenshot of WASABI view of the MSA and phylogenetic tree on the Consurf web server. The option to select sub tree and use it for a new ConSurf run is demonstrated. The subtree is marked by the red circle, and the selected sequences are highlighted with the black rectangle.

Example

As an example we provide the main output of a ConSurf run for the N-terminal region of the GAL4 transcription factor in yeast (PDB ID: 3COQ, chain A and B) in complex with its DNA recognition site( Figure 4). The analysis revealed, as expected, that the functional regions of this protein are highly conserved. For example, all the cysteines that form the Zn(2)-C6 DNA binding domain (CYS11, CYS14, CYS21, CYS28, CYS31, CYS38;(Pan and Coleman, 1990)) were assigned the highest conservation scores. Likewise, PRO26, which is known to be central for DNA binding (Johnston, 1987) is also highly conserved according to our analysis. In addition, other amino acid residues, which are in contact with the DNA (i.e. GLN9, LYS17, LYS18, LYS20, ARG15, LYS23; (Marmorstein et al., 1992)) are relatively conserved. ConSurf was also applied to nucleic acid sequences from yeast, which are the known binding sites of GAL4 and their adjacent neighborhood (Figure 4). As anticipated, the analysis revealed that the consensus pattern CGG-N11-CCG typical to GAL4 binding site is highly conserved. An extended full ConSurf analysis of this example is available in the 'GALLERY' section the ConSurf web site.

Figure 4: A ConSurf analysis for the GAL4 transcription factor and its DNA binding site. The 3D structure of the N-terminal region of the GAL4 transcription factor in yeast bound to the DNA is presented using a space-filled model. The amino-acids and the nucleotides are colored by their conservation grades using the color-coding bar, with turquoise-through-maroon indicating variable-through-conserved. Positions, for which the inferred conservation level was assigned with low confidence, are marked with light yellow. The figure reveals that the functionally important regions on both the DNA and the protein are highly conserved. The run was carried out using PDB code 3COQ and the figure was generated using the PyMol (DeLano, 2008) script output by ConSurf.

ConSurf vs estimates of the evolutionary conservation without explicit use of phylogenetic relationships

We conducted simulations to examine the accuracy of ConSurf's estimate of the evolutionary rates of the amino acid sites in comparison to entropy-based methods of estimating conservation/divergence, which do not take explicitly into account the phylogenetic relationships among the homologues. We synthesized a benchmark set of multiple sequence alignment, prepared using preset evolutionary rates, and examined the capacity of Rate4Site (Mayrose et al., 2004), the algorithm used in ConSurf, to reproduce these rates. In comparison, we examined the accuracy or reconstructing the rates using entropy based methods of estimating the evolutionary conservation, implemented by (Capra and Singh, 2007).

We compiled a set of multiple sequence alignments (MSAs) of homologous proteins by simulating the evolutionary process using INDELible (Fletcher and Yang, 2009) and preset evolutionary trees. Our set included 25 replicas of 16, 32, 64, 128, 256, and 512 taxa, i.e., a total of 150 trees. These ultrametric trees were generated using Mesquite (Maddison and Maddison, 2015), following a birth-death process with the default birth rate of 0.3 and death rate of 0.1. The tree height of each of these trees was then adjusted to be 0.6. Protein MSA was simulated for each phylogenetic tree by INDELible, using the following parameters: with amino-acid replacements following the LG matrix; Power law insertion/deletion length distribution (a=1.7) with maximum indel length of 30 amino acids; insertion rate which equals to deletion rate and equals to 0.01 relative to average substitution rate of 1; evolutionary rates for each position were drawn from 16 categories discrete gamma distribution with alpha=0.3; The length of the root sequence was set to be 500 amino acids. Overall, these parameters resulted in MSAs which looks biologically reasonable based on visual comparison of the alignments' total length, number and length of indels. The average length of the simulated MSAs were 567.32, 591.92, 663.08, 776.8, 984.64, 1295.92 amino acids for the 16, 32, 64 and 128, 256 and 512 taxa datasets, respectively. In addition we simulated sequences with rate equals to the Shannon entropy scores calculated for the default ConSurf MSA generated for the influenza neuraminidase (PDB ID: 2HU4 chain A) and using phylogenetic tree inferred from the MSA using RAxML. The simulated dataset and all calculated results can be downloaded here.

Each simulated MSA was next provided as input for each of the conservation inference methodologies. It is noteworthy that in rate4site, the phylogenetic tree was inferred from the MSA (rather than using the tree that was used in the simulations). Spearman correlation between the true evolutionary rate used during the simulation of each position and the inferred conservation score was calculated for each of the methods. The results are shown in table 1.

The data in the table 1 shows that Rate4Site is far superior to all examined alternatives; the Spearman correlations between the true (simulated) rates and those inferred by Rate4Site are significantly higher than these of the other entropy based methods. The last row in the table includes comparison to rates derived from a set of highly similar sequences (the influenza neuraminidase), as a particularly difficult test case. Rate4Site is slightly more accurate than the alternatives also in this case.