Abstract Summary: We introduce an algorithm that uses the information gained from simultaneous consideration of an entire group of related proteins to create multiple structure alignments (MSTAs). Consistency-based alignment (CBA) first harnesses the information contained within regions that are consistently aligned among a set of pairwise superpositions in order to realign pairs of proteins through both global and local refinement methods. It then constructs a multiple alignment that is maximally consistent with the improved pairwise alignments. We validate CBA's alignments by assessing their accuracy in regions where at least two of the aligned structures contain the same conserved sequence motif. Results: CBA correctly aligns well over 90% of motif residues in superpositions of proteins belonging to the same family or superfamily, and it outperforms a number of previously reported MSTA algorithms. Availability: CBA is available at and the source code is freely available at Contact: [email protected] Supplementary information: Supplementary data are available at Bioinformatics online.

Protein interactions typically arise from a physical interaction of one or more small sites on the surface of the two proteins. Identifying these sites is very important for drug and protein design. In this paper, we propose a computational method based on probabilistic relational model that at- tempts to address this task using high-throughput protein interaction data and a set of short sequence motifs. We learn the model using the EM algorithm, with a branch-and-bound algorithm as an approximate infer- ence for the E-step. Our method searches for motifs whose presence in a pair of interacting proteins can explain their observed interaction. It also tries to determine which motif pairs have high affinity, and can therefore lead to an interaction. We show that our method is more accurate than others at predicting new protein-protein interactions. More importantly, by examining solved structures of protein complexes, we find that 2/3 of the predicted active motifs correspond to actual interaction sites. 1 Introduction Many cellular functions are carried out through physical interactions between proteins. Discovering the protein interaction map can therefore help to better understand the work- ings of the cell. Indeed, there has been much work recently on developing high-throughput methods to produce a more complete map of protein-protein interactions [1, 2, 3]. Interactions between two proteins arise from physical interactions between small re- gions on the surface of the proteins [4] (see Fig. 2(b)). Finding interaction sites is an important task, which is of particular relevance to drug design. There is currently no high- throughput experimental method to achieve this goal, so computational methods are re- quired. Existing methods either require solving a protein's 3D structure (e.g., [5]), and therefore are computationally very costly and not applicable on a genome-wide scale, or use known interaction sites as training data (e.g., [6]), which are relatively scarce and hence have poor coverage. Other work focuses on refining the highly noisy high-throughput inter- action maps [7, 8, 9], or on assessing the confidence levels of the observed interactions [10]. In this paper, we propose a computational method for predicting protein interactions P1 P2 P5 d a d A A A A A P a d b d d 5 P1 A A A A A ab db ad dd bd B S B B S Bab B S ab ab b ab db c d b T T T P P4 12 15 25 2 P O O O 3 b (a) (b) Figure 1: (a) Simple illustration of our assumptions for protein-protein interactions. The small elements denote motif occurrences on proteins, with red denoting active and gray denoting inactive motifs. (b) A fragment of our probabilistic model, for the proteins P1, P2, P5. We use yellow to denote an assignment of the value true, and black to denote the value false; full circles denote an assignment observed in the data, and patterned circles an assignment hypothesized by our algorithm. The dependencies involving inactive motif pairs were removed from the graph because they do not affect the rest of the model. and the sites at which the interactions take place, which uses as input only high-throughput protein-protein interaction data and the protein sequences. In particular, our method as- sumes no knowledge of the 3D protein structure, or of the sites at which binding occurs. Our approach is based on the assumption that interaction sites can be described using a limited repertoire of conserved sequence motifs [11]. This is a reasonable assumption since interaction sites are significantly more conserved than the rest of the protein surface [12]. Given a protein interaction map, our method tries to explain the observed interactions by identifying a set of sites of motif occurrence on every pair of interacting proteins through which the interaction is mediated. To understand the intuition behind our approach, con- sider the example of Fig. 1(a). Here, the interaction pattern of the protein P1 can best be explained using the motif pair a, b, where a appears in P1 and b in the proteins P2, P3, P4 but not in P5. By contrast, the motif pair d, b is not as good an explanation, because d also appears in P5, which has a different interaction pattern. In general, our method aims to identify motif pairs that have high affinity, potentially leading to interaction between protein pairs that contain them. However, a sequence motif might be used for a different purpose, and not give rise to an active binding site; it might also be buried inside the protein, and thus be inaccessible for interaction. Thus, the appearance of an appropriate motif does not always imply interaction. A key feature of our approach is that we allow each motif occurrence in a protein to be either active or inactive. Interactions are then induced only by the interactions of high- affinity active motifs in the two proteins. Thus, in our example, the motif d in p2 is inactive, and hence does not lead to an interaction between p2 and p4, despite the affinity between the motif pair c, d. We note that Deng et al. [8] proposed a somewhat related method for genome-wide analysis of protein interaction data, based on protein domains. However, their method is focused on predicting protein-protein interactions and not on revealing the site of interaction, and they do not allow for the possibility that some domains are inactive. Our goal is thus to identify two components: the affinities between pairs of motifs, and the activity of the occurrences of motifs in different proteins. Our algorithm addresses this problem by using the framework of Bayesian networks [13] and probabilistic relational models [14], which allows us to handle the inherent noise in the protein interaction data and the uncertain relationship between interactions and motif pairs. We construct a model encoding our assumption that protein interactions are induced by the interactions of active motif pairs. We then use the EM algorithm [15], to fill in the details of the model, learning both the motif affinities and activities from the observed data of protein-protein interactions and protein motif occurrences. We address the computational complexity of the E-step in these large, densely connected models by using an approximate inference procedure based on branch-and-bound. We evaluated our model on protein-protein interactions in yeast and Prosite motifs [11]. As a basic performance measure, we evaluated the ability of our method to predict new protein-protein interactions, showing that it achieves better performance than several other models. In particular, our results validate our assumption that we can explain interactions via the interactions of active sequence motifs. More importantly, we analyze the ability of our method to discover the mechanism by which the interaction occurs. Finally, we examined co-crystallized protein pairs where the 3D structure of the interaction is known, so that we can determine the sites at which the interaction took place. We show that our active motifs are more likely to participate in interactions. 2 The Probabilistic Model The basic entities in our probabilistic model are the proteins and the set of sequence motifs that can mediate protein interactions. Our model therefore contains a set of protein entities P = {P1, . . . , Pn}, with the motifs that occur in them. Each protein P is associated with the set of motifs that occur in it, denoted by P.M . As we discussed, a key premise of our approach is that a specific occurrence of a sequence motif may or may not be active. Thus, each motif occurrence a P.M is associated with a binary-value variable P.Aa, which takes the value true if Aa is active in protein P and false otherwise. We structure the prior probability P (P.Aa = true) = min{0.8, 3+0.1|P.M| }, to capture our intuition that |P.M | the number of active motifs in a protein is roughly a constant fraction of the total number of motifs in the protein, but that even proteins with few motifs tend to have at least some number of active motifs. A pair of active motifs in two proteins can potentially bind and induce an interaction between the corresponding proteins. Thus, in our model, a pair of proteins interact if each contains an active motif, and this pair of motifs bind to each other. The probability with which two motifs bind to each other is called their affinity. We encode this assumption by including in our model entities Tij corresponding to a pair of proteins Pi, Pj. For each pair of motifs a Pi.M and b Pj.M , we introduce a variable Tij.Aab, which is a deterministic AND of the activity of these two motifs. Intuitively, this variable represents whether the pair of motifs can potentially interact. The probability with which two active motif occurrences bind is their affinity. We model the binding event between two motif occurrences using a variable Tij.Bab, and define: P (Tij.Bab = true | Tij.Aab = true) = ab and P (Tij.Bab = true | Tij.Aab = false) = 0, where ab is the affinity between motifs a and b. This model reflects our assumption that two motif occurrences can bind only if they are both active, but their actual binding probability depends on their affinity. Note that this affinity is a feature of the motif pair and does not depend on the proteins in which they appear. We must also account for interactions that are not explained by our set of motifs, whether because of false positives in the data, or because of inadequacies of our model or of our motif set. Thus, we add a spurious binding variable Tij.S, for cases where an interaction between Pi and Pj exists, but cannot be explained well by our set of active motifs. The probability that a spurious binding occurs is given by P (Tij.S = true) = S. Finally, we observe an interaction between two proteins if and only if some form of binding occurs, whether by a motif pair or a spurious binding. Thus, we define a variable Tij.O, which represents whether protein i was observed to interact with protein j, to be a deterministic OR of all the binding variables Tij.S and Tij.Bab. Overall, Tij.O is a noisy-OR [13] of all motif pair variables Tij.Aab. Note that our model accounts for both types of errors in the protein interaction data. False negatives (missing interactions) in the data are addressed through the fact that the presence of an active motif pair only implies that binding takes place with some probability. False positives (wrong interactions) in the data are addressed through the introduction of the spurious interaction variables. The full model defines a joint probability distribution over the entire set of attributes: P (P.A, T.A, T.B, T.S, T.O) = P (P i aP i.Aa) i .M P (T aP ij .Aab | Pi.Aa, Pj .Ab)P (Tij .Bab | Tij .Aab) i .M,bPj .M ij P (Tij.S)P (Tij.O | Tij.B, Tij.S) where each of these conditional probability distributions is as specified above. We use to denote the entire set of model parameters {a,b}a,b {S}. An instantiation of our probabilistic model is illustrated in Fig. 1(b). 3 Learning the Model We now turn to the task of learning the model from the data. In a typical setting, we are given as input a protein interaction data set, specifying a set of proteins P and a set of observed interacting pairs T.O. We are also given a set of potentially relevant motifs, and the occurrences of these motifs in the different proteins in P. Thus, all the variables except for the O variables are hidden. Our learning task is thus twofold: we need to infer the values of the hidden variables, both the activity variables P.A, T.A, and the binding variables T.B, T.S; we also need to find a setting of the model parameters , which specify the motif affinities. We use a variant of the EM algorithm [15] to find both an assignment to the parameters , and an assignment to the motif variables P.A, which is a local maximum of the likelihood function P (T.O, P.A | ). Note that, to maximize this objective, we search for a MAP assignment to the motif activity variables, but sum out over the other hidden variables. This design decision is reasonable in our setting, where determining motif activities is an important goal; it is a key assumption for our computational procedure. As in most applications of EM, our main difficulty arises in the E-step, where we need to compute the distribution over the hidden variables given the settings of the observed variables and the current parameter settings. In our model, any two motif variables (both within the same protein and across different proteins) are correlated, as there exists a path of influence between them in the underlying Bayesian network (see Fig. 1(c)). These cor- relations make the task of computing the posterior distribution over the hidden variables intractable, and we must resort to an approximate computation. While we could apply a general purpose approximate inference algorithm such as loopy belief propagation [16], such methods may not converge in densely connected model such as this one, and there are few guarantees on the quality of the results even if they do converge. Fortunately, our model turns out to have additional structure that we can exploit. We now describe an ap- proximate inference algorithm that is tailored to our model, and is guaranteed to converge to a (strong) local maximum. Our first observation is that the only variables that correlate the different protein pairs Tij are the motif variables P.A. Given an assignment to these activity variables, the net- work decomposes into a set of independent subnetworks, one for each protein pair. Based on this observation, we divide our computation of the E-step into two parts. In the first, we find an assignment to the motif variables in each protein, P.A; in the second, we com- pute the posterior probability over the binding motif pair variables T.B, T.S, given the assignment to the motif variables. We begin by describing the second phase. We observe that, as all the motif pair vari- ables, T.A, are fully determined by the motif variables, the only variables left to reason about are the binding variables T.B and T.S. The variables for any pair Tij are inde- pendent of the rest of the model given the instantiation to T.A and the interaction evi- dence. That fact, combined with the noisy-OR form of the interaction, allows us to com- pute the posterior probability required in the E-step exactly and efficiently. Specifically, the computation for the variables associated with a particular protein pair Tij is as fol- lows, where we omit the common prefix Tij to simplify notation. If Aab = false, then P (Bab = true | Aab = false, O, ) = 0. Otherwise, if Aab = true, then P (B P (B ab | A, )P (O | Bab = true, A, ) ab = true | A, O, ) = . P (O | A, ) The first term in the numerator is simply the motif affinity ab; the second term is 1 if O = true and 0 otherwise. The numerator can easily be computed as P (O | A, ) = 1 - (1 - S) (1 - A ab). The computation for P (S) is very similar. a,b =true We now turn to the first phase, of finding a setting to all of the motif variables. Un- fortunately, as we discussed, the model is highly interconnected, and a finding an optimal joint setting to all of these variables P.A is intractable. We thus approximate finding this joint assignment using a method that exploits our specific structure. Our method iterates over proteins, finding in each iteration the optimal assignment to the motif variables of each protein given the current assignment to the motif activities in the remaining proteins. The process repeats, iterating over proteins, until convergence. As we discussed, the likelihood of each assignment to Pi.A can be easily computed using the method described above. However, the computation for each protein is still ex- ponential in the number of motifs it contains, which can be large (e.g., 15). However, in our specific model, we can apply the following branch-and-bound algorithm (similar to an approach proposed by Henrion [17] for BN2O networks) to find the globally optimal as- signment to the motif variables of each protein. The idea is that we search over the space of possible assignments Pi.A for one that maximizes the objective we wish to maximize. We can show that if making a motif active relative to one assignment does not improve the objective, it will also not improve the objective relative to a large set of other assignments. More precisely, let f (Pi.A) = P (Pi.A, P-i.A|O, ) denote the objective we wish to maximize, where P-i.A is the fixed assignment to motif variables in all proteins except Pi. Let Pi.A-a denote the assignment to all the motif variables in Pi except for Aa. We compute the ratio of f after we switch Pi.Aa from false to true. Let ha(Pj) = (1 - P ab) denote the probability that motif a does not bind with j .Ab =true any active motif in Pj. We can now compute: f (P g i.Aa = true, Pi.A-a) a(Pi.A-a) = = f (Pi.Aa = false, Pi.A-a) 1 - g 1 - (1 - S)ha(Pj) hb(Pj) h a=b,Pi.Ab=true a(Pj ) (1) 1 - (1 - S) h a=b,P b(Pj ) 1jn 1jn i .Ab =true Tij .O=false Tij .O=true where g is the prior probability for a motif in protein Pi to be active. Now, consider a different point in the search, where our current motif activity assign- ment is Pi.A-a, which has all the active motifs in Pi.A-a and some additional ones. The first two terms in the product of Eq. (1) are the same for a(Pi.A-a) and a(Pi.A-a). For the final term (the large fraction), one can show using some algebraic manipulation that this term in a(Pi.A-a) is lower than that for a(Pi.A-a). We conclude that a(Pi.A-a) a(Pi.A-a), and hence that: f (Pi.Aa = true, Pi.A-a) f (P 1 i.Aa = true, Pi.A-a) 1. f (Pi.Aa = false, Pi.A-a) f (Pi.Aa = false, Pi.A- ) a It follows that, if switching motif a from inactive to active relative to Pi.A decreases f , it will also decrease f if we have some additional active motifs. We can exploit this property in a branch-and-bound algorithm in order to find the glob- ally optimal assignment Pi.A. Our algorithm keeps a set V of viable candidates for motif assignments. For presentation, we encode assignments via the set of active motifs they contain. Initially, V contains only the empty assignment {}. We start out by considering motif assignments with a single active motif. We put such an assignment {a} in V if its f -score is higher than f ({}). Now, we consider assignments {a, b} that have two active motifs. We consider {a, b} only if both {a} and {b} are in V . If so, we evaluate its f -score, and add it to V if this score is greater than that of {a} and {b}. Otherwise, we throw it away. We continue this process for all assignments of size k: For each assignment with active motif set S, we test whether S - {a} V for all a S; if we compare f (S) to each f (S - {a}), and add it if it dominates all of them. The algorithm terminates when, from some k, no assignment of size k is saved. To understand the intuition behind this pruning procedure, consider a candidate assign- ment {a, b, c, d}, and assume that {a, b, c} V , but {b, c, d} V . In this case, we must have that {b, c} V , but adding d to that assignment reduces the f -score. In this case, as shown by our analysis, adding d to the superset {a, b, c} would also reduce the f -score. This algorithm is still exponential in worst case. However, in our setting, a protein with many motifs has a low prior probability that each of them is active. Hence, adding new motifs is less likely to increase the f -score, and the algorithm tends to terminate quickly. As we show in Section 4, this algorithm significantly reduces the cost of our procedure. Our E-step finds an assignment to P.A which is a strong local optimum of the ob- jective function max P (P.A | T.O, ): The assignment has higher probability than any assignment that changes any of the motif variables for any single protein. For that assign- ment, our algorithm also computes the distribution over all of the binding variables, as described above. Using this completion, we can now easily compute the (expected) suffi- cient statistics for the different parameters in the model. As each of these parameters is a simple binomial distribution, the maximum likelihood estimation in the M-step is entirely standard; we omit details.

