Protein structure similarity from principle component correlation analysis

Materials

A total of fifty-six protein structures, grouped into 6 different sets according to CATH [23, 24] are used to test our algorithms. Proteins in structure set I belong to the "mainly alpha" class, including mostly apoptosis regulators in the BCL-x_L super family as well as others with remote conformational resemblance; all have the "Orthogonal Bundle" architecture. The atomic coordinates were retrieved from PDB with accession codes 1A4F, 1A6G, 1COL (A), 1DDB (A), 1F16 (A), 1G5M (A), 1GJH (A), 1MAZ, 1MDT (A), and 2BID (A), where (A) means chain A. Set II is also "mainly alpha" and has the same architecture as Set I, including structures 1CK7 (A), 1CXW (A), 1E8B (A), 1E88 (A), 1J7M (A), 1KS0 (A), 1PDC, and 2FN2. However this set consists of DNA helicase domains that have vastly different topology from Set I. Set II is used here to test the ability of our method to separate proteins that are in the same class of secondary structure but have different topologies. Set III belongs to the "mainly beta" class and has the barrel architecture, consisting of acid protease structures 1A5T, 1BVS (A), 1CUK, 1DV (A), 1F4I (A), 1G4A (E), 1G41 (A), 1HJP, 1IM2 (A), and 1JR3 (E). Set IV consists of the "alpha/beta" class proteins with the roll architecture, including structures 1FM0 (D), 1D4B (A), 1C78 (A), 1LM8 (B), 1NDD (A), 1UBQ, 1IBQ (A), and 1IP9 (A). The structures in set IV all have the Ubiquitin-like topology. Set V consists of the "mainly alpha" with the Alpha/alpha barrel architecture, including 1C82 (A), 1CB8 (A), 1EGU (A), 1F1S (A), 1F9G (A), 1HM2 (A), 1HM3 (A), 1HMU (A), 1HMW (A), 1HV6 (A), 1I8Q (A), and 1QAZ (A). The structures in Set V all have the Glycosyltransferase topology. Set VI consists of the "mainly beta" with the ribbon architecture, including 1AIW, 1E6N (A), 1E6P (A), 1E6R (A), 1E6Z (A), 1E15 (A), 1ED7 (A), and 1GOI (A). The structures in Set VI have the Seminal Fluid Protein PDC-109 (domain B).

Clustering of structurally similar proteins by SMEC method

One of the goals of this study is to compare and identify structurally or topologically similar proteins. In other words, given a new experimentally determined protein structure, the proposed method is expected to rapidly place the structure into a group of structurally or topologically similar proteins in the database, thereby aiding in correlating topological similarity with functional similarity. To illustrate the application of the SMEC approach, we compute the scaled eigenvalues of PD and PID interaction matrices (Section Methods). Figure 2a shows the plot of scaled λ₂ versus λ₁, calculated using the PD matrix, for all proteins in the four data sets. Figure 2b shows the plot of λ₁ of PID matrix versus that of PD matrices. The different symbols represent different structural groups. These plots were used to resolve clusters of structurally similar structures.

Pair-Wise structural comparison by PCC method

Principle component analysis of secondary interaction matrix

Assuming a protein having n secondary fragments denoted by h₁, h₂,..., h _n , and the number of residues in each secondary structure denoted by l₁, l₂,..., l _n , respectively, the total number of residues belonging to secondary structures is given by $N={\displaystyle \sum _{i=1}^n{l}_i}$. The invariant relation between a pair of secondary elements (h _i , h _j ) is described by a block matrix F(h _i , h _j ), in which the individual matrix elements represent a particular relation between residues of the two secondary structures. Since h _i has l_i residues (denoted by $c_i^1$, $c_i^2$,..., $c_i^{l_i}$), and h _j has l_j residues (denoted by $c_j^1$, $c_j^2$,..., $c_j^{l_j}$), the elements of the l _i × l _j F block matrix, g($c_i^u$, $c_j^v$), are defined as

