Sequence context-specific profiles for homology searching

Results

We first show that amino acid substitution scores are directly related to pairwise amino acid mutation probabilities and sequence-profile pseudocounts. We can therefore derive sequence context-specific amino acid similarity scores from context-specific mutation probabilities. These mutation probabilities can be predicted with a probabilistic model by using a large library of sequence-profile windows representing very specific local sequence contexts.

Any matrix of substitution scores S(x, y) describing the similarity between amino acids x and y can be written in the form (31) S(x, y) = const × log [P(x, y)/P(x)P(y)], where P(x, y) is the probability that x and y occur aligned to each other in an alignment of homologous sequences, and P(x) and P(y) are the background probabilities of x and y to occur in representative sequences (whether aligned or unaligned). This can also be written as a log odds score, S(x, y) = log [P(y|x)/P(y)], where P(y|x) = P(x, y)/P(x) is the conditional probability of y given x, i.e., the probability for amino acid x to mutate into y. If y occurs more often in positions aligned with an x (described by P(y|x)) than what would be expected by chance (described by P(y)), then the score is positive, otherwise negative.

We next explore the connection of mutation probabilities P(y|x) with sequence-profile pseudocounts. A sequence profile is a matrix p(i, y) that succinctly represents a multiple alignment of homologous sequences: p(i, y) is the frequency of amino acid y in column i of the multiple alignment. The profile describes what amino acids are likely to occur in related sequences at each position, or, in other words, the probability of a residue at position i to mutate into amino acid y. A single sequence (x_i) can be turned into a sequence profile by adding artificial mutations (i.e., pseudocounts) with the method of substitution matrix pseudocounts (11, 12): p(i, y) = P(y|x_i). Here, P(y|x_i) are the conditional probabilities giving rise to substitution matrix S(x, y). The profile-to-sequence score of column i of this single-sequence profile p with residue y_j of a sequence (y_j) is

Hence, substitution matrix scores can be seen as a special case of profile-to-sequence scores, where the profile is generated from one of the sequences by using substitution matrix pseudocounts.

Fig. 1A illustrates the equivalence of sequence-to-sequence and profile-to-sequences scoring with the alignment matrix of 2 zinc-finger sequences (x_i) and (y_j). The query profile resulting from the artificial mutations is illustrated as a histogram, in which the bar heights are proportional to the corresponding amino acid probabilities p(i, y). The score of each matrix cell (i, j) can be interpreted in 2 ways: either as sequence-to-sequence score S(x_i, y_j) between residues x_i and y_j, or as profile-to-sequence score S(p(i,·),y_j) between profile column p(i,·) and residue y_j.

In the above schemes, the expected mutation probabilities P(y/x_i) at position i depend only on the single amino acid x_i. However, the sequence context X_i, defined below, contains much more information than just residue x_i itself about what amino acids to expect in related sequences. If we were able to calculate a context-specific mutation probability P(y|X_i), we could define a score in a way analogous to Eq. 1, but by using a context-specific profile p_cs(i, y) = P(y|X_i) instead of P(y|x_i).

The context X_i is defined as the window of l residues surrounding x_i, i.e., X_i = (x_i−d, …, x_i+d) with l = 2d + 1. To predict the mutation probabilities for each position i, we compare its sequence window X_i with a precomputed library of K context profiles, p₁, …, p_K, each of length l. The context-specific mutation probability P(y|X_i), i.e., the probability of observing amino acid y in a homologous sequence given context X_i, will be calculated by a weighted mixing of the amino acids in the central columns of the most similar context profiles (Fig. 1B). To derive the weight of each profile p_k, we first need the probability P(X_i|p_k) that the sequence window X_i is emitted by profile p_k, which is equal to the product of probabilities for x_i+j (j ∈ {−d, …, d}) being emitted by profile column p_k(j,·): P(X_i|p_k) = ∏_{j = −d}^d p_k(j, x_i+j). Because the inner positions in the window will be most informative to predict the amino acid distribution for the central residue, we can refine the above formula by defining coefficients w_j, which weight the contribution of each window position: P(X_i|p_k) ∝ ∏_{j = −d}^d p_k(j, x_i+j)^w_j. The values of w_j are parameterized by w_center and β (see Fig. 1C). (For i within d residues from either end of (x_i) the product runs only over those j for which x_i+j is defined.)

Next, we need to know the probability P(p_k|X_i) that profile p_k was the one that emitted X_i. Using Bayes' theorem, we find