Abstract Summary: 3motif is a web application that visually maps conserved sequence motifs onto three-dimensional protein structures in the Protein Data Bank (PDB; Berman et al., Nucleic Acids Res., 28, 235–242, 2000). Important properties of motifs such as conservation strength and solvent accessible surface area at each position are visually represented on the structure using a variety of color shading schemes. Users can manipulate the displayed motifs using the freely available Chime plugin. Availability: http://motif.stanford.edu/3motif/ Contact: [email protected] Supplementary information: http://motif.stanford.edu/3motif/supplementary/ * To whom correspondence should be addressed. † Present adress: Google, Inc., 2400 Bayshore Parkway, Mountain View, CA~94043, USA

Classic techniques for simulating molecular motion, such as the Monte Carlo and molecular dynamics methods, generate individual motion pathways one at a time and spend most of their time trying to escape from the local minima of the energy landscape of a molecule. Their high computational cost prevents them from being used to analyze many pathways. We introduce Stochustic Roadmap Sirrrcllation (SRS), a new approach for exploring the kinetics of molecular motion by simultaneously examining multiple pathways encoded compactly in a graph, called a roadmap. A roadmap is computed by sampling a molecule's conformation space at random. The computation does not suffer from the localminima problem encountered with existing methods. Each path in the roadmap represents a potential motion pathway and is associated with a probability indicating the likelihood that the molecule follows this pathway. By viewing the roadmap as a Markov chain, we can efficiently compute properties of molecular motion over the entire molecular energy landscape. We also prove that, in the limit, SRS converges to the same distribution as Monte Carlo simulation. To test the effectiveness of our approach, we apply it to the computation of the transmission coefficients for protein folding, an important order parameter that measures the kinetic distance of a protein's conformation to its native state Our computational studies show that SRS obtains more accurate results and achieves several orders- of- magnitude reduction in computation time, compared with Monte Carlo simulatio.