$g(c_i^u,c_j^v)=\{\begin{array}{cc}d(c_i^u,c_j^v)& i\ne j\\ 0& i=j\end{array},\left(1\right)$

where 1 ≤ u ≤ l _i , 1 ≤ v ≤ l _j , and d($c_i^u$, $c_j^v$) is a real number representing an arbitrary invariant relation between residues of h _i and h _j . Note this approach allows the definition of d($c_i^u$, $c_j^v$) to be rather arbitrary. The full interaction matrix of a protein structure is square and symmetric and is defined as

$\widehat{I}={\left[\begin{array}{cccc}0& F(h_1,h_2)& \cdots & F(h_1,h_n)\\ F(h_2,h_1)& 0& \cdots & F(h_2,h_n)\\ ⋮& ⋮& \ddots & ⋮\\ F(h_n,h_1)& F(h_n,h_2)& \cdots & 0\end{array}\right]}_{N\times N}\left(2\right)$

The principle components of the interaction matrix is then obtained by orthogonal decomposition as shown below:

$\widehat{I}=E^T\left[\begin{array}{c}{\lambda }_1\\ {\lambda }_2\\ \ddots \\ {\lambda }_N\end{array}\right]E\left(3\right)$

where λ₁ ≥ λ₂ ≥ ⋯ ≥ λ _N are the sorted eigenvalues, the corresponding eigenvectors are e₁, e₂,..., e _N , and E = [e₁, e₂,..., e _N ] is an invertible matrix. Generally, the maximum eigenvalue, λ₁, and its corresponding eigenvector in N-dimensional space encode the most dominant features in the structure and therefore can be effectively used to directly compare structures, as well as to identify the less obvious topological features common to the proteins. Since the eigenvalues depend largely on the dimension of interaction matrix, they are divided by the matrix size N, a treatment similar to the scaling of writhing numbers in the SGM method (Rogen P. and Fain B., 2003). In a relatively crude analysis, λ₁ can be directly compared to infer structural similarity. This method is referred here as the Scaled Maximum Eigenvalue Comparison (SMEC).

In addition to the maximum eigenvalues, their corresponding eigenvectors can also be used to correlate similar structures. Particularly for pair-wise structure comparison, degree of similarity can be more accurately measured by comparing both eigenvalue and eigenvector. Since proteins are generally not of the same length, their eigenvectors cannot be directly correlated due to different dimensionality. Therefore, a "sliding window" approach is employed to correlate the smaller protein to all matching segments (length-wise) in the larger protein. Let us consider two proteins, A and B, having N and M secondary structure residues, respectively, and N ≤ M. For the protein having shorter secondary segments, λ^A and e^A are respectively the maximum eigenvalue and its corresponding N-dimensional eigenvector. For the protein with more secondary structure residues, M-N+1 interaction matrices are decomposed, where (λ^B₁, e^B₁) represent the principle components of the interaction matrix constructed from secondary structure residues 1 ... N, (λ^B₂, e^B₂) are from secondary structure residues 2 ... N+1, and so on. To quantify structural similarity, we define a difference metric, R, between Î of protein A and Î of the j th matching segment of protein B as

$R_j=\left|\right|e^A-e_j^B\left|\right|\left|{\lambda }^A-{\lambda }_j^B\right|\text{ , }1\le j\le M-N+1.\left(4\right)$

Obviously, smaller R _j indicates better correlation or higher degree of structural similarity. The overall difference between the two proteins is defined as

R = min(R₁, R₂,..., R_M-N+1).     (5)

The minimum of R₁, R₂, ..., R_M-N+1is used here to measure similarity because this potentially allows mapping a smaller structure onto a homologous domain within a larger protein. This method is called the Principle Component Correlation (PCC) analysis.

Defining the matrix elements

The definition of block matrix elements, d($c_i^u$, $c_j^v$), depends on the desired structural features to be extracted. In the current study, we focus structural comparison on protein backbone conformation. Clearly the simplest invariant describing the backbone conformation is the Euclidian distance between a pair of C^α atoms from two different secondary segments. Formally, the elements are defined as d($c_i^u$, $c_j^v$) = ||$c_i^u$ - $c_j^v$|| where $c_i^u$ and $c_j^v$ are the coordinates of the two C^α atoms of residues u of h_i and v of h_j, respectively. For conciseness, we name the interaction matrix so defined as the Pair-wise Distance (PD) matrix. For illustration purpose, the interaction matrix for the structure of Pb1, Domain of Bem1P (PDB accession code 1IP9), is shown in Fig. 1. This structure, consisting of two α helices and four β strands (Fig. 1a), is used here to provide distances between all pairs of C_α atoms in the six secondary elements (Fig. 1b).

Furthermore, two variations of the PD matrix definition are explored in attempt to provide a better resolution in structural comparison and classification. Since physical energy of interaction between a pair of atoms typically increase monotonically as the inverse of their separation, inverse of distance is used to mimic physical interactions between secondary elements. Here the elements of F(h _i , h _j ) are defined as

$d(c_i^u,c_j^v)=\{\begin{array}{cc}\frac{1}{\left|\left|c_i^u-c_j^v\right|\right|},& \left|\left|c_i^u-c_j^v\right|\right|\ge u_0\\ \frac{1}{u_0}& \left|\left|c_i^u-c_j^v\right|\right|<u_0\end{array}\left(6\right)$

where u₀ represent a hard-sphere boundary below which the interaction is constant. In this study, we arbitrarily set u₀ to 3Å. This definition is referred as Pair-wise Inverse Distance (PID) matrix.

Another variation of the PD matrix definition is to take into account the N – C terminal sense, in attempt to further emphasize protein topological features. For a secondary element, h _i , its direction vector v _i is defined by two points in Cartesian space: the center of mass of the five consecutive N-terminal C^α and the center of mass of the five consecutive C-terminal C^α atoms. Given a pair of secondary elements h _i and h _j , the new matrix elements are defined as

d($c_i^u$, $c_j^v$)' = d($c_i^u$, $c_j^v$)sgn(v _i ·v _j )     (7)

where sgn(x) is a symbol function which is 1 when x ≥ 0 and -1 when x < 0. This variation is referred as Pair-wise Distance with Sense (PDS) matrix in this study.

Linking/Writhing numbers

To evaluate the ability of PCC analysis in extracting pure topological features, the linking and writhing numbers, which are good measures of global topology, are also calculated for the four sets of structures for comparison. The linking number of two curves is defined by the Că lugă reanu-Fuller-White formula [25–27]: Lk = Wr + Tw, where the linking number Lk counts the sum of signed crossings between the ribbon's two boundary curves, the writhing number Wr counts the sum of signed self-crossings of the curve, averaged over all projection directions [28], and Tw is the twist number.Lk is an invariant to any smooth deformation that avoids self-intersections [29], and it is also independent of projection direction. Wr and Tw are invariant to some transformations, such as rigid body motions. Here we compute the writhing numbers using the Scaled Gauss Metric (SGM) approach previously described by Rogen and Fain [22].

Given two curves c₁ and c₂, which are two closed non-intersecting curves in 3-dimentional space, and define e(s, t) = (c₂(t) - c₁(s))/||c₂(t) - c₁(s)||, where ||·|| denotes the Euclidean norm. For two closed curves, the vector field e(s, t) is doubly periodic. Such mappings have an integer-valued degree that is invariant under topological deformations. The linking number of two curves is further defined as

$Lk(c_1,c{}_2)=\frac{1}{4\pi }{\displaystyle {\int }_{c_1}{\displaystyle {\int }_{c_2}\left[e,e_s,e_t\right]}}dsdt=\frac{1}{4\pi }{\displaystyle {\int }_{c_1}{\displaystyle {\int }_{c_2}\frac{\left(c_1^{\text{'}}(s)\times c_2^{\text{'}}(t)\right)\cdot \left(c_1(s)-c_2(t)\right)}{{\Vert c_1(s)-c_2(t)\Vert }^3}}}dsdt\left(8\right)$

where e _s and e _t are the tangents of e(s, t) at point (s, t), as well as $c_1^{\text{'}}$(s) and $c_2^{\text{'}}$(t) are the tangents along the c₁ and c₂ at s and t. Note that here e _s , e _t , $c_1^{\text{'}}$(s), and $c_2^{\text{'}}$(t) are vectors. Define w(s, t) = (c₁(t) - c₁(s)/||c₁(t) - c₁(s)||. The writhing number for a single curve c₁ is defined as

$Wr(c_1)=\frac{1}{4\pi }{\displaystyle {\int }_{c_1}{\displaystyle {\int }_{c_1}\left[w,w_s,w_t\right]}}dsdt=\frac{1}{4\pi }{\displaystyle {\int }_{c_1}{\displaystyle {\int }_{c_2}\frac{\left(c_1^{\text{'}}(s)\times c_1^{\text{'}}(t)\right)\cdot \left(c_1(s)-c_1(t)\right)}{{\Vert c_1(s)-c_1(t)\Vert }^3}}}dsdt\left(9\right)$

where w _s and w _t are the tangent of w(s, t) at point (s, t). Writhing number is not invariant under general smooth deformations such as translations, rotations, re-parameterizations, and dilations (Murasugi, 1996). Since the backbone of a protein is a polygonal curve, the writhing number of c₁(t) can be calculated by

$Wr(c_1)={\displaystyle \sum _{0<i_1<i_2<N}W(i_1,i_2),W(i_1,i_2)=\frac{1}{2\pi }{\displaystyle {\int }_{i_1=s}^{s+1}{\displaystyle {\int }_{i_2=t}^{t+1}w(s,t)dsdt}}}\left(10\right)$

where W(i₁, i₂) is the writhing number between the i₁ th and the i₂th segment; s and t denote two different C^α atoms, and N is the total number of C^α atoms. The SGM method is defined as the normalized writhing number, namely, Wr is divided by N [22]. The absolute difference between their writhing numbers is used to infer topological similarity.