P(p_k) is the Bayesian prior probability for profile p_k, determined in the process of computing the profile library [supporting information (SI) Appendix]. It quantifies the probability that a sequence window is emitted by profile p_k prior to knowing that sequence window. P(X_i) = Σ_k P(X_i|p_k) P(p_k) is a normalization constant.

We can now calculate the context-specific mutation probabilities P(y|X_i) by mixing the amino acid distributions p_k(0, y) from the central columns of all K profiles with weights P(p_k|X_i):

Normalizing over all 20 amino acids yields the expected mutation probability P(y|X_i). To have more flexibility in adjusting the diversity of the context-specific profile p_cs(i,·), we mutate only a fraction τ ∈ [0,1] of (x_i) while leaving a fraction 1 − τ unchanged:

Here, δ_{xi, y} = 1 if x_i = y and 0 otherwise. In principle, τ needs to be optimized depending on the evolutionary distance over which homologous sequences are to be found, in a similar way as the substitution matrix with optimum diversity might be chosen. In practice, we have found that, as with substitution matrices, a single diversity works well for the entire range of evolutionary distances (SI Appendix).

Fig. 1B illustrates the calculation of expected mutation probabilities P(y|X_i) for a cysteine residue (highlighted in yellow) at position i belonging to a zinc-finger motif. Three profiles similar to the sequence window X_i (red box) are shown, whose central columns contribute to the context-specific sequence profile p(i, y) = P(y|X_i) at position i with weights P(p_k|X_i) of 7%, 60%, and 3%, respectively. With the resulting profile (Lower), a profile-to-sequence search can be performed, e.g., by using PSI-BLAST, which is equivalent to a sequence search with context-specific amino acid similarity scores (Eq. 1). In this example, the context-specific scheme recognizes the sequence context of the cysteine and correctly assigns a zinc-finger profile a high weight, resulting in a highly conserved cysteine.

The context-specificity paradigm is not restricted to sequences but applies equally well to sequence profiles or profile hidden Markov models (HMMs) (Materials and Methods). It can therefore be used in profile-to-sequence (9, 32, 33) and profile-to-profile (8, 10, 34–37) comparison, for example.

Our method CS-BLAST for context-specific protein sequence searching is a simple extension of BLAST. First, a context-specific sequence profile is generated for the query sequence as described. This step is very fast. Then PSI-BLAST is jump-started with this profile. PSI-BLAST is extended to the context-specific case in an analogous way (CSI-BLAST) (Materials and Methods).

Benchmark

The homology detection performance of our context-specific method CS-BLAST and standard NCBI BLAST is evaluated on a benchmark dataset derived from SCOP 1.73 (38), filtered to a maximum pairwise sequence identity of 20% (SCOP20, 6,616 domains). SCOP is a database of protein domains with known structure, hierarchically ordered by class, fold, superfamily, and family. Following a standard procedure, we consider all domains from the same superfamily to be homologous (true positives) and all pairs from different SCOP folds to be nonhomologous (false positives). Domain pairs from the same fold but different superfamilies are ignored.

We randomly assign members of every fifth fold in SCOP20 to the optimization set (1,329 domains), the others to the test set (5,287 domains). By using the optimization set, we determined the best values for the pseudocount admixture (τ = 0.9) and the window weights (w_center = 1.6, β = 0.85). The values for the window length (l = 13) and the context library size (K = 4,000) are a trade-off between sensitivity and time efficiency (see SI Appendix).

We perform an all-against-all comparison of the test-set domains and count the true and false positive hits at various E-value thresholds (Fig. 2A). To avoid a few large families from dominating the benchmark, we weight each true and false positive pair with 1/(size of SCOP family of first domain). Compared with NCBI BLAST (version 2.2.19, BLOSUM62, default parameters), CS-BLAST detects 139% more homologs at a cumulative error rate of 20%, 138% more at 10%, and, for the easiest cases at 1% error rate, still 96% more. To get an idea of the upside potential when parameters are trained on a larger set, we optimized w_center, β, and τ directly on the test set (red broken trace). These parameters (w_center = 1.3, β = 0.9, and τ = 0.95) are used in the official version of CS-BLAST.

