The goal of the bootstrap phase is to find an optimal alignment of small substructures of 3 nodes, called seeds. A seed is represented by a triple of residues A = (A₁, A₂, A₃).

The set of possible candidates for the initial alignment consists of all seeds satisfying the following conditions:

Feasible candidates are seeds satisfying both (a) and (b). If one or more proteins contain no feasible candidates, the search stops and a new iteration of PROPOSAL begins.

Once a set of suitable candidates is generated, PROPOSAL tries to construct an optimal initial alignment through Gibbs Sampling on top of a Monte Carlo Markov Chain (MCMC). In the MCMC each state represents an alignment of N seeds, one from each protein structure.

Starting from a random initial state (i.e., a random initial alignment), the sampling method iteratively performs a transition from a state of the chain to another, by replacing a randomly chosen seed of the current alignment with a feasible candidate of the same protein, according to a properly defined transition probability distribution. When Gibbs sampling stops, the last current alignment is returned. If the sampling procedure is iterated a sufficient number of times, it converges to a local optimum solution.

A critical task is to establish when Gibbs sampling can be stopped. The procedure ends when the alignment of seeds does not change. Let P=(N−1N)i be the probability that a protein structure is never selected in i consecutive iterations of Gibbs sampling. The number of iterations of Gibbs sampling is determined by the following parameter k: (2)k=max{k′:(N−1N)k′>α} where α is a user-defined probability threshold. If the alignment does not change for k consecutive iterations, the Gibbs sampling is stopped. The lower is α, the more precise and slower will be the sampling procedure. Therefore, α represents a trade-off between accuracy and speed of PROPOSAL.

The transition probability is defined on top of a similarity score, based on the distances between the residues of the seeds. Let Dist(R₁, R₂) be the euclidean distance between the two residues R₁ and R₂ of a 3D structure. Given two seeds A = (A₁, A₂, A₃) and B = (B₁, B₂, B₃), we define the pairwise distance between A and B as:

Now, let S = {S₁, S₂, …, S_N} be the alignment of seeds at the i-th iteration of Gibbs sampling and suppose we have to replace S_j by a feasible candidate X of the same protein. The similarity score of X is defined as the inverse of the product of all pair distances between X and the seeds of the current alignment (except S_j):

The transition probability is then computed by normalizing such similarity scores in [0,1].

In the extension phase, the alignment of residues is extended up to size w by iteratively adding N residues to the current alignment, one from each protein.

Suppose that we start from a substructure alignment of size w′<w. The goal is to find an optimal alignment of N residues R₁, R₂,…, R_N, one for each protein, and add such residues to the substructure alignment. R_i must be at distance at most equal to 10 Å from residues in the corresponding current aligned substructure. At the end of this process, the alignment size will be w′+1.

Each extension step is performed through a Gibbs sampling strategy similar to the one used during the bootstrap phase. In the extension phase the similarity score takes into account:

Let SA = {SA₁, SA₂, …, SA_N} be the current alignment of size w′, where each SA_i = {R_i,1, R_i,2,…, R_i,w′} is a set of residues, and let A^m = {A^m₁, A^m₂,…, A^m_N} be the alignment of candidate residues at the generic m-th iteration of Gibbs sampling.

Suppose we replace A^m_j with a candidate residue X. First, we define a similarity score, SimSymb(X) which evaluates the similarity between the symbol of X and the symbols of residues in A^m (except A^m_j): (5)SimSymb(X)=∏k=1,i≠jN SIMMATRIX(X,Akm) where SimMatrix(X,A^m_k) is a BLOSUM similarity score between X and A^m_k.

Then, we define another similarity function, SimDist(X): (6)SimDist(X)=1∏k=1,k≠jNPairDist(X,Akm) where PairDist(X) is defined as follow:

Finally, the similarity score of X, Sim(X), is the product of SimSymb(X) and SimDist(X). Again, the transition probability of X is the normalization of Sim(X) in [0,1].

The goal of the refinement phase is to increase the quality of the discovered alignment. An alignment of residues is iteratively removed from the current alignment of substructures and replaced with a new one. The number of iterations is bounded by a user-defined parameter called IterRefine. According to our experimental results (Section 3.2), a good accuracy can be achieved with relatively small values of such parameter (e.g., 10).

The replaced alignment is chosen according to a Badness function defined below.

Let SA = {SA₁, SA₂, …, SA_N} be the final alignment of size w, where each SA_i = {R_i,1, R_i,2,…, R_i,w} is a set of residues. We can view the alignment SA as a matrix R[N,w], where each column represents an alignment of residues and R[i,j] is the j-th aligned residue of the i-th substructure. Our final goal is to compute a Badness score for each column of SA and remove the column that maximizes the Badness score function from SA.

First, given two aligned residues R[i,k] and R[j,k], we define the function PairDistAligned as follows:

The Badness of a generic column k is:

