Searching Protein 3-D Structures in Linear Time
Tetsuo Shibuya
Human Genome Center, Institute of Medical Science, University of Tokyo
4-6-1 Shirokanedai, Minato-ku, Tokyo 108-8639, Japan
tshibuya@hgc.jp

Abstract. Finding similar structures from 3-D structure databases ofproteins is becoming more and more important issue in the post-genomic
molecular biology. To compare 3-D structures of two molecules, biologists
mostly use the RMSD (root mean square deviation) as the similarity
measure. We propose new theoretically and practically fast algorithms
for the fundamental problem of ﬁnding all the substructures of structuresin a structure database of chain molecules (such as proteins), whose
RMSDs to the query are within a given constant threshold. We ﬁrst
propose a breakthrough linear-expected-time algorithm for the problem,
while the previous best-known time complexity was O(N log m), where
N is the database size and m is the query size. For the expected time
analysis, we propose to use the random-walk model (or the ideal chain
model) as the model of average protein structures. We furthermore propose a series of preprocessing algorithms that enable faster queries. We
checked the performance of our linear-expected-time algorithm throughcomputational experiments over the whole PDB database. According to
the experiments, our algorithm is 3.6 to 28 times faster than previously
known algorithms for ordinary queries. Moreover, the experimental results support the validity of our theoretical analyses.

1

Introduction

3-D structure database searching of molecules, especially proteins, plays a very
important role in molecular biology [2,8,10]. For example, if we have proteins
whose structures are known but their functions are unknown, we may be able
to predict their functions by searching for similar structures whose functions are
known, as structurally similar proteins tend to have similar functions. Moreover,more and more protein structures are solved today with the aid of state-of-theart technologies such as nuclear magnetic resonance (NMR) techniques, as seenin the rapid growth of the PDB database [3]. Thus, faster searching techniques
are seriously needed for the molecular structure databases.
A protein is a chain of amino acids. Thus, its structure can be represented bya sequence of 3-D coordinates, each of which corresponds to the position of a
speciﬁed atom (the Cα atom is usually used) of each amino acid. Such molecules
are called chain molecules. There are also many other important chain moleculesin living cells, such as DNAs, RNAs and glycans. The RMSD (root mean square
S. Batzoglou (Ed.): RECOMB 2009, LNCS 5541, pp. 1–15, 2009.
c Springer-Verlag Berlin Heidelberg 20092

T. Shibuya

deviation) [1,7,12,13,16,19] is the fundamental measure to determine the geometric similarity between two same-length sequences of 3-D coordinates. It has been
widely used not only for molecular structure comparison, but also for various
problems in various ﬁelds, such as computer vision and robotics. It is deﬁned as
the square root of the minimum value of the average squared distance between
each pair of corresponding atoms, over all the possible rotations and translations (see section 2.2). In this paper, we consider one of the most fundamental
RMSD-related problems as follows.
Problem. Given a structure database D of chain molecules and a query structure Q, ﬁnd all the substructures of the structures in D whose RMSDs to Q are
at most a given ﬁxed threshold c, without considering any insertions or deletions.
In general, c should be set to a ﬁxed constant proportional to the distance
between two adjacent atoms of the chain molecules. In the case of proteins,the distance between two adjacent Cα atoms is around 3.8Å, while two protein
structures are said to be similar to each other if their RMSD is smaller than 1
or 2Å.
Our results. The best-known worst-case/expected time complexity of the problem was O(N log m) [16,19], where N is the database size (i.e., the sum of thelengths of all the structures in the database) and m is the query size. We propose the ﬁrst linear-expected-time (i.e., O(N )) algorithm. To analyze the expected time of the algorithm, we give an assumption that the structures in the
database follow a model called the random-walk model (see section 2.4). We also
propose several preprocessing algorithms that enable faster queries. We ﬁrst propose an O(N√ log N )-time and O(N )-space preprocessing algorithm that enables
O(m+N/ m)-expected-time query, for queries of a ﬁxed length. We next extend
it to an O(N log2 N )-time and√ O(N log N )-space preprocessing algorithm thatenables the same O(m + N/ m)-expected-time query, for queries of arbitrary
lengths. We also propose an O(N log N )-time and O(N )-space preprocessing algorithm that enables O( √Nm + m log(N/m))-expected-time query, for queries ofarbitrary lengths.
We also examine the performance of our algorithms by computational experiments on the whole PDB database. Our linear-time algorithm is much fasterthan any of previous algorithms, i.e., 3.6 to 28 times faster to search for substructures whose RMSDs are at most 1Å. Moreover, no inconsistency is observed
between the theoretical results and the experimental results, which means ourrandom-walk assumption is very reasonable for analyses of algorithms for protein
structure database search.
The organization of this paper. In section 2, we describe the basic deﬁnitions√
and related previous work as preliminaries. In section 3, we propose an O(N m)algorithm for the problem above, where N is the database size and m is thelength of the query. We next improve it to obtain the linear-time algorithm
in section 4. In section 5, we further extend our algorithms for faster queries
after preprocessing. In section 6, we examine the performance of our algorithms
against the PDB database. In section 7, we conclude our results.

Searching Protein 3-D Structures in Linear Time

2

Preliminaries