Probabilistic roadmaps are an effective tool to compute the connectivity of the collision-free subset of high-dimensional robot configuration spaces. This paper extends them to capture the pertinent features of continuous functions over high-dimensional spaces. We focus here on computing energetically favorable motions of bio-molecules. A molecule is modeled as an articulated structure moving in an energy field. The set of all its 3D placements is the molecule's conformational space, over which the energy field is defined. A probabilistic conformational roadmap (PCR) tries to capture the connectivity of the low-energy subset of a conformational space, in the form of a network of weighted local pathways. The weight of a pathway measures the energetic difficulty for the molecule to move along it. The power of a PCR derives from its ability to compactly encode a large number of energetically favorable molecular pathways, each defined as a sequence of contiguous local pathways. This paper describes general techniques to compute and query PCRs, and presents implementations to study ligand-protein binding and protein folding.

MOTIVATION: We present techniques for increasing the speed of sequence analysis using scoring matrices. Our techniques are based on calculating, for a given scoring matrix, the quantile function, which assigns a probability, or p, value to each segmental score. Our techniques also permit the user to specify a p threshold to indicate the desired trade-off between sensitivity and speed for a particular sequence analysis. The resulting increase in speed should allow scoring matrices to be used more widely in large-scale sequencing and annotation projects. RESULTS: We develop three techniques for increasing the speed of sequence analysis: probability filtering, lookahead scoring, and permuted lookahead scoring. In probability filtering, we compute the score threshold that corresponds to the user-specified p threshold. We use the score threshold to limit the number of segments that are retained in the search process. In lookahead scoring, we test intermediate scores to determine whether they will possibly exceed the score threshold. In permuted lookahead scoring, we score each segment in a particular order designed to maximize the likelihood of early termination. Our two lookahead scoring techniques reduce substantially the number of residues that must be examined. The fraction of residues examined ranges from 62 to 6%, depending on the p threshold chosen by the user. These techniques permit sequence analysis with scoring matrices at speeds that are several times faster than existing programs. On a database of 12 177 alignment blocks, our techniques permit sequence analysis at a speed of 225 residues/s for a p threshold of 10-6, and 541 residues/s for a p threshold of 10-20. In order to compute the quantile function, we may use either an independence assumption or a Markov assumption. We measure the effect of first- and second-order Markov assumptions and find that they tend to raise the p value of segments, when compared with the independence assumption, by average ratios of 1.30 and 1.69, respectively. We also compare our technique with the empirical 99. 5th percentile scores compiled in the BLOCKSPLUS database, and find that they correspond on average to a p value of 1.5 x 10-5. AVAILABILITY: The techniques described above are implemented in a software package called EMATRIX. This package is available from the authors for free academic use or for licensed commercial use. The EMATRIX set of programs is also available on the Internet at http://motif.stanford.edu/ematrix.