Once the column with the highest Badness score is removed, a new single extension step is performed through the Gibbs sampling procedure described in Section 2.2.

The alignments produced by PROPOSAL are sorted according to the average RMSD across all possible pairs of structures. This sorted list is finally post-processed to filter highly overlapping alignments. Let SAⁱ = {SAⁱ₁, SAⁱ₂,…, SAⁱ_N} be the local alignment of rank i in the sorted list. We define Perc(SAⁱ_k) as the percentage of residues in the substructure SAⁱ_k observed in the previous i − 1 alignments, and Perc(SAⁱ) as the average value of Perc(SAⁱ_k) across all the aligned substructures. If Perc(SAⁱ) is above a given threshold Overlap, the alignment is discarded.

The dataset is divided into four categories, depending on similarity degree and sequence length. Table 1 synthesizes the features of each family with respect to the number of proteins, the average sequence length and the average similarity.

To evaluate the reliability of PROPOSAL we considered different values of w, depending on proteins sequence similarity. We chose w = 10 for the CheY-related proteins' family, w = 12 for the Ferritin family, w = 15 for the Plastocyanin proteins, and w = 20 for the TIM Barrel family. In all experiments, we set AvgOverlap = 50% to reduce the final set of alignments. Table 2 gives the running time of PROPOSAL and the RMSD of the best alignments.

The best alignments have been generated through the 2D alignment of their contact maps. A protein contact map is a 2D matrix storing the distances between all possible amino acids pairs of a 3D protein structure. It is represented as a graph where nodes are amino acids and edges connect nodes having a distance less than a fixed cut-off, usually 7–12 Å. A contact map is a signature of a protein structure with respect to its 3D coordinates (Vassura et al., 2008).

Figures 2–5 show the 10 Å cut-off contact map alignments. It can be seen that a good structural correspondence between proteins is guaranteed even when the value of w increases. In most cases the absence of few edges or the presence of new links between nodes are due to pairs of residues whose distance is very close to the cut-off.

We analyzed label similarity of the four best alignments, by building the sequence logos (Crooks et al., 2004) of mapped residues (Figures 6–9). Each position contains a graphical representation of the frequencies of residues in that position within the final mapping. Amino acids are represented with different colors, depending on their chemical properties: basic residues (K, R, H) are colored in blue, the acidic ones (D, E) in purple, the neutral ones (Q, N, P, S, C) in green, the hydrophobic ones (V, L, I, W, F, M, Y) in orange, and the remaining ones (G, T, A) in red.

Sequence logos reflect the average sequence similarity of proteins within each family: Plastocyanin and TIM Barrel proteins show the best label correspondence. The alignment of Ferritin proteins is quite interesting, since the structural similarity is high, the average LRMSD is very low (0.428 Å, Table 2), but the corresponding sequence logo shows remarkable dissimilarities between mapped residues. This is an example confirming that protein structural similarity and protein sequence similarity are not always related.

Next, we investigated the effects of varying PROPOSAL parameters. The default values are N = 6, w = 15, α = 0.05, and IterRefine = 10. First, we analyse how parameters influence the running time (Figure 10) by varying one parameter and leaving the rest unchanged. Figure 10A depicts the running time varying the number N of structures. Figure 10B deals with the effect of varying w from 1 to 20. Figure 10C reports the PROPOSAL behavior with α ranging from 0.01 to 0.30. Finally, in Figure 10D different values of IterRefine (from 1 to 30) are considered. As expected, when N and w grow and alpha decreases, the running time goes up. Such a trend is even more evident in the TIM Barrel family which has the highest average protein sequence length and similarity.

Figure 11 shows the influence of α and IterRefine on the global accuracy of PROPOSAL. We measured the average RMSD over all the computed alignments. In Figure 11A alpha varies from 0.01 to 0.30 and IterRefine is set to 10, while in Figure 11B iterRefine varies from 1 to 30 and α is set to 0.05. Default values (w = 15 and N = 6) were assigned. As expected, the best performance of our method are obtained with low values of α and high values of IterRefine. However, if we also consider the influence of such parameters on running time (in particular the IterRefine parameter), the best trade-off between speed and accuracy can be achieved with 0.01 ≤ α ≤ 0.1 and IterRefine = 10.

As far as we are concerned, PROPOSAL is the first algorithm proposed for multiple local alignments of protein structures. On the other hand a few existing tools can solve the pairwise local structure alignment problem (Holm and Park, 2000; Xie and Bourne, 2008; Konc and Janezic, 2010). According to the experiment results reported in Konc and Janezic (2010) and Moll et al. (2011), ProBiS and SMAP seem to be the best existing pairwise local structure alignment methods.

In order to compare PROPOSAL with ProBiS and SMAP, we run all the algorithms on a properly defined dataset of pairwise alignments.