2.1

Notations and Deﬁnitions

3

A chain molecule is represented like S = {s1 , s2 , . . . , sn }, where si denotesthe 3-D coordinates of the i-th atom. The length n of S is denoted by |S|.
A structure S[i..j] = {si , si+1 , . . . , sj } (1 ≤ i ≤ j ≤ n) is called a substructure of S. R · S denotes the structure S rotated by the rotation matrix R, i.e.,
R · S = {Rs1 , Rs2 , . . . , Rsn }. v t denotes the transpose of the vector v and AT
denotes the transpose of the matrix A. trace(A) denotes the trace of the matrix
A. |v| denotes the norm of the vector v. 0 denotes the zero vector. x denotes
the expected value of x. var(x) denotes the variance of x. P rob(X) denotes theprobability of X.
In the rest of this paper, we consider that the target database D consistsof one long structure P = {p1 , p2 , . . . , pN }, and we let Q = {q1 , q 2 , . . . , q m }denote the query structure, where m is supposed to be smaller than N . Our
problem is to ﬁnd all the positions i such that the RMSD (see section 2.2 for its
deﬁnition) between P[i..i + m − 1] and Q is at most a given ﬁxed threshold c. An
ordinary database may contain more than one structure, but the problem against
such databases can be reduced to the problem against databases with only one
structure, by concatenating all the database structures into one structure andignoring substructures that cross over the boundaries of concatenated structures.
2.2

RMSD: The Root Mean Square Deviation