A general framework is presented for analyzing multiple protein structures. A family of related protein structures may be analyzed using statistical regression methods. The analysis requires alternating steps of finding correspondences among the protein structures and superimposing the corresponding landmarks. The superposition step may be performed using either affine or orthogonal transformations, thereby allowing protein structures to undergo either pure rotations or rotation plus shear operations. Regression analysis permits a separate weight for each position, allowing one to emphasize particular segments of a protein structure or to compensate for variances that differ at various positions in a structure. In addition, a novel method is introduced for finding an initial correspondence, based on matching discrete curvatures along the protein backbone. Another novel method is introduced for obtaining gap functions that adapt to the given data, thereby making dynamic programming meth...

Medical informatics is defined largely by its host disciplines in clinical and biological medicine, and to project the agenda for informatics into the next decade, the health community must envision the broad context of biomedical research. This paper is a sketch of this vision, taking into account pressures from changes in the U.S. health care system, the need for more objective information on which to base health care decisions, and the accelerating progress and clinical impact of genomics research. The lessons of modern genomics research demonstrate the power of computing and communication tools to facilitate rapid progress through the adoption of open community standards for information exchange and collaboration. While aspects of this vision are speculative, it seems clear that the core agenda for informatics must be the development of interoperating systems that can facilitate the secure gathering, interchange, and analysis of high-quality information and can gain leverage from worldwide collaboration in advancing and applying new medical knowledge.