To assess the alignment quality, we compare predicted sequence alignments to gold-standard structural alignments generated by TM-align (39). We start by randomly picking up to 10 domain pairs from each family in SCOP 1.73, requiring a maximum sequence identity of 30%, and aligning each pair with TM-align. Those domains that are not well superposable (TM-align score < 0.6) are discarded. This results in 11,457 domain pairs from 5,747 different domains. With each of the 5747 domains we perform a CS-BLAST and NCBI BLAST search against a database consisting of all domains belonging to the same family as the query domain and evaluate the quality of the predicted alignments for those pairs with a structural reference alignment. The alignment quality is assessed by 2 standard performance measures: Alignment sensitivity is the fraction of structurally aligned residue pairs that are correctly predicted, i.e., pairs correctly aligned/pairs structurally alignable. Alignment precision is defined as the fraction of aligned residue pairs in the predicted alignment that are correct, i.e., pairs correctly aligned/pairs aligned. Fig. 2B plots alignment sensitivity and precision for various sequence identity bins. CS-BLAST is able to improve the BLAST results over the entire range of sequence identities, especially for the difficult alignments. Very similar results are obtained when reference alignments are generated with DALI (40).

Another critical aspect for database search tools is the reliability of the reported E-values. The E-value of a match is an estimate of the number of chance hits to be expected with a score better than that of the database match. We check the reliability of CS-BLAST E-values by using the all-against-all searches of Fig. 2A. We count the number of false positives at a given E-value threshold, which, together with the size of the benchmark database, allows us to derive the actual E-value. Fig. 2C plots the actual against the reported E-value. NCBI BLAST's reported E-values are nearly identical to the observed ones. CS-BLAST E-values are too optimistic by a factor of ≈3–5, e.g., a reported E-value of 10⁻³ corresponds to an E-value of 5 × 10⁻³. Considering that this deviation is quite small and that it changes little with E-value, it should be easy to accommodate in practice.

Finally, we evaluate the homology detection performance of CSI-BLAST, the context-specific version of PSI-BLAST, on the benchmark of Fig. 2A. Because, typically, PSI-BLAST searches are done with a large sequence database, such as the nonredundant protein database (NR) at NCBI (41), to build diverse profiles, only the last search is performed against our benchmark database; all previous iterations use the full NR database (E-value inclusion thresholds set to 1 × 10⁻³ for PSI-BLAST and 2 × 10⁻⁴ for CSI-BLAST). Fig. 2D plots true positives versus false positives detected by PSI-BLAST and CSI-BLAST after up to 5 search iterations. (The trace for CSI-BLAST with 5 iterations has been omitted because it does not significantly improve over 3 iterations anymore. The traces for one iteration are the same as in Fig. 2A.) Remarkably, 2 iterations of CSI-BLAST are more sensitive than 5 rounds of standard PSI-BLAST (≈15% more homologs detected). This result surprised us. We had expected that context-specific profiles would only marginally improve sensitivity over standard sequence profiles, because profiles already contain family- and position-specific mutation rates. However, the lead of CS-BLAST over BLAST is even extended in absolute terms after the second iteration, demonstrating that the context profiles contain local information from analogous sequences (i.e., with similar sequence context) that is partially independent of information from the homologous sequences in the profile.

Example: Activation Domain of SOX-9

Fig. 1B gave an example in which the context-specific method led to above-average conservation of Zn-finger cysteines. In practice, it will be equally important to be able to guess which residues are conserved less than average. As an example, Fig. 3 presents profiles of a region from the activation domain of human SOX-9 transcription factor, generated with substitution matrix pseudocounts (Left) and context-specific pseudocounts (Right). Because this region is natively disordered, its sequence is only very weakly conserved. The substitution matrix method assigns the same amino acid distribution to the prolines as it would to a proline in a globular domain. The context-specific method, however, mixes the pseudocounts mainly from contexts that are also disordered, weakly conserved, and have a similar, biased amino acid distribution. Therefore, its profile exhibits below-average conservation of prolines, alanines, and glutamines while having higher overall probabilities for these residues.

Discussion

Sequence context is much more powerful than a single residue in predicting which amino acids that particular residue is likely to mutate into (Fig. 2 A and B). Because this context information is as easy to get as the sequence itself, it is surprising that sequence context is practically never exploited. The main reason seems to be the focus of past research on structural context, with its limitation to proteins of known structures (17, 18, 20, 21, 27). Another reason may be the challenge to develop sequence context-specific methods that can compete with traditional context-free methods such as BLAST and PSI-BLAST in speed and usability (26, 27). We have shown how context-specific pseudocounts can be used in combination with existing profile-based methods to extend residue-centered sequence comparison to the context-specific case, without loss of speed or usability.