The RMSD (root mean square deviation) [1,7,12,13,16,19] between two 3-D co. , sn } and T = {t1 , t2 , . . . , tn } is deﬁned as theordinate sequences S = {s1 , s2 , . .nminimum value of ER,v (S, T) = n1 i=1 |si − (R · ti + v)|2 over all the possible rotation matrices R and translation vectors v. Let RM SD(S, T) denote theminimum value, and let R̂(S, T) and v̂(S, T) denote the rotation matrix andthe translation vector that minimizes ER,v (S, T). RM SD(S, T), R̂(S, T) andv̂(S, T) can be computed in O(n) time as follows.
If the rotation matrix R is ﬁxed, ER,v (S, T) is known to be minimized when
the centroid (center of mass) of R · T is translated to the centroid of S by thetranslation vector v, regardless of what the rotation matrix R is. It means thatv̂(S, T) can be computed in linear time if we are given R̂(S, T). Moreover, it alsomeans that the problem of computing the RMSD can be
reduced to a problem
of ﬁnding R (i.e., R̂(S, T)) that minimizes ER(S, T) = ni=1 |si − R · ti |2 , bytranslating both S and T so that both of their centroids are moved to the originof the coordinates, which can be done in linear time.
After translating both structures so that both of their centroids are moved
1
Let
to thenorigin, wet can compute R̂(S, T) in linear time as follows [1,12,13].

J = i=1 si · ti . Clearly, J can be computed in O(n) time. Then ER (S, T) can
1Here, S and T are the translated structures which are diﬀerent from the original S
and T in the original problem (that optimizes both rotation and translation).

4

T. Shibuya
be described as ni=1 (si t si + ti t ti )−2·trace(R·J), and trace(R·J) is maximized
when R = V U T , where U ΛV is the singular value decomposition (SVD) of J.
Thus R̂(S, T) can be obtained from J in constant time, as J is a 3 × 3 matrix
and the SVD can be computed in O(d3 ) time for a d × d matrix [11]. Note that
there are degenerate cases where det(V U T ) = −1, which means that V U T is areﬂection matrix. See [1,7,12] for the details of the degenerate cases. Finally, we
can compute the RMSD in linear time once we have obtained R̂(S, T). In total,
we can compute the RMSD in O(n) time.
2.3

Previous Best-Known Searching Algorithms

According to the previous section, we can compute the RMSD between any
substructure P[i..i + m − 1] and the query Q in O(m) time. Consequently, we
can solve our problem in O(N m) time by checking the RMSDs between thequery and all the O(N ) substructures of length m in the database.
Schwartz and Sharir [16] proposed a more sophisticated approach for the problem that solves it in O(N log N ) time, based on the convolution technique using
the FFT (fast Fourier transform) [5]. Shibuya [19] also proposed a diﬀerent algorithm with the same time complexity, also based on the convolution technique.
These algorithms are not faster than the naive algorithm when N  m. But
this time bound can be easily improved to O(N log m) as follows. Break P intoO(N/m) substructures of length m + O(m) each of which overlaps with its adjacent fragment with overlap length m − 1. Then our problem can be solved inO(N log m) time by applying the above O(N log N )-time algorithm against eachfragment. The expected time complexity of these algorithms are all the same astheir worst-case time complexity, and no algorithm with better expected timecomplexity is known. But note that the above FFT-based O(N log m)-time algorithm is not practically faster than the naive O(N m)-time algorithm in casem is not large enough, and it is rarely used in practice.
For the problem, a linear-size indexing data structure called the geometricsuﬃx tree [17,18] is known to enable practically faster query than the above
algorithms. But its worst-case query time complexity is still O(N m), while weneed O(N 2 ) time to construct the data structure. In fact, there have been no
known indexing algorithms whose theoretical query time complexity is smaller
than the above O(N log m) bound.
2.4

The Random-Walk Model for Chain Molecule Structures

The random-walk model for chain molecule structures is a simple but usefulmodel for analyzing their behavior [4,6,9,15]. The model is also called the freelyjointed chain model or the ideal chain model. In the model, we assume that thestructure of a chain molecule is constructed as a result of a random walk in 3-D
space. It is useful in various analyses in molecular physics, as it reﬂects propertiesof structures of real chain molecules very well [4].
Consider a chain molecule S = {s0 , s2 , . . . , sn } of length n + 1, in which the
distance between two adjacent atoms is ﬁxed to some constant . Note that

Searching Protein 3-D Structures in Linear Time5the length between two adjacent Cα atoms in a protein structure is constantly
3.8Å, as mentioned in section 1. In the random-walk model, a bond between
two adjacent atoms, i.e., bi = si+1 − si , is considered as a random vector thatsatisﬁes |bi | = , and bi is independent from bj for any i and j (i = j). If n
is large enough, the distribution of the end-to-end vector sn − s0 is known toconverge to the Gaussian distribution in 3-D space, in which sn − s0  = 0
and |sn − s0 |2  = n · 2 . In the distribution, the probability (or probability
density) that sn −s0 is located at some position (x, y, z) is Wn, (x, y, z)dxdydz =
2
2
2
2
3
3
2 −3(x +y +z )/2n dxdydz.
( 2πn2) e
In our theoretical analysis in later sections, we will give an assumption thatthe structures in the target database follow the random-walk model. Our experimental results on the PDB database in section 6 show high consistency with
our theoretical analyses based on the model.

3

An O(N√m) Algorithm

In this section, we propose an algorithm that will be the basis for our algo√rithms in later sections. The algorithm given in section 3.2 achieves O(N m)time under the random-walk assumption, by ﬁltering out most of the dissimilar
substructures (i.e., substructures with RMSDs larger than the given thresholdc) before computing the actual RMSDs, based on a measure described in section 3.1. Its time complexity will be analyzed in section 3.3.
3.1

An Eﬃciently-Computable Lower Bound for the RMSD

We ﬁrst propose a nontrivial, but easily-computable lower bound for the
RMSD between any two structures with the same length. Let Ulef t denote {u1 , u2 , . . . , uk/2 } and Uright denote {uk/2+1 , uk/2+2 , . . . , u2·k/2 }for a structure U = {u1 , u2 , . . . , uk }. Let G(U) denote the centroid (cenkter of mass) of the structure U, i.e., G(U) = k1 i=1 ui . Let F (U) denote|G(Ulef t ) − G(Uright )|/2, which means the half of the distance between the centroids of Ulef t and Uright . For any two structures S and T with the same lengthn, we deﬁne D(S, T) as |F (S)− F (T)| if n is an even integer. If n is an odd

integer, we deﬁne D(S, T) as n−1
n |F (S) − F (T)|. From now on, we prove that
D(S, T) is a lower bound of the RMSD between S and T.
Let si = si − G(S), and ti = ti − G(T). In case n is an even integer, we provethat D(S, T) is always smaller than or equal to RM SD(S, T), as follows:

 n
n
1 
1 
|si − R̂(S, T) · ti |2 ≥|s − R̂(S, T) · ti |RM SD(S, T) = 
n i=1
n i=1 i1  
1|{s − R̂(S, T) · ti }| + |n i=1 i
n
n/2

≥

n


{si − R̂(S, T) · ti }|

i=n/2+1

6

T. Shibuya= |G(Slef t ) − G(S) − R̂(S, T) · {G(Tlef t ) − G(T)}|/2+|G(Sright ) − G(S) − R̂(S, T) · {G(Tright ) − G(T)}|/2= |G(Slef t ) − G(Sright ) − R̂(S, T) · {G(Tlef t ) − G(Tright )}|/2≥ |{|G(Slef t ) − G(Sright )| − |G(Tlef t ) − G(Tright )|}|/2= D(S, T).

(1)

Note that we used the fact that G(S) = {G(Slef t ) + G(Sright )}/2 and G(T) ={G(Tlef t ) + G(Tright )}/2 in the above.
In case n is an odd integer, we prove the same as follows:n−1
RM SD(S, T) ≥RM SD(S− , T− )nn−1
D(S− , T− ) = D(S, T),≥(2)nwhere S− denotes {s1 , s2 , . . . , sn−1 }, and T− denotes {t1 , t2 , . . . , tn−1 }.
If the above lower bound can be computed very eﬃciently, we may solve
our problem (presented in section 1) more eﬃciently by ﬁltering out hopelessly
dissimilar substructures before computing the actual RMSD value. In fact, we
can compute the above lower bound D(P[i..i + m − 1], Q) for all the positions
i in linear time, as follows. For any m, we can compute G(P[i..i + m/2 − 1])for all the positions i such that 1 ≤ i ≤ N − m/2 + 1 in O(N ) time, as1(pi−1 − pi+m/2−1 ).
G(P[i..i + m/2 − 1] = G(P[i − 1..i + m/2 − 2]) − m/2
It means that F (P[i..i + m − 1], Q) and consequently D(P[i..i + m − 1], Q) can
be computed for all the positions i such that 1 ≤ i ≤ N − m + 1, also in O(N )
time.
3.2

The Algorithm

Our basic algorithm is simple. It uses the above lower bound to ﬁlter out some(hopefully most) of the substructures in the database before the time-consuming
RMSD computation, as follows.
Algorithm 1
1 Compute Di = D(P[i..i + m − 1], Q)
for all i such that 1 ≤ i ≤ N − m + 1.
2 for (all i such that 1 ≤ i ≤ N − m + 1) {3
if (Di ≤ c) {4
if (RM SD(P[i..i + m − 1], Q) ≤ c)5
{ output(“P[i..i + m − 1] is similar to the query Q.”) }6}
7 }Searching Protein 3-D Structures in Linear Time7

The above algorithm is valid, i.e., it enumerates all the positions of the substructures whose RMSDs to the query are at most c, because Di = D(P[i..i+m−1], Q)
is always smaller than or equal to RM SD(P[i..i + m − 1], Q). Let the number oftimes of the RMSD computation in line 4 be N  (≤ N ). Then, the time complexity of the above algorithm is O(N + N  m), as the line 1 of the algorithm requires
only O(N ) time according to the discussion
in section 3.1. In the next section,√
we will prove that N   is in O(N/ m), if the structures in the database follow
the random-walk model.
3.3

Computational Time Analysis

Consider a structure S = {s1 , s2 , . . . , s2n } of length 2n that follows the randomwalk model. In this section, we let the distance between two adjacent atoms( in section 2.4) be 1 without loss of generality, i.e., we consider the distance
betweentwo adjacent atoms as the unit of distance. Then si can be represented
i−1
as s1 + j=1 bj , where bi is an independent random vector that satisﬁes |bi | = 1.
Let H(S) = G(Slef t ) − G(Sright ). Notice that F (S) = |H(S)/2|. Then thefollowing equation holds:
H(S) =

n
i−1
2n
i−1


1 
1
(s1 +
bj ) −
(s1 +
bj )
n i=1
n i=n+1
j=1
j=1

= −{

n

i

n
i=1

· bi +2n−1


2n − i· bi }.[1mm]
n
i=n+1

(3)

Let bi denote ni · bi if i ≤ n and 2n−in · bi if i > n. Then H(S) can be described2n
as i=1 bi . Let zi denote the z coordinate of bi and zi denote the z coordinate ofbi . It is easy to see that zi  = 0 and var(zi ) = 1/3, as bi is a random vector that

2+δ
2n
2n

, where δ issatisﬁes |bi | = 1. Let Mn = i=1 |zi − zi |2+δ /i=1 var(zi )some positive constant. According to Lyapunov’s central limit theorem [14], the2ndistribution of i zi converges to the Gaussian distribution, if Mn converges
to 0 as n grows up to inﬁnity for some δ such that δ > 0. It can be proved asfollows:


2+δ
2+δ
 2n
 2n
2n
2n




|zi |2+δ / |zi |2≤
|zi |2 / |zi |2
Mn =
i=1

={

i=1

2n


|zi |2 }−δ/2 = {

i=1

i=1
n


i=1
2n−1


1
·{
i2 +
(2n − i)2 }}−δ/2
3n2 i=1i=n+1

1 −δ/22}={ n+→ 0 (n → ∞).[2mm]
(4)
9
9n
2n
Hence, we conclude that i zi converges to the Gaussian distribution. It
also means that H(S) converges to the Gaussian distribution in 3-D space if n

8

T. Shibuya

grows up to inﬁnity, as the same discussion can be done for the other two axes
(x and y). The variance of H(S) is:var(H(S)) = |H(S)|2  − |H(S)|2 = |H(S)|2 
= |

n

i
i=1

=

n

· bi +2n−1


2n − i· bi |2 
n
i=n+1

n
2n−1

1
2
1  2
2
≈ n, [1mm]{i+(2n − i)2 } = n +n2 i=1
3
3n
3
i=n+1

(5)

as bi · bj  = 0 if i = j. Moreover, it is easy to see that H(S) = 0. Thus
the distribution of H(S) is the same as the distribution of random walks oflength 2n/3. Hence the probability distribution of H(S) is Zn (x, y, z)dxdydz =
2
2
2
3
9
( 4πn
) 2 e−9(x +y +z )/4n dxdydz. Consequently, the probability (or probability
239) 2 e−9r /4n dr. Integrating
density) that |H(S)| = r is Zn (r)dr = 4πr2 ( 4πnyZn (r)dr, we obtain P rob(x ≤ |H(S)| ≤ y) = r=x Zn (r)dr. Zn (r) takes the√2√
maximum value at rmax = 3 n and Zn (rmax ) = 6e−1 / √πn. Thus P rob(x ≤−1|H(S)| ≤ y) is at most (y − x) · Zn (rmax ) = 6e (y − x)/ πn for any x and y(x < y).
Therefore, for any structure T such that |T| = |S| = 2n, the probability
P rob(|D(S, T)| ≤ c) = P rob(F (T) − c ≤ F (S) ≤ F (T) + c) = P rob(2√ · F (T) −
−1
(r
)
=
24c
·
e
/
πn, which2c ≤ |H(S)|
≤
2
·
F
(T)
+
2c))
is
at
most
4
·
c
·
Z
n max
√
is in O(1/ n) as c is a ﬁxed constant. Notice that there is no assumption on the
structure T in this analysis.
≈ 1, the probability P rob(Di ≤ c) in the line 3 of√the algorithm in section 3.2 is in O(1/ m) if P follows the random-walk model,√
no matter what the query structure Q is. It means that N   is in O(N/ m).
Therefore, we conclude that√ the expected time complexity of the algorithm isO(N + N   · m) = O(N m), under the assumption that the structures inthe database follow the random-walk model.2 Note that the worst-case timecomplexity of the algorithm is still O(N m) as N  can be in O(N ) at worst, but
it should be rare under the random-walk assumption.
m−1
mConsequently, as

4

The Linear-Time Algorithm

We next improve the algorithm 1 to obtain better expected time complexity.
From the deﬁnition of the RMSD, we can deduce that√
RM SD(S, T) ≥ [2mm] t · {(RM SD(Slef t , Tlef t ))2 +≥2√(RM SD(Sright , Tright ))2 }1/2
t · {(D(Slef t , Tlef t ))2 + (D(Sright , Tright ))2 }1/2 ,(6)

The same discussion can be done in case the query structures, instead of the databasestructures, follow the random-walk model.

Searching Protein 3-D Structures in Linear Time9√where t = |Slef t |/|S|√= |Sright |/|S| ≈ 1/2. Let Dlef t (S, T) = t · D(Slef t , Tlef t )and Dright (S, T) = t · D(Sright , Tright ). The expression (6) can also be used
as a lower bound of the RMSD for another valid ﬁltering algorithm, as follows:
Algorithm 2
1 For all i such that 1 ≤ i ≤ N − m + 1, computeDi = {(Dlef t (P[i..i + m − 1], Q))2 +(Dright (P[i..i + m − 1], Q))2 }1/2 .
2 for (all i such that 1 ≤ i ≤ N − m + 1) {3
if (Di ≤ c) {4
if (RM SD(P[i..i + m − 1], Q) ≤ c)5
{ output(“P[i..i + m − 1] is similar to the query Q.”) }6}
7 }The only diﬀerence from the Algorithm 1 is the lower bound Di used in the line
1. Note that the time complexity of the line 1 is still O(N ).
If Di ≤ c in line 3, both Dlef t (P[i..i+ m− 1], Q) and Dright (P[i..i+ m− 1], Q)
must also be at most c. According to the discussion in section 3.3, the two probabilities P rob(Dlef t (P[i..i+ m− 1], Q) ≤√c) and P rob(Dright (P[i..i+ m− 1], Q) ≤c) are independent and both in O(1/ m), under the assumption that P fol
lows the√ random-walk model. Thus, P rob(Di ≤ c) in line 3 must be inO((1/ m)2 ) = O(1/m). Therefore the expected number of RMSD computations in line 4 should be in only O(N/m), and consequently the expected timecomplexity spent in the lines 4–6 of the above algorithm is O(N ). Thus, the total
expected time complexity of the algorithm 2 is O(N ) under the random-walkassumption. Note that the worst-case time complexity is still O(N m), but itshould be very rare under the random-walk assumption.

5

Faster Queries after Preprocessing

5.1

Preprocessing for Queries of a Fixed Length

From now on, we further improve the query time complexity by allowing preprocessing on the database structure P. First, we consider the case where
each query has the same length m. According to section 3.1, we can compute
F (P[i..i + w − 1]) for all i in O(N ) time for a ﬁxed value of w. Let Lw be thesorted list of i according to the value of F (P[i..i + w − 1]), which can be obtained in O(N log N ) time. By doing a binary search on Lw , we can ﬁnd all thei such that x ≤ F (P[i..i + w − 1]) ≤ y in O(log N + occ) time for any x and y,where occ is the number of the outputs. Hence, we can list all the i such that
D(S, P[i..i + w − 1]) ≤ c in O(log N + occ) time for any structure S of length
w by utilizing Lw , where c is some constant and occ is the number of outputs,
as F (S) − c ≤ F (P[i..i + w − 1]) ≤ F (S) + c if D(S, P[i..i + w − 1]) ≤ c. Our
preprocessing algorithm in this section computes F (P[i..i + m − 1]) for all i and10

T. Shibuya

sorts them to obtain Lm , where m = m/3 , which can be done in O(N log N )time in total.
Consider yet another lower bound for the RMSD, as follows:3

3
m 

· { (D(P[i + (j − 1) · m ..i + j · m − 1],
Di =
m
j=1

≤

Q[1 + (j − 1) · m ..j · m ]))2 }1/2
3m 
· { (RM SD(P[i + (j − 1) · m ..i + j · m − 1],m
j=1

Q[1 + (j − 1) · m ..j · m ]))2 }1/2
≤ RM SD(P[i..i + m − 1], Q).

(7)

Notice that D(P[i + (j − 1) · m ..i + j · m − 1], Q[1 + (j − 1) · m ..j · m ]) ≤m
c m
≤ c. Let Xj be the set of all positions i such that for any j, if Dim
D(P[i + (j − 1) · m ..i + j · m − 1], Q[1 + (j − 1) · m ..j · m ]) ≤ c m
 (for
1 ≤ j ≤ 3). By using Lm , we can ﬁnd all i ∈ Xj in O(log N + |Xj |) time forany of j = 1, 2, 3 after we have computed F (Q[1 + (j − 1) · m ..j · m ]). Notethat F (Q[1 + (j − 1) · m ..j · m ]) for all of j = 1, 2, 3 can be computed in O(m)time. Let X be the set of common integers of the three sets X1 , X2 , and X3 .
X can be obtained in expected O(|X1 | + |X2 | + |X3 |) time by using a hashingtechnique. Notice that the positions i such that RM SD(P[i..i + m − 1], Q) ≤ cmust be included in X. For substructures at the positions i ∈ X, we ﬁnallyhave to compute the RMSD to check whether the actual RMSD is at most c, if
Di ≤ c. It can be done
√ in at most O(m · |X|) time.
|Xj | is in O(N/ m) under the assumption that P follows the random-walkaremodel. Moreover, as the structures P[i+(j−1)·m ..i+j·m −1]√ with diﬀerent j1.5
independent random walks, |X| is estimated as O(N/( m)3 ) = O(N/m√ ).

Thus the total expected query timecomplexity
utilizing
L
is
O(m
+
N/
m+
m
√
m · N/m1.5 + log N ) = O(m + N/ m).
5.2

Preprocessing for Queries of Arbitrary Lengths

We next consider queries of arbitrary lengths. For such queries, consider computing Lw for all w such that w is a power of 2, i.e., representable as 2d for some
integer d. They can be obtained in O(N log2 N ) time, as the number of diﬀerent
w is in O(log N ). Let m be the largest power of 2 such that 3m ≤ m. Thenthe inequality (7) also holds for this case. The only diﬀerence is that m is somepower of 2 that satisﬁes m/6 < m ≤ m/3, while m = m/3 in the previous
section. Thus, according to the same discussion as in the√ previous section, we
obtain the same query time complexity, i.e., O(m + N/ m). A problem is that3

We can consider another linear-expected time algorithm based on the Algorithms 1
and 2 by replacing Di or Di with Di . We will examine the performance of the algorithm (denoted as A3 in section 6) through computational experiments in section 6.

Searching Protein 3-D Structures in Linear Time11

the algorithm requires O(N log N ) space to store all the Lw , which might beundesired for huge databases.
5.3

Preprocessing with Linear Space for Queries of Arbitrary
Lengths

In this section, we propose another preprocessing algorithm that uses only O(N )space for queries of arbitrary lengths. Consider dividing P into substructures oflength 2d for each d such that 1 ≤ d ≤ log2 N . By doing so, we get substructures
Pk,d = P[(k − 1) · 2d + 1..k · 2d ] (1 ≤ k ≤ N/2d ) for each d. There are only O(N )number of substructures denoted as Pk,d , even if we enumerate all the possible
k and d. G(Pk,1 ) (see section 3 for its deﬁnition) can be computed in constanttime for each k. Moreover, G(Pk,d ) = {G(P2k−1,d−1 ) + G(P2k,d−1 )}/2. Thus,
we can compute G(Pk,d ) for all the possible k and d such that 1 ≤ k ≤ N/2d
and 1 ≤ d ≤ log2 N in O(N ) time by dynamic programming. Consequently, allof the F (Pk,d ) values can also be computed in O(N ) time. For each d (1 ≤ d ≤log2 N ), let Kd be the sorted list of integers k (1 ≤ k ≤ N/2d) according to the
F (Pk,d ) values. Kd can be computed in O((N log N )/2d ) time. Our preprocessing
algorithm in this section computes all these F (Pk,d ) and Kd for all the k and d,which can be done in O(N log N ) time in total.
By doing a binary search on Kd , we can ﬁnd all the k such that x ≤ F (Pk,d ) ≤y for any x, y, and d, in O(log(N/2d) + occ) time, where occ is the number ofthe outputs. Hence, if we are given any structure S of length 2d and the value ofF (S), we can list all the k such that D(S, Pk,d ) ≤ c in O(log(N/2d ) + occ) time,
as F (S) − c ≤ F (Pk,d ) ≤ F (S) + c iﬀ D(S, Pk,d ) ≤ c. Let dQ be the largest
d such that, for any i, there exists some k such that Pk,d , Pk+1,d and Pk+2,dare substructures of P[i..i + m − 1]. Explicitly, dQ = log2 (m + 1) − 2. Let
wQ = |Pk,dQ | = 2dQ . Notice that wQ > m/8. Let IpQ be a set of integers whose
remainder is p− 1 when divided by wQ (1 ≤ p ≤ wQ ), i.e., integers representableas p + j · wQ + 1 with some integer j.
Now consider comparing the query Q and substructures P[i..i + m − 1]such that i ∈ IpQ . There are O(N/wQ ) = O(N/m) such substructures. Let
kQ = (i − 1)/wQ  + 1. Then, PkQ ,dQ , PkQ +1,dQ , and PkQ +2,dQ are substructures of P[i..i + m − 1]. Let PQ,i,j = PkQ +j−1,dQ for j = 1, 2, and 3. Let
Qp,j = Q[p + (j − 1) · wQ ..p + j · wQ − 1] for j = 1, 2, and 3, and let
31Di∗ = 2√
{ j=1 (D(PQ,i,j , Qp,j ))2 }1/2 . Then, Di∗ is also a lower bound of2
RM SD(P[i..i + m − 1], Q), as shown in the following inequality:
1 (8)Di∗ ≤ √ { (D(PQ,i,j , Qp,j ))2 }1/2 < RM SD(P[i..i + m − 1], Q).
2 2 j=1√
Notice that D(PQ,i,j , Qp,j ) < 2 2c for any j, if Di∗ < √c. Given the value ofF (Qp,j ), we can list all i ∈ IpQ such that D(PQ,i,j , Qp,j ) < 2 2c in O(log(N/m)+occ) time for any j, by a binary search on the list KdQ , where occ is the number ofthe i to be listed. Let the list be Yj (j = 1, 2, 3). Note that F (Qp,j ) for all j and p(1 ≤ j ≤ 3, 1 ≤ p ≤ wQ ) can be computed in O(m) time in total.
3

12

T. Shibuya

Then the same discussion as in section 5.1 can√be done. According to thediscussions in section 3.3, |Yj | is in O((N/m)/ m) = O(N/m1.5 ), under
the random-walk assumption on P. The set of the start positions i ∈ IpQ
of similar (i.e., the corresponding RMSD is at most c) substructures must
be included in all of the three lists: Y1 , Y2 and Y3 . Thus, the next thing
to do is to choose the common positions from the three lists. By using ahashing technique, it can be done in O(|Y1 | + |Y2 | + |Y3 |) time, which isO(N/m1.5 ) under the random-walk assumption. Let Y denote the list of thepositions commonly listed in Y1 , Y2 and Y3 . As PQ,i,1 , PQ,i,2 and P√Q,i,3 are
independent random walks, |Y | is estimated to be in O((N/m)/( m)3 ) =O(N/m2.5 ). For substructures at the positions i ∈ Y , we ﬁnally have tocompute the RMSD to check whether the actual RMSD is at most c, if
Di∗ < c. It takes at most O(m · |Y |) = O(N/m1.5 ) expected time under
the random-walk assumption. Thus the total computational time to enumerate all the positions i of similar substructures such that i ∈ IpQ is O(N/m1.5 +log(N/m)).
To enumerate all the positions of similar structures, we execute the above
for all p (1 ≤ p ≤ wQ ). Thus the total√ expected query time complexity isO(m+wQ ·(N/m1.5 +log(N/m))) = O(N/ m+m log(N/m)) under the randomwalk assumption.

6

Computational Experiments on the PDB Database

We did computational experiments using the whole PDB database [3] of thedate September 5th, 2008. Note that more detailed results will be given in thefull version of this paper. The database contains 52,821 entries, which include
244,719 chains of proteins. The total number of amino acids of all the chainsis 38,267,694. We used the Cα coordinates as the representative coordinatesof each amino acid. We did experiments of searching queries of 10 diﬀerent
lengths. In each experiment, we selected 100 substructures of each speciﬁed
length randomly from the whole database, as sample queries. We used 1 CPUof 1200MHz UltraSPARC III Cu on a SunFire 15K super computer for eachexperiment.
Table 1 shows the results. Each column shows the result of the queries ofeach speciﬁed length. The ‘#Substructures’ row shows the number of substructures of each speciﬁed length in the PDB database. The ‘#Hits’ row
shows the average number of hits, i.e., the average number of substructures
whose RMSDs to queries are at most 1Å, among the results of the 100
random queries of each speciﬁed length. The ‘A1’, ‘A2’, ‘A3’, ‘Naive’, and‘FFT’ rows
√ show the average computation time for the 5 algorithms: A1:
the O(N m)-time algorithm proposed in section 3, A2: the linear-time algorithm proposed in section 4, A3: another linear-time algorithm that uses thelower bound Di proposed in section 5.1, Naive: the previously known standard O(N m)-time algorithm, and FFT: the previously known O(N log N )-timealgorithm based on the FFT (introduced in section 2.3). The algorithms A1,

Searching Protein 3-D Structures in Linear Time

13

Table 1. Experimental results over the PDB database for various-length queries
Query length
20
40
60
80
100
#Substructures 33,722,208 29,299,006 25,273,633 21,692,707 18,634,620#Hits
15,093.64
38.07
27.36
32.90
28.61
A1 (sec)
119.87
98.94
98.94
92.43
85.80
A2 (sec)
117.39
58.86
44.01
36.41
36.25
A3 (sec)
151.52
74.56
33.63
25.54
20.46
Naive (sec)423.52447.01450.39442.13428.06
FFT (sec)
551.94531.92505.52463.01425.77
Query length
120
140
160
180
200
#Substructures 16,134,096 14,084,515 12,362,509 10,884,548 9,559,056
#Hits
27.26
27.71
16.01
17.70
23.21
A1 (sec)
75.58
71.22
59.47
63.48
59.98
A2 (sec)
32.84
30.39
27.27
23.12
25.71
A3 (sec)
17.30
15.85
14.25
12.78
12.91
Naive (sec)415.24395.54378.87361.43342.50FFT (sec)
399.83367.76330.57307.89293.03A2 and A3 are all the same algorithms, except for the lower bounds used inthem.
In the experiments, we achieved 3.6 to 28 times speed-up against any of theprevious algorithms for any-length queries, if we choose to use the lower boundDi when the query is longer than 40, and choose Di otherwise. Moreover, theexperiments show that our linear-expected-time algorithms run actually in linear
time on the PDB database, i.e., the algorithm is not inﬂuenced by the diﬀerenceof query lengths.4 It means our random-walk assumption is very reasonable foranalyses of protein structure databases.

7

Concluding Remarks

We proposed the ﬁrst linear-expected-time algorithm for searching similar substructures from structure databases, based on the RMSD measure. Moreover, we
proposed several preprocessing algorithms that enable theoretically even faster
queries. The performance of our algorithms is examined by computational experiments on the whole PDB database.
As for the future work, it would be very interesting to apply our techniques against protein alignment problems that consider insertions and deletions,though it is known to be theoretically much more diﬃcult. Another challenging
task would be to design a deterministically linear-time algorithm for our problem. It is also very interesting to extend our techniques against similar problems
4

Notice that the numbers of substructures of some speciﬁed length in the databasedecreases as the length increases.

14

T. Shibuya

in higher dimensions. Moreover, our technique should be applicable to many
problems in other research ﬁelds such as robotics and computer vision.

Acknowledgement
The author would like to thank Jesper Jansson, Gregory Kucherov, and Kunihiko Sadakane for invaluable discussions. This work was partially supportedby the Grant-in-Aid for Young Scientist (B) No. 20700264 from the Ministry of
Education, Culture, Sports, Science and Technology of Japan. The author used
the super computer system of the Human Genome Center, Institute of Medical
Science, University of Tokyo.

References
1. Arun, K.S., Huang, T.S., Blostein, S.D.: Least-squares ﬁtting of two 3-D point sets.
IEEE Trans. Pattern Anal. Machine Intell. 9, 698–700 (1987)
2. Aung, Z., Tan, K.-L.: Rapid retrieval of protein structures from databases. Drug
Discovery Today 12, 732–739 (2007)
3. Berman, H.M., Westbrook, J., Feng, Z., Gilliland, G., Bhat, T.N., Weissig, H.,
Shindyalov, I.N., Bourne, P.E.: The protein data bank. Nucl. Acids Res. 28, 235–
242 (2000)
4. Boyd, R.H., Phillips, P.J.: The Science of Polymer Molecules: An Introduction Concerning the Synthesis. In: Structure and Properties of the Individual Molecules
That Constitute Polymeric Materials. Cambridge University Press, Cambridge(1996)
5. Cooley, J.W., Tukey, J.W.: An algorithm for the machine calculation of complex
Fourier series. Math. Comput. 19, 297–301 (1965)
6. de Gennes, P.-G.: Scaling Concepts in Polymer Physics. Cornell University Press(1979)
7. Eggert, D.W., Lorusso, A., Fisher, R.B.: Estimating 3-D rigid body transformations: a comparison of four major algorithms. Machine Vision and Applications 9,272–290 (1997)
8. Eidhammer, I., Jonassen, I., Taylor, W.R.: Structure comparison and structure
patterns. J. Comput. Biol. 7(5), 685–716 (2000)
9. Flory, P.J.: Statistical Mechanics of Chain Molecules. Interscience, New York (1969)
10. Gerstein, M.: Integrative database analysis in structural genomics. Nat. Struct.
Biol., 960–963 (2000)
11. Golub, G.H., Van Loan, C.F.: Matrix Computation, 3rd edn. John Hopkins University Press (1996)
12. Kabsch, W.: A solution for the best rotation to relate two sets of vectors. Acta
Cryst. A32, 922–923 (1976)
13. Kabsch, W.: A discussion of the solution for the best rotation to relate two sets ofvectors. Acta Cryst. A34, 827–828 (1978)
14. Kallenberg, O.: Foundations of Modern Probability. Springer, Heidelberg (1997)
15. Kramers, H.A.: The behavior of macromolecules in inhomogeneous ﬂow. J. Chem.
Phys. 14(7), 415–424 (1946)

Searching Protein 3-D Structures in Linear Time

15

16. Schwartz, J.T., Sharir, M.: Identiﬁcation of partially obscured objects in two andthree dimensions by matching noisy characteristic curves. Intl. J. of Robotics Res. 6,
29–44 (1987)
17. Shibuya, T.: Geometric suﬃx tree: a new index structure for protein 3-D structures.
In: Lewenstein, M., Valiente, G. (eds.) CPM 2006. LNCS, vol. 4009, pp. 84–93.
Springer, Heidelberg (2006)
18. Shibuya, T.: Preﬁx-shuﬄed geometric suﬃx tree. In: Ziviani, N., Baeza-Yates, R.
(eds.) SPIRE 2007. LNCS, vol. 4726, pp. 300–309. Springer, Heidelberg (2007)
19. Shibuya, T.: Eﬃcient substructure RMSD query algorithms. J. Comput.
Biol. 14(9), 1201–1207 (2007)