We have developed a knowledge-based simulation of DNA metabolism that accurately predicts the actions of enzymes on DNA under a large number of environmental conditions. Previous simulations of enzyme systems rely predominantly on mathematical models. We use a frame-based representation to model enzymes, substrates and conditions. Interactions between these objects are expressed using production rules and an underlying truth maintenance system. The system performs rapid inference and can explain its reasoning. A graphical interface provides access to all elements of the simulation, including object representations and explanation graphs. Predicting enzyme action is the first step in the development of a large knowledge base to envision the metabolic pathways of DNA replication and repair.

We have increased the sensitivity of DNA and protein sequence database searches by allowing similar but non-identical amino acids or nucleotides to match. In addition, one can match k-tuples or words instead of matching individual residues in order to speed the search. A matching matrix species which k-tuples match each other. The matching matrix can be calculated from a similarity matrix of amino acids and a threshold of similarity required for matching. This permits amino acid similarity matrices or replacement matrices (PAM matrices) to be used in the first step of a sequence comparison rather than in a secondary scoring phase. The concept of matching non-identical k-tuples also increases the power of DNA database searches. For example, a matrix that specifies that any 3-tuple in a DNA sequence can match any other 3-tuple encoding the same amino acid permits a DNA database search using a DNA query sequence for regions that would encode a similar amino acid sequence.