First of all, we collected a set of 346 non-redundant literature derived small query motifs (having 4-6 residues), taken from CSA (Catalytic Site Atlas) (Furnham et al., 2013). CSA is a database of hand-annotated entries, containing enzyme active sites (i.e., a set of residues thought to be directly involved in the reaction catalyzed by an enzyme). The complete list of these motifs and the corresponding PDB structures is available in the supplementary material Table S1.

Then, we used LabelHash, which is the state-of-the-art tool for substructure matching, to search for a match between each query motif and the rest of the dataset. Finally we selected all matches with RMSD ≤ 1.5 Å. This resulted in a final reference dataset of 6380 pairwise alignments (the dataset is available in the supplementary material Table S2).

The dataset has many highly dissimilar pairs of proteins. In order to analyse the sequence similarity between the 6380 couples of proteins with the lowest RMSD alignments, we run BLAST and considered the percentage of residues with positive matches in the shortest sequence. We call PPos the latter measure. Among the 6380 couples, 3835 (≃ 60%) have PPos < 5% and 6173 (≃ 97%) have PPos < 15%.

For each couple, we run PROPOSAL with no overlapping filter (AvgOverlap = 100%) and w equals to the number of residues of the query motif. We ran SMAP and ProBiS with default parameter values.

We analyzed the performance of the three methods on the 6380 pairwise alignments, by taking into account three parameters:

We analyzed the average values of these parameters by considering different ranges of PPos similarities. All results are plotted in Figure 12.

PROPOSAL exhibits the highest QMC for highly dissimilar proteins, while for medium and high PPos similarities ProBiS is the best method (Figure 12A). However, in all the tested instances PROPOSAL yields the lowest average RMSD with respect to both ProBiS and SMAP. Furthermore, the difference between RMSDs tends to increase as long as PPos decreases (Figure 12B). We also notice that the average QMC and RMSD of PROPOSAL alignments are approximately constant for all values of PPos, while ProBiS and SMAP seem to be quite sensitive to protein similarity.

Finally, ProBiS is by far the fastest algorithm for all possible ranges of PPos similarity values (Figure 12C), while PROPOSAL and SMAP have similar running times (except for 80% ≤ PPos ≤ 100%, where SMAP is faster). It is worth noting that our method has been designed for solving the multiple alignment problem, while ProBiS and SMAP have been efficiently implemented for comparing pairs of protein structures. Moreover, PROPOSAL and SMAP have been implemented in Java, while ProBiS has been written in C++. Interestingly, our method is faster when PPos ranges from 10 to 30%. However, when proteins are very dissimilar, the convergence of Gibbs sampling in the bootstrap phase may be slower. On the other hand, when proteins are very similar PROPOSAL performs more extension and refinement phases, producing more feasible alignments. A similar trend holds for ProBiS, where the best performance is obtained when PPos ranges from 15 to 60%.

In the last case study, we run PROPOSAL on a different set of 172 motifs taken from CSA to test the capability of our method to detect known conserved binding sites in the multiple case (see supplementary material Table S3).

The dataset has been built by selecting literature derived motifs of proteins belonging to fully qualified EC classes with at most 25 elements. This resulted in a final set of 172 motifs, spanning 162 distinct EC classes.

EC class (Webb, 1993) is a code having the format “EC” followed by four numbers separated by periods. It denotes the type of reaction catalyzed by an enzyme. An EC class is fully qualified if all four numbers are specified (e.g., 1.1.1.149 is fully qualified, while 1.1.1 or 1.1 are not).

For each EC family, we run PROPOSAL on the set of protein structures belonging to that family. We fixed w equals to the number of residues in the corresponding motif and AvgOverlap = 100% (i.e., no overlapping filter). The remaining parameters were set up to the default values.

We filtered out all alignments with average RMSD above 1 Å, taking for each query motif the local alignment with maximum QMC. In case of ties on QMC, the alignment with minimum average RMSD was chosen. PROPOSAL successfully completed all the alignments in about 29 h, with an average QMC of 50.08% and average running time of 10 min. In Table 3 we report motifs with highest QMC and the RMSD of the corresponding alignment (see supplementary material Table S3 for the complete list of results). Results clearly show the ability to identify known motifs from scratch. Out of 172 motifs, 24 have QMC ≥ 75% and 126 have specificity ≥ 50%.

In Torrance et al. (2005); Moll et al. (2010), authors observed that the EC-class coverage of a motif has not been considered for the design of CSA. Consequently, some motifs may be not conserved across all proteins in an EC class. This may be the origin of failures of PROPOSAL on the alignment tasks with QMC <50%). In some cases CSA motifs could contain one or more residues with few global matches. Moreover, two motifs could match mutually exclusive sets of proteins within the corresponding EC class. These cases may cause a drastic increase of average RMSD for that specific motif. Examples of such CSA motifs are reported in Moll et al. (2010). In order to overcome these problems, methods like Geometric Sieving (Chen et al., 2007) can be applied to refine a given motif and increase sensitivity while keeping high specificity values.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.