As examples, we have built context-specific versions of BLAST and PSI-BLAST that considerably improve their performance at very little runtime overhead. For a typical protein of length L = 250 and a library size of K = 4,000, the computation of the context-specific profile requires ≈1 s CPU time. Also, runtime scales favorably, T ∝ KlL (SI Appendix, Fig. S1). (Note that HMMSUM's runtime scales as T ∝ K²L, which places a strict practical limit on the number K of states/contexts in HMMSUM.) Because the output of CS-BLAST and CSI-BLAST is generated by the BLAST and PSI-BLAST programs themselves, users do not have to get accustomed to different command line options or output formats, and updates to the BLAST package will directly benefit the context-specific versions. The only caveat is that E-values need to be corrected by a factor of 3 to 5 (Fig. 2C). We expect CS-BLAST to be useful to find homologs for singleton sequences, because, for these, the lack of homologs precludes the use of profile-to-sequence search methods such as PSI-BLAST.

A pleasant surprise is the extent of improvements of sequence profiles through context-specific pseudocounts (Fig. 2D), even though profiles already contain evolutionary information on position- and family-specific mutation probabilities. Hence, the information from locally similar, analogous sequences that are contained in the context profiles is at least partly orthogonal to the evolutionary information in the homologous sequences that contribute to the sequence profiles. Consequently, we can expect improvements when applying the new paradigm to the pairwise comparison of sequence profiles (34–36) and profile HMMs (10, 37), or to hierarchical multiple sequence alignment programs (8, 42).

It is possible to extend Dirichlet mixture pseudocounts (29, 30) to the context-specific case. This would yield an alternative formulation of context-specific sequence comparison that is worth exploring. In that scheme, the context library would have K metaprofiles, i.e., multicolumn pseudocount priors. Each metaprofile would consist of l Dirichlet distributions and would be able to emit a profile with l columns. An advantage over the presented scheme might be that the diversity of each column in the metaprofiles is encoded by one additional parameter per column (the sum of all pseudocounts in a column), which might lead to better modeling of the profile contexts.

The paradigm presented here should be easily transferable to nucleotide sequences. The application to noncoding regions such as promoter regions and regions harboring putative noncoding RNAs (ncRNAs) is of particular interest. The low information content of nucleotide sequences and the often weak overall conservation in these regions render alignments between related species difficult, whereas reliable alignments offer enormous potential to identify functional regions (such as cis-regulatory elements or ncRNAs) through their interspecies conservation (see, e.g., ref. 43).

In summary, the paradigm of sequence context specificity offers greatly improved sensitivity and alignment quality in protein sequence comparison and is likely to hold similar advantages for nucleotide sequences. We believe that these advantages are sufficient to warrant a paradigm shift in biological sequence comparison, alignment, and molecular evolution from amino acid- and nucleotide-centric to context-specific methods.

Materials and Methods

Generalization to Sequence Profiles.

To apply the paradigm to sequence profiles and profile HMMs, we show how to generalize the calculation of pseudocounts from the single sequence case in Eq. 2 to the case of sequence alignments, from which the profile is derived. In analogy to the sequence context X_i, we define the context of the query alignment at position i as Q_i = (c_q(i − d,·), …, c_q(i + d,·)), where c_q(j, x) are the counts of amino acid x at position j of the query alignment. These counts are obtained from the sequence profile q(j, x) by multiplying with the effective number of sequences N_q(j) at position j in the query alignment: c_q(j, x) = N_q(j) q(j, x) (see SI Appendix for details). We now merely need to show how to generalize P(X_i|p_k) to P(Q_i|p_k), because all other transformations leading to Eq. 2 remain essentially unchanged. To derive P(Q_i|p_k), we model the amino acid counts c_q(i) with multinomial distributions. Because N_q(j) can be real-valued, however, we replace the factorials in the multinomial distribution by Gamma functions (n! = Γ(n + 1))

Note that, because the factor containing the Gamma functions does not depend on k, it will cancel out during the normalization of P(p_k|Q_i) (Eq. 2, SI Appendix, Eq. S9). Similar to PSI-BLAST (9), we choose the pseudocount admixture τ in Eq. 3 depending on the diversity of the query alignment, τ = a(b + 1)/(b + N_q(i)), where a = 0.9 and b = 12.0 have been determined on the training set as described in SI Appendix.

Supplementary Material