1 // $Id: ssm_superpose.cpp $
2 // =================================================================
3 //
4 // CCP4 SSM Library: Protein structure superposition based on
5 // SSM algorithm: E. Krissinel & K. Henrick (2004) Acta Cryst. D60,
6 // 2256-2268.
7 //
8 // Copyright (C) Eugene Krissinel 2002-2013.
9 //
10 // This library is free software: you can redistribute it and/or
11 // modify it under the terms of the GNU Lesser General Public
12 // License version 3, modified in accordance with the provisions
13 // of the license to address the requirements of UK law.
14 //
15 // You should have received a copy of the modified GNU Lesser
16 // General Public License along with this library. If not, copies
17 // may be downloaded from http://www.ccp4.ac.uk/ccp4license.php
18 //
19 // This program is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU Lesser General Public License for more details.
23 //
24 // =================================================================
25 //
26 // 03.02.14 <-- Date of Last Modification.
27 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
28 // -----------------------------------------------------------------
29 //
30 // **** Module : ssm_spose <implementation>
31 // ~~~~~~~~~
32 // **** Functions : SuperposeCalphas ( superposing protein structures )
33 // ~~~~~~~~~~~
34 //
35 // E. Krissinel 2002-2014
36 //
37 // =================================================================
38 //
39
40
41 #include <string.h>
42 #include <math.h>
43
44 #include "mmdb2/mmdb_math_bfgsmin.h"
45 #include "ssm_superpose.h"
46
47 // =================================================================
48
49
CompatibleSSE(RSpAtom a)50 bool ssm::SpAtom::CompatibleSSE ( RSpAtom a ) {
51 if (sse==V_HELIX) return (a.sse==V_HELIX);
52 return (a.sse!=V_HELIX);
53 }
54
55
Copy(RSectionDist D)56 void ssm::SectionDist::Copy ( RSectionDist D ) {
57 dist = D.dist;
58 rmsd = D.rmsd;
59 cosine = D.cosine;
60 na = D.na;
61 core_pos1 = D.core_pos1;
62 core_pos2 = D.core_pos2;
63 core_e1 = D.core_e1;
64 core_e2 = D.core_e2;
65 pos1 = D.pos1;
66 pos2 = D.pos2;
67 e1 = D.e1;
68 e2 = D.e2;
69 sse1 = D.sse1;
70 sse2 = D.sse2;
71 }
72
73
74 // ==================================================================
75
Transform(mmdb::mat44 & T)76 void ssm::SSEDesc::Transform ( mmdb::mat44 & T ) {
77 x1 = T[0][0]*xs1 + T[0][1]*ys1 + T[0][2]*zs1 + T[0][3];
78 y1 = T[1][0]*xs1 + T[1][1]*ys1 + T[1][2]*zs1 + T[1][3];
79 z1 = T[2][0]*xs1 + T[2][1]*ys1 + T[2][2]*zs1 + T[2][3];
80 x2 = T[0][0]*xs2 + T[0][1]*ys2 + T[0][2]*zs2 + T[0][3];
81 y2 = T[1][0]*xs2 + T[1][1]*ys2 + T[1][2]*zs2 + T[1][3];
82 z2 = T[2][0]*xs2 + T[2][1]*ys2 + T[2][2]*zs2 + T[2][3];
83 }
84
CalcScore(RSSEDesc D)85 void ssm::SSEDesc::CalcScore ( RSSEDesc D ) {
86 // Don't change this procedure without appropriate changing
87 // ChooseFirstRotation().
88 mmdb::realtype dx,dy,dz;
89
90 dx = x1 - D.x1;
91 dy = y1 - D.y1;
92 dz = z1 - D.z1;
93 score = sqrt ( dx*dx + dy*dy + dz*dz );
94
95 dx = x2 - D.x2;
96 dy = y2 - D.y2;
97 dz = z2 - D.z2;
98 score += sqrt ( dx*dx + dy*dy + dz*dz );
99
100 dx = ((x1+x2) - (D.x1+D.x2))/2.0;
101 dy = ((y1+y2) - (D.y1+D.y2))/2.0;
102 dz = ((z1+z2) - (D.z1+D.z2))/2.0;
103 score += 2.0*sqrt(dx*dx+dy*dy+dz*dz);
104
105 score /= 16.0;
106 D.score = score;
107
108 }
109
110
GetVectorScore(mmdb::rvector r1,mmdb::rvector v1,mmdb::rvector r2,mmdb::rvector v2)111 mmdb::realtype GetVectorScore ( mmdb::rvector r1, mmdb::rvector v1,
112 mmdb::rvector r2, mmdb::rvector v2 ) {
113 // r1, r2: SSE origins
114 // v1, v2: SSE directions and lengths
115 mmdb::vect3 r0,v0;
116 mmdb::realtype vl2, t11,t12, t21,t22, dt1,dt2;
117 mmdb::realtype dr, dr11,dr12,dr21,dr22, dr2,dl2;
118 int i;
119
120 vl2 = 0.0;
121 t11 = 0.0;
122 t12 = 0.0;
123 t21 = 0.0;
124 t22 = 0.0;
125 for (i=0;i<3;i++) {
126 r0[i] = (r1[i]+r2[i])/2.0; // medium r1 and r2
127 v0[i] = (v1[i]+v2[i])/2.0; // medium v1 and v2
128 vl2 += v0[i]*v0[i]; // length of the medium
129 t11 += v0[i]*(r1[i]-r0[i]); // projection of r1 to medium
130 t12 += v0[i]*(r1[i]+v1[i]-r0[i]); // projection of r1+v1 to medium
131 t21 += v0[i]*(r2[i]-r0[i]); // projection of r2 to medium
132 t22 += v0[i]*(r2[i]+v2[i]-r0[i]); // projection of r2+v2 to medium
133 }
134
135 if (vl2<=0.0) return -1.0; // no score
136 if ((t12<t11) || (t22<t21)) return -2.0;
137
138 t11 /= vl2;
139 t12 /= vl2;
140 t21 /= vl2;
141 t22 /= vl2;
142
143 dr11 = 0.0;
144 dr12 = 0.0;
145 dr21 = 0.0;
146 dr22 = 0.0;
147 for (i=0;i<3;i++) {
148 dr = r1[i]-r0[i]; dr11 += dr*dr;
149 dr = r1[i]+v1[i]-r0[i]; dr12 += dr*dr;
150 dr = r2[i]-r0[i]; dr21 += dr*dr;
151 dr = r2[i]+v2[i]-r0[i]; dr22 += dr*dr;
152 }
153 dr11 -= t11*t11*vl2;
154 dr12 -= t12*t12*vl2;
155 dr21 -= t21*t21*vl2;
156 dr22 -= t22*t22*vl2;
157
158 dr2 = mmdb::RMax ( mmdb::RMax(dr11,dr12),mmdb::RMax(dr21,dr22) );
159
160 dt1 = mmdb::RMax ( t12-t11,t22-t21 );
161 dt2 = mmdb::RMax ( mmdb::RMax(t11,t12),mmdb::RMax(t21,t22) ) -
162 mmdb::RMin ( mmdb::RMin(t11,t12),mmdb::RMin(t21,t22) );
163
164 dl2 = dt2-dt1;
165 dl2 = dl2*dl2*vl2;
166
167 return dr2 + dl2;
168
169 }
170
Cosine(RSSEDesc D)171 mmdb::realtype ssm::SSEDesc::Cosine ( RSSEDesc D ) {
172 mmdb::realtype dx1,dy1,dz1, dx2,dy2,dz2, d1,d2,r;
173
174 dx1 = x2 - x1;
175 dy1 = y2 - y1;
176 dz1 = z2 - z1;
177 dx2 = D.x2 - D.x1;
178 dy2 = D.y2 - D.y1;
179 dz2 = D.z2 - D.z1;
180
181 d1 = dx1*dx1 + dy1*dy1 + dz1*dz1;
182 d2 = dx2*dx2 + dy2*dy2 + dz2*dz2;
183 r = d1*d2;
184
185 if (r>0.0) return (dx1*dx2+dy1*dy2+dz1*dz2)/sqrt(r);
186 else return 1.0;
187
188 }
189
Copy(RSSEDesc D)190 void ssm::SSEDesc::Copy ( RSSEDesc D ) {
191 x1 = D.x1; y1 = D.y1; z1 = D.z1;
192 x2 = D.x2; y2 = D.y2; z2 = D.z2;
193 xs1 = D.xs1; ys1 = D.ys1; zs1 = D.zs1;
194 xs2 = D.xs2; ys2 = D.ys2; zs2 = D.zs2;
195 score = D.score;
196 Qscore = D.Qscore;
197 Rscore = D.Rscore;
198 Xscore = D.Xscore;
199 pos = D.pos;
200 len = D.len;
201 pend = D.pend;
202 type = D.type;
203 classID = D.classID;
204 m = D.m;
205 match = D.match;
206 }
207
208
209 // ==================================================================
210
Superpose()211 ssm::Superpose::Superpose() {
212 InitSuperpose();
213 mmdb::GetMatrixMemory ( A ,3,3,1,1 );
214 mmdb::GetMatrixMemory ( U ,3,3,1,1 );
215 mmdb::GetMatrixMemory ( V ,3,3,1,1 );
216 mmdb::GetVectorMemory ( W ,3,1 );
217 mmdb::GetVectorMemory ( RV1,3,1 );
218 }
219
~Superpose()220 ssm::Superpose::~Superpose() {
221 FreeMemory();
222 mmdb::FreeMatrixMemory ( A ,3,1,1 );
223 mmdb::FreeMatrixMemory ( U ,3,1,1 );
224 mmdb::FreeMatrixMemory ( V ,3,1,1 );
225 mmdb::FreeVectorMemory ( W ,1 );
226 mmdb::FreeVectorMemory ( RV1,1 );
227 if (selString1) delete[] selString1;
228 if (selString2) delete[] selString2;
229 selString1 = NULL;
230 selString2 = NULL;
231 }
232
233
234
InitSuperpose()235 void ssm::Superpose::InitSuperpose() {
236
237 Rmsd0 = 3.0; // quality optimization parameter
238 minContact = 3.0; // minimal Calpha-pair contact parameter
239 maxContact = 5.0; // maximal Calpha-pair contact parameter
240 maxRMSD = 15.0; // maximal RMSD allowed
241 minCosine = 0.7; // minimum cosine between co-directional SSEs
242 SSEweight = 3.0; // currently does not work
243 sseGray = 4; // gray zone on the ends of SSEs allowed for
244 // matching to non-SSE atoms
245 shortSect1 = 1; // sections shorter than or equal to
246 // shortSect are removed
247 shortSect2 = 2;
248 iterMax = 1000;
249 iterMin = 10;
250 minQStep = 0.00001; // minimal quality improvement that counts
251 maxHollowIt = 10; // maximal allowed number of consequtive
252 // iterations without quality improvement
253 allowMC = false; // do not allow for misconnections
254
255 selInclHnd1 = 0; // selection handle for included Calpha1
256 selInclHnd2 = 0; // selection handle for included Calpha2
257
258 driverID = 0; // ID of the used Superpose driver
259
260 selString1 = NULL; // optional selection string for 1st structure
261 selString2 = NULL; // optional selection string for 2nd structure
262
263 a1 = NULL;
264 a2 = NULL;
265 nres1 = 0;
266 nres2 = 0;
267 nalgn = 0;
268 ngaps = 0;
269 nmd = 0;
270 rmsd_achieved = -1.0;
271 Q_achieved = -1.0;
272 nmisdr = 0;
273 seqIdent = 0.0;
274 ncombs = 1.0;
275
276 Calpha1 = NULL;
277 Calpha2 = NULL;
278
279 SSED1 = NULL;
280 SSED2 = NULL;
281 FH1 = NULL;
282 FS1 = NULL;
283 FH2 = NULL;
284 FS2 = NULL;
285 nSSEs1 = 0;
286 nSSEs2 = 0;
287 nFH1 = 0;
288 nFS1 = 0;
289 nFH2 = 0;
290 nFS2 = 0;
291
292 SDist = NULL;
293 SDistAlloc = 0;
294
295 }
296
FreeMemory()297 void ssm::Superpose::FreeMemory() {
298 int i;
299
300 if (a1) {
301 delete[] a1;
302 a1 = NULL;
303 }
304 if (a2) {
305 delete[] a2;
306 a2 = NULL;
307 }
308 nres1 = 0;
309 nres2 = 0;
310 nalgn = 0;
311 ngaps = 0;
312 rmsd_achieved = -1.0;
313 Q_achieved = -1.0;
314
315 if (SSED1) delete[] SSED1;
316 if (SSED2) delete[] SSED2;
317 SSED1 = NULL;
318 SSED2 = NULL;
319
320 mmdb::FreeVectorMemory ( FH1,1 );
321 mmdb::FreeVectorMemory ( FS1,1 );
322 mmdb::FreeVectorMemory ( FH2,1 );
323 mmdb::FreeVectorMemory ( FS2,1 );
324 nSSEs1 = 0;
325 nSSEs2 = 0;
326 nFH1 = 0;
327 nFS1 = 0;
328 nFH2 = 0;
329 nFS2 = 0;
330
331 if (SDist) {
332 for (i=0;i<SDistAlloc;i++)
333 if (SDist[i]) delete SDist[i];
334 delete[] SDist;
335 SDist = NULL;
336 }
337 SDistAlloc = 0;
338
339 }
340
SetAllowMC(bool allowMisconnections)341 void ssm::Superpose::SetAllowMC ( bool allowMisconnections ) {
342 allowMC = allowMisconnections;
343 }
344
SetIterationLimits(int iter_max,int iter_min,int max_hollow)345 void ssm::Superpose::SetIterationLimits ( int iter_max, int iter_min,
346 int max_hollow ) {
347 iterMax = iter_max;
348 iterMin = iter_min;
349 maxHollowIt = max_hollow;
350 }
351
SetCaSelections(mmdb::cpstr selection1,mmdb::cpstr selection2)352 void ssm::Superpose::SetCaSelections ( mmdb::cpstr selection1,
353 mmdb::cpstr selection2 ) {
354 mmdb::CreateCopy ( selString1,selection1 );
355 mmdb::CreateCopy ( selString2,selection2 );
356 }
357
GetTMatrix(mmdb::mat44 & TMat)358 void ssm::Superpose::GetTMatrix ( mmdb::mat44 & TMat ) {
359 int i,j;
360 for (i=0;i<4;i++)
361 for (j=0;j<4;j++)
362 TMat[i][j] = TMatrix[i][j];
363 }
364
GetTMatrix()365 mmdb::mat44 * ssm::Superpose::GetTMatrix() {
366 return &TMatrix;
367 }
368
GetRMSD()369 mmdb::realtype ssm::Superpose::GetRMSD() {
370 if (a1 && a2) return rmsd_achieved;
371 else return -1.0;
372 }
373
GetNAlign()374 int ssm::Superpose::GetNAlign() {
375 if (a1 && a2) return nalgn;
376 else return -1;
377 }
378
379
GetSuperposition(mmdb::ivector & Ca1,mmdb::rvector & dist1,int & nCa1,mmdb::ivector & Ca2,int & nCa2,mmdb::mat44 & TMat,mmdb::realtype & rmsdAchieved,int & nAligned,int & nGaps,mmdb::realtype & seqIdentity,int & nMisD,mmdb::realtype & nCombs)380 void ssm::Superpose::GetSuperposition ( mmdb::ivector & Ca1 ,
381 mmdb::rvector & dist1,
382 int & nCa1 ,
383 mmdb::ivector & Ca2 ,
384 int & nCa2 ,
385 mmdb::mat44 & TMat ,
386 mmdb::realtype & rmsdAchieved,
387 int & nAligned ,
388 int & nGaps ,
389 mmdb::realtype & seqIdentity ,
390 int & nMisD ,
391 mmdb::realtype & nCombs ) {
392 //
393 // Ca1[i]>=0 gives the index of Calpha of 2nd structure,
394 // superposed on ith Calpha of 1st structure (0<=i<nres1).
395 // If ith Calpha was not superposed then Ca1[i]<0.
396 //
397 // For superposed atoms
398 // Ca2[Ca1[i]] = i
399 // Ca1[Ca2[j]] = j
400 //
401 // Ca1 and Ca2 are dynamically allocated (indexed from 0 on)
402 // and should be disposed by the application. If Ca1 and/or Ca2
403 // are not initially set NULL, GetSuperposition attempts to
404 // dispose them first.
405 //
406 int i,j;
407 mmdb::FreeVectorMemory ( Ca1 ,0 );
408 mmdb::FreeVectorMemory ( dist1,0 );
409 mmdb::FreeVectorMemory ( Ca2 ,0 );
410 if (a1 && a2) {
411 mmdb::GetVectorMemory ( Ca1 ,nres1,0 );
412 mmdb::GetVectorMemory ( dist1,nres1,0 );
413 mmdb::GetVectorMemory ( Ca2 ,nres2,0 );
414 for (i=0;i<nres1;i++) {
415 Ca1[i] = a1[i].c0;
416 if (Ca1[i]>=0) dist1[i] = sqrt(a1[i].dist0);
417 else dist1[i] = -1.0;
418 }
419 for (i=0;i<nres2;i++)
420 Ca2[i] = a2[i].c0;
421 for (i=0;i<4;i++)
422 for (j=0;j<4;j++)
423 TMat[i][j] = TMatrix[i][j];
424 nCa1 = nres1;
425 nCa2 = nres2;
426 rmsdAchieved = rmsd_achieved;
427 nAligned = nalgn;
428 nGaps = ngaps;
429 seqIdentity = seqIdent;
430 nMisD = nmd;
431 nCombs = ncombs;
432 } else {
433 for (i=0;i<4;i++) {
434 for (j=0;j<4;j++)
435 TMat[i][j] = 0.0;
436 TMat[i][i] = 1.0;
437 }
438 nCa1 = 0;
439 nCa2 = 0;
440 rmsdAchieved = -1.0;
441 nAligned = 0;
442 nGaps = 0;
443 seqIdentity = 0.0;
444 nMisD = 0;
445 nCombs = 1.0;
446 }
447 }
448
449
GetCalphas1(mmdb::PPAtom & Calpha,int & numRes)450 void ssm::Superpose::GetCalphas1 ( mmdb::PPAtom & Calpha,
451 int & numRes ) {
452 Calpha = Calpha1;
453 numRes = nres1;
454 }
455
GetCalphas2(mmdb::PPAtom & Calpha,int & numRes)456 void ssm::Superpose::GetCalphas2 ( mmdb::PPAtom & Calpha,
457 int & numRes ) {
458 Calpha = Calpha2;
459 numRes = nres2;
460 }
461
462
GetSSEDesc1(RPSSEDesc sseDesc,int & numSSEs)463 void ssm::Superpose::GetSSEDesc1 ( RPSSEDesc sseDesc, int & numSSEs ) {
464 int i;
465 if (sseDesc) delete[] sseDesc;
466 sseDesc = new SSEDesc[nSSEs1];
467 for (i=0;i<nSSEs1;i++)
468 sseDesc[i].Copy ( SSED1[i] );
469 numSSEs = nSSEs1;
470 }
471
GetSSEDesc1()472 ssm::PSSEDesc ssm::Superpose::GetSSEDesc1() {
473 return SSED1;
474 }
475
GetSSEDesc2(RPSSEDesc sseDesc,int & numSSEs)476 void ssm::Superpose::GetSSEDesc2 ( RPSSEDesc sseDesc, int & numSSEs ) {
477 int i;
478 if (sseDesc) delete[] sseDesc;
479 sseDesc = new SSEDesc[nSSEs2];
480 for (i=0;i<nSSEs2;i++)
481 sseDesc[i].Copy ( SSED2[i] );
482 numSSEs = nSSEs2;
483 }
484
GetSSEDesc2()485 ssm::PSSEDesc ssm::Superpose::GetSSEDesc2() {
486 return SSED2;
487 }
488
GetSuperposedSSEs(mmdb::ivector v1,mmdb::ivector v2,int & nSupSSEs)489 void ssm::Superpose::GetSuperposedSSEs ( mmdb::ivector v1,
490 mmdb::ivector v2,
491 int & nSupSSEs ) {
492 int i;
493 nSupSSEs = 0;
494 for (i=0;i<nSSEs1;i++)
495 if (SSED1[i].match>0) {
496 nSupSSEs++;
497 v1[nSupSSEs] = i+1;
498 v2[nSupSSEs] = SSED1[i].match;
499 }
500 }
501
CalculateTMatrix()502 int ssm::Superpose::CalculateTMatrix() {
503 mmdb::realtype det,B;
504 int i,j,k;
505
506 det = A[1][1]*A[2][2]*A[3][3] +
507 A[1][2]*A[2][3]*A[3][1] +
508 A[2][1]*A[3][2]*A[1][3] -
509 A[1][3]*A[2][2]*A[3][1] -
510 A[1][1]*A[2][3]*A[3][2] -
511 A[3][3]*A[1][2]*A[2][1];
512
513 mmdb::math::SVD ( 3,3,3,A,U,V,W,RV1,true,true,i );
514
515 if (i!=0) {
516 for (j=0;j<4;j++) {
517 for (k=0;k<4;k++)
518 TMatrix[j][k] = 0.0;
519 TMatrix[j][j] = 1.0;
520 }
521 return 1;
522 }
523
524 if (det<=0.0) {
525 k = 0;
526 B = mmdb::MaxReal;
527 for (j=1;j<=3;j++)
528 if (W[j]<B) {
529 B = W[j];
530 k = j;
531 }
532 for (j=1;j<=3;j++)
533 V[j][k] = -V[j][k];
534 }
535
536 for (j=1;j<=3;j++)
537 for (k=1;k<=3;k++) {
538 B = 0.0;
539 for (i=1;i<=3;i++)
540 B += U[j][i]*V[k][i];
541 TMatrix[j-1][k-1] = B;
542 }
543
544 TMatrix[3][0] = 0.0;
545 TMatrix[3][1] = 0.0;
546 TMatrix[3][2] = 0.0;
547 TMatrix[3][3] = 1.0;
548
549 return 0;
550
551 }
552
553
554 DefineStructure(SMinTFunc)
555
556 struct SMinTFunc {
557 mmdb::vect3 * r1;
558 mmdb::vect3 * v1;
559 mmdb::vect3 * r2;
560 mmdb::vect3 * v2;
561 int nV;
562 void Init ();
563 void Allocate ( int nVectors );
564 void Dispose ();
565 };
566
Init()567 void SMinTFunc::Init() {
568 r1 = NULL; v1 = NULL;
569 r2 = NULL; v2 = NULL;
570 nV = 0;
571 }
572
Allocate(int nVectors)573 void SMinTFunc::Allocate ( int nVectors ) {
574 Dispose();
575 nV = nVectors;
576 if (nV>0) {
577 r1 = new mmdb::vect3[nV];
578 v1 = new mmdb::vect3[nV];
579 r2 = new mmdb::vect3[nV];
580 v2 = new mmdb::vect3[nV];
581 }
582 }
583
Dispose()584 void SMinTFunc::Dispose() {
585 if (r1) delete[] r1;
586 if (v1) delete[] v1;
587 if (r2) delete[] r2;
588 if (v2) delete[] v2;
589 Init();
590 }
591
MinTFunc(void * Data,int N,mmdb::rvector X,mmdb::realtype & F)592 void MinTFunc ( void * Data, int N, mmdb::rvector X,
593 mmdb::realtype & F ) {
594 UNUSED_ARGUMENT(N);
595 PSMinTFunc MTFD;
596 mmdb::mat33 erm;
597 mmdb::vect3 r1,v1;
598 mmdb::realtype rr,vv, score;
599 int i,j,k;
600
601 MTFD = PSMinTFunc(Data);
602
603 mmdb::GetEulerRotMatrix ( erm,X[4],X[5],X[6] );
604
605 F = 0.0;
606 for (i=0;i<MTFD->nV;i++) {
607 for (j=0;j<3;j++) {
608 rr = X[j+1];
609 vv = 0.0;
610 for (k=0;k<3;k++) {
611 rr += erm[j][k]*MTFD->r1[i][k];
612 vv += erm[j][k]*MTFD->v1[i][k];
613 }
614 r1[j] = rr;
615 v1[j] = vv;
616 }
617 score = GetVectorScore ( r1,v1,MTFD->r2[i],MTFD->v2[i] );
618 if (score>0.0) F += score;
619 }
620 }
621
SuperposeSSGraphs(PGraph G1,mmdb::ivector F1,PGraph G2,mmdb::ivector F2,int matchlen)622 int ssm::Superpose::SuperposeSSGraphs ( PGraph G1, mmdb::ivector F1,
623 PGraph G2, mmdb::ivector F2,
624 int matchlen ) {
625 mmdb::math::PBFGSMin BFGS;
626 SMinTFunc SMTF;
627 PVertex Vx;
628 //mat33 erm,trm;
629 mmdb::vect3 v1,v2;
630 mmdb::realtype T[7],TypT[7];
631 mmdb::realtype B,B1, x01,y01,z01, x02,y02,z02, mass,mass1,mass2;
632 int i,j,k,l, nE1,nE2;
633 bool AddEdges1,AddEdges2;
634
635 nE1 = G1->GetNofEdges();
636 if (!nE1) G1->BuildGraph();
637
638 nE2 = G2->GetNofEdges();
639 if (!nE2) G2->BuildGraph();
640
641 for (j=1;j<=3;j++)
642 for (k=1;k<=3;k++)
643 A[j][k] = 0.0;
644
645 for (i=1;i<=matchlen;i++) {
646 Vx = G1->GetGraphVertex ( F1[i] );
647 Vx->GetDirection ( v1 );
648 B = Vx->GetMass();
649 Vx = G2->GetGraphVertex ( F2[i] );
650 Vx->GetDirection ( v2 );
651 B += Vx->GetMass();
652 for (j=1;j<=3;j++)
653 for (k=1;k<=3;k++)
654 A[j][k] += B*v1[k-1]*v2[j-1];
655 }
656
657 AddEdges1 = true;
658 for (i=1;(i<matchlen) && AddEdges1;i++) {
659 Vx = G1->GetGraphVertex ( F1[i] );
660 for (l=i+1;(l<=matchlen) && AddEdges1;l++) {
661 B = Vx->GetCosine ( G1->GetGraphVertex(F1[l]) );
662 if (fabs(B)<0.8) AddEdges1 = false;
663 }
664 }
665
666 AddEdges2 = true;
667 for (i=1;(i<matchlen) && AddEdges2;i++) {
668 Vx = G2->GetGraphVertex ( F2[i] );
669 for (l=i+1;(l<=matchlen) && AddEdges2;l++) {
670 B = Vx->GetCosine ( G2->GetGraphVertex(F2[l]) );
671 if (fabs(B)<0.8) AddEdges2 = false;
672 }
673 }
674
675 if (AddEdges1 || AddEdges2)
676 for (i=1;i<matchlen;i++)
677 for (l=i+1;l<=matchlen;l++)
678 if (G1->GetEdgeDirection(F1[i],F1[l],v1) &&
679 G2->GetEdgeDirection(F2[i],F2[l],v2)) {
680 B = (G1->GetMass(F1[i]) + G1->GetMass(F1[l]) +
681 G2->GetMass(F2[i]) + G2->GetMass(F2[l]))/2.0;
682 for (j=1;j<=3;j++)
683 for (k=1;k<=3;k++)
684 A[j][k] += B*v1[k-1]*v2[j-1];
685 }
686
687
688 if (CalculateTMatrix()) return 1;
689
690
691 // 9. Add translation
692 x01 = 0.0; y01 = 0.0; z01 = 0.0; mass1 = 0.0;
693 x02 = 0.0; y02 = 0.0; z02 = 0.0; mass2 = 0.0;
694 for (i=1;i<=matchlen;i++) {
695 Vx = G1->GetGraphVertex ( F1[i] );
696 mass = Vx->GetMass();
697 Vx->GetPosition ( v1 );
698 x01 += v1[0]*mass;
699 y01 += v1[1]*mass;
700 z01 += v1[2]*mass;
701 mass1 += mass;
702 Vx = G2->GetGraphVertex ( F2[i] );
703 mass = Vx->GetMass();
704 Vx->GetPosition ( v2 );
705 x02 += v2[0]*mass;
706 y02 += v2[1]*mass;
707 z02 += v2[2]*mass;
708 mass2 += mass;
709 }
710 x01 /= mass1; y01 /= mass1; z01 /= mass1;
711 x02 /= mass2; y02 /= mass2; z02 /= mass2;
712
713 // 10. Optimize translation
714
715 for (i=0;i<7;i++) {
716 T [i] = 0.0;
717 TypT[i] = 1.0;
718 }
719
720 SMTF.Init();
721 SMTF.Allocate ( matchlen );
722 // Bring both structures' mass centers into (0,0,0) and
723 // superpose
724 for (i=0;i<matchlen;i++) {
725 Vx = G1->GetGraphVertex ( F1[i+1] );
726 v1[0] = Vx->GetX1();
727 v1[1] = Vx->GetY1();
728 v1[2] = Vx->GetZ1();
729 v2[0] = Vx->GetX2() - v1[0];
730 v2[1] = Vx->GetY2() - v1[1];
731 v2[2] = Vx->GetZ2() - v1[2];
732 v1[0] -= x01;
733 v1[1] -= y01;
734 v1[2] -= z01;
735 for (j=0;j<3;j++) {
736 B = 0.0;
737 B1 = 0.0;
738 for (k=0;k<3;k++) {
739 B += TMatrix[j][k]*v1[k];
740 B1 += TMatrix[j][k]*v2[k];
741 }
742 SMTF.r1[i][j] = B;
743 SMTF.v1[i][j] = B1;
744 }
745 Vx = G2->GetGraphVertex ( F2[i+1] );
746 SMTF.r2[i][0] = Vx->GetX1();
747 SMTF.r2[i][1] = Vx->GetY1();
748 SMTF.r2[i][2] = Vx->GetZ1();
749 SMTF.v2[i][0] = Vx->GetX2() - SMTF.r2[i][0];
750 SMTF.v2[i][1] = Vx->GetY2() - SMTF.r2[i][1];
751 SMTF.v2[i][2] = Vx->GetZ2() - SMTF.r2[i][2];
752 SMTF.r2[i][0] -= x02;
753 SMTF.r2[i][1] -= y02;
754 SMTF.r2[i][2] -= z02;
755 }
756
757 BFGS = new mmdb::math::BFGSMin();
758
759 BFGS->SetMinFunction ( &SMTF,MinTFunc );
760
761 BFGS->BFGS_Driver ( 6,T,TypT,B,i, // MinN,x0,TypX,FuncValue,TermCode
762 0,100,1.0, // Digits,ItnLmt,TypF
763 0.0,0.0,mmdb::MaxReal, // GrdTol,StpTol,MaxStp
764 false,NULL,NULL // Hess,LowLimit,TopLimit
765 );
766
767 delete BFGS;
768
769 SMTF.Dispose();
770
771 /*
772 GetEulerRotMatrix ( erm,T[4],T[5],T[6] );
773
774 for (i=0;i<3;i++)
775 for (j=0;j<3;j++) {
776 B = 0.0;
777 for (k=0;k<3;k++)
778 B += erm[i][k]*TMatrix[k][j];
779 trm[i][j] = B;
780 }
781 */
782
783 TMatrix[0][3] = x02 - TMatrix[0][0]*x01 - TMatrix[0][1]*y01 -
784 TMatrix[0][2]*z01 - T[1];
785 TMatrix[1][3] = y02 - TMatrix[1][0]*x01 - TMatrix[1][1]*y01 -
786 TMatrix[1][2]*z01 - T[2];
787 TMatrix[2][3] = z02 - TMatrix[2][0]*x01 - TMatrix[2][1]*y01 -
788 TMatrix[2][2]*z01 - T[3];
789
790 if (!nE1) G1->ReleaseEdges();
791 if (!nE2) G2->ReleaseEdges();
792
793 return 0;
794
795 }
796
797
SelectCalphas(mmdb::PManager MMDB,PGraph G,mmdb::PPAtom & Calpha,PSpAtom & a,int & nres,int & selHnd,int selInclHnd,mmdb::cpstr selString)798 void ssm::Superpose::SelectCalphas ( mmdb::PManager MMDB, PGraph G,
799 mmdb::PPAtom & Calpha, PSpAtom & a,
800 int & nres, int & selHnd,
801 int selInclHnd,
802 mmdb::cpstr selString ) {
803 int i;
804
805 if (a && (driverID!=2)) {
806 delete[] a;
807 a = NULL;
808 }
809
810 if ((selHnd<=0) || (driverID!=2)) {
811 G->SelectCalphas ( MMDB,selHnd,selString );
812 MMDB->GetSelIndex ( selHnd,Calpha,nres );
813 }
814
815 if (nres>0) {
816 if (!a) {
817 a = new SpAtom[nres];
818 for (i=0;i<nres;i++) {
819 strcpy ( a[i].chID,Calpha[i]->GetChainID() );
820 a[i].sse = V_UNKNOWN;
821 a[i].c0 = -1;
822 a[i].dist = -1.0;
823 a[i].dist0 = -1.0;
824 if (selInclHnd>0)
825 a[i].excluded = !Calpha[i]->isInSelection ( selInclHnd );
826 else a[i].excluded = false;
827 }
828 } else {
829 for (i=0;i<nres;i++) {
830 a[i].c0 = -1;
831 a[i].dist = -1.0;
832 a[i].dist0 = -1.0;
833 if (selInclHnd>0)
834 a[i].excluded = !Calpha[i]->isInSelection ( selInclHnd );
835 else a[i].excluded = false;
836 }
837 }
838 }
839
840 }
841
842
MapSSEs(mmdb::PPAtom Calpha,PSpAtom a,int nres,PGraph G,RPSSEDesc SSED,int & nSSEs)843 void ssm::Superpose::MapSSEs ( mmdb::PPAtom Calpha, PSpAtom a, int nres,
844 PGraph G, RPSSEDesc SSED,
845 int & nSSEs ) {
846 // Returns vector SSED.pos containing the origin of SSE sequences in
847 // array a, and vector SSED.len containing the SSEs' lengths.
848 // nSSEs returns the number of SSEs.
849 PVertex V;
850 mmdb::ChainID chID;
851 int initSeqNum,endSeqNum,seqNum;
852 mmdb::InsCode initICode,endICode;
853 mmdb::pstr iCode;
854 int i,j,vtype,k, p1,p2, c1,c2;
855
856 if (SSED && (driverID!=2)) {
857 delete[] SSED;
858 SSED = NULL;
859 }
860
861 nSSEs = G->GetNofVertices();
862 if (nSSEs<=0) return;
863
864 if (!SSED) {
865
866 SSED = new SSEDesc[nSSEs];
867
868 for (i=1;i<=nSSEs;i++) {
869 G->GetVertexRange ( i,chID,initSeqNum,initICode,
870 endSeqNum,endICode );
871 vtype = G->GetVertexType(i);
872 p1 = -1;
873 p2 = -1;
874 c1 = -1;
875 c2 = -1;
876 for (j=0;(j<nres) && ((p2<0) || (p1<0));j++) {
877 if (!strcmp(chID,Calpha[j]->GetChainID())) {
878 if (c1<0) c1 = j;
879 c2 = j;
880 seqNum = Calpha[j]->GetSeqNum ();
881 iCode = Calpha[j]->GetInsCode();
882 if ((p1<0) && (initSeqNum==seqNum) &&
883 (!strcmp(initICode,iCode))) p1 = j;
884 if ((p2<0) && (endSeqNum==seqNum) &&
885 (!strcmp(endICode,iCode))) p2 = j;
886 }
887 }
888 k = i-1;
889 if ((p1>=0) && (p2<0)) p2 = c2;
890 if ((p1<0) && (p2>=0)) p1 = c1;
891 if ((p1>=0) && (p2>=0)) {
892 if (p2<p1) {
893 j = p1;
894 p1 = p2;
895 p2 = j;
896 }
897 SSED[k].pos = p1;
898 SSED[k].pend = p2;
899 SSED[k].len = p2-p1+1;
900 for (j=p1+sseGray;j<=p2-sseGray;j++)
901 a[j].sse = vtype;
902 } else {
903 SSED[k].pos = -1;
904 SSED[k].pend = -1;
905 SSED[k].len = 0;
906 }
907 V = G->GetGraphVertex ( i );
908 if (V) {
909 SSED[k].x1 = V->GetX1();
910 SSED[k].x2 = V->GetX2();
911 SSED[k].y1 = V->GetY1();
912 SSED[k].y2 = V->GetY2();
913 SSED[k].z1 = V->GetZ1();
914 SSED[k].z2 = V->GetZ2();
915 } else {
916 SSED[k].x1 = 0.0;
917 SSED[k].x2 = 0.0;
918 SSED[k].y1 = 0.0;
919 SSED[k].y2 = 0.0;
920 SSED[k].z1 = 0.0;
921 SSED[k].z2 = 0.0;
922 }
923 SSED[k].xs1 = SSED[k].x1;
924 SSED[k].xs2 = SSED[k].x2;
925 SSED[k].ys1 = SSED[k].y1;
926 SSED[k].ys2 = SSED[k].y2;
927 SSED[k].zs1 = SSED[k].z1;
928 SSED[k].zs2 = SSED[k].z2;
929 SSED[k].type = vtype;
930 SSED[k].classID = G->GetVertexClass(i);
931 SSED[k].score = mmdb::MaxReal;
932 SSED[k].Qscore = 0.0;
933 SSED[k].Rscore = 0.0;
934 SSED[k].Xscore = 0.0;
935 SSED[k].m = -1;
936 SSED[k].match = -1;
937 }
938
939 } else {
940
941 for (i=0;i<nSSEs;i++) {
942 SSED[i].score = mmdb::MaxReal;
943 SSED[i].Qscore = 0.0;
944 SSED[i].Rscore = 0.0;
945 SSED[i].Xscore = 0.0;
946 SSED[i].m = -1;
947 SSED[i].match = -1;
948 }
949
950 }
951
952 }
953
954
IdentifyUnmatchedSSEs(mmdb::ivector & FH,int & nFH,mmdb::ivector & FS,int & nFS,mmdb::ivector F,int mlen,PGraph G)955 void ssm::Superpose::IdentifyUnmatchedSSEs (
956 mmdb::ivector & FH, int & nFH,
957 mmdb::ivector & FS, int & nFS,
958 mmdb::ivector F, int mlen,
959 PGraph G ) {
960 int i,j,k,nSSEs;
961
962 mmdb::FreeVectorMemory ( FH,1 );
963 mmdb::FreeVectorMemory ( FS,1 );
964
965 nSSEs = G->GetNofVertices();
966 if (nSSEs<=0) return;
967
968 mmdb::GetVectorMemory ( FH,nSSEs,1 );
969 mmdb::GetVectorMemory ( FS,nSSEs,1 );
970 nFH = 0;
971 nFS = 0;
972 for (i=1;i<=nSSEs;i++) {
973 k = 0;
974 for (j=1;(j<=mlen) && (!k);j++)
975 if (F[j]==i) k = j;
976 if (!k) {
977 if (G->GetVertexType(i)==V_HELIX) FH[++nFH] = i;
978 else FS[++nFS] = i;
979 }
980 }
981
982 }
983
984
CalcDistance(int SSE1,int SSE2,RSectionDist D)985 void ssm::Superpose::CalcDistance ( int SSE1, int SSE2,
986 RSectionDist D ) {
987 // Calculates the minimal distance between SSEs number SSE1
988 // and SSE2 from pos1 to pos1+len1-1 and pos2 to pos2+len2-1,
989 // as given in the SSE description arrays, SSED1 and SSED2,
990 // respectively. The distance is returned in D together with
991 // other parameters needed for 3D alignment of the SSEs.
992 mmdb::realtype d,d1;
993 int minAlign, pos1,pos2, len1,len2, i,j,k,l, i1,i2;
994 int p1,p2;
995
996
997 // 1. Initial preparations
998
999 i1 = SSE1-1;
1000 i2 = SSE2-1;
1001 pos1 = SSED1[i1].pos;
1002 pos2 = SSED2[i2].pos;
1003 len1 = SSED1[i1].len;
1004 len2 = SSED2[i2].len;
1005
1006 // store indices of SSEs for future references
1007 D.sse1 = SSE1;
1008 D.sse2 = SSE2;
1009
1010 if ((len1<=0) || (len2<=0)) {
1011 D.dist = mmdb::MaxReal;
1012 D.rmsd = mmdb::MaxReal;
1013 D.cosine = -1.0;
1014 D.core_pos1 = -1;
1015 D.core_pos2 = -1;
1016 D.core_e1 = -1;
1017 D.core_e2 = -1;
1018 D.na = 0;
1019 D.pos1 = -1;
1020 D.pos2 = -1;
1021 D.e1 = -1;
1022 D.e2 = -1;
1023 return;
1024 }
1025
1026 // we require that at least minAlign Calphas from each SSE
1027 // are corresponded
1028 if (SSED1[i1].type==V_HELIX) minAlign = 4; //hx_min_len;
1029 else minAlign = 3; //sd_min_len;
1030
1031 minAlign = mmdb::IMin ( minAlign,len1 );
1032 minAlign = mmdb::IMin ( minAlign,len2 );
1033
1034 // calculate cosine between SSE main vectors
1035 D.cosine = SSED1[i1].Cosine ( SSED2[i2] );
1036
1037
1038 // 2. Calculate the square distance matrix
1039
1040 i1 = pos1;
1041 for (i=0;i<len1;i++) {
1042 i2 = pos2;
1043 for (j=0;j<len2;j++)
1044 AD[i][j] = Calpha1[i1]->GetDist2 ( Calpha2[i2++] );
1045 i1++;
1046 }
1047
1048
1049 // 3. Look for minAlign-long section with minimal rmsd
1050
1051 d = mmdb::MaxReal;
1052 p1 = -1;
1053 p2 = -1;
1054
1055 for (i=0;i<=len1-minAlign;i++) {
1056 l = mmdb::IMin(len1-i,len2) - minAlign;
1057 for (j=0;j<=l;j++) {
1058 d1 = 0.0;
1059 for (k=j;k<j+minAlign;k++)
1060 d1 += AD[i+k][k];
1061 if (d1<d) {
1062 d = d1;
1063 p1 = i + j;
1064 p2 = j;
1065 }
1066 }
1067 }
1068
1069 for (j=0;j<=len2-minAlign;j++) {
1070 l = mmdb::IMin(len2-j,len1) - minAlign;
1071 for (i=0;i<=l;i++) {
1072 d1 = 0.0;
1073 for (k=i;k<i+minAlign;k++)
1074 d1 += AD[k][j+k];
1075 if (d1<d) {
1076 d = d1;
1077 p1 = i;
1078 p2 = j+i;
1079 }
1080 }
1081 }
1082
1083 D.core_pos1 = pos1 + p1;
1084 D.core_pos2 = pos2 + p2;
1085 D.core_e1 = D.core_pos1 + minAlign - 1;
1086 D.core_e2 = D.core_pos2 + minAlign - 1;
1087
1088 if (p1>=0) D.na = minAlign;
1089 else D.na = 0;
1090
1091 D.dist = d;
1092 D.rmsd = d/minAlign;
1093
1094 // 4. Expand the initial alignment section into both directions
1095
1096 l = mmdb::IMin ( p1,p2 );
1097 D.pos1 = D.core_pos1 - l;
1098 D.pos2 = D.core_pos2 - l;
1099
1100 l = mmdb::IMin ( pos1+len1-D.core_e1,pos2+len2-D.core_e2 ) - 1;
1101 D.e1 = D.core_e1 + l;
1102 D.e2 = D.core_e2 + l;
1103
1104 }
1105
1106
AlignSSEs(RSectionDist D,int unmap)1107 void ssm::Superpose::AlignSSEs ( RSectionDist D, int unmap ) {
1108 // Makes correspondence between sections of atoms in arrays Calpha1
1109 // and Calpha2 from pos1 to pos1+len1-1 and pos2 to pos2+len2-1,
1110 // respectively using data obtained from calcDistance.
1111 // The correspondence is returned as pairs c2[c1[i]]=i and
1112 // c1[c2[j]]=j for atoms i and j; c1[k] and c2[m] for
1113 // non-corresponding atoms k and m are not changed.
1114 int i1,i2, c1;
1115
1116 if (D.na>0) {
1117
1118 i1 = D.pos1;
1119 i2 = D.pos2;
1120 c1 = D.core_pos1;
1121 if (unmap!=UNMAP_NO) c1 = (c1+D.core_e1)/2;
1122
1123 while (i1<c1) {
1124 a1[i1].c = i2;
1125 a2[i2].c = i1;
1126 if (i1>D.pos1) {
1127 a1[i1].unmap1 = D.pos1; // do not unmap
1128 a1[i1].unmap2 = i1-1; // until any of
1129 a2[i2].unmap1 = D.pos2; // previous atoms
1130 a2[i2].unmap2 = i2-1; // is unmapped
1131 }
1132 a1[i1].dist = Calpha1[i1]->GetDist2 ( Calpha2[i2] );
1133 a2[i2].dist = a1[i1].dist;
1134 i1++;
1135 i2++;
1136 }
1137
1138 if (unmap==UNMAP_NO)
1139 while (i1<=D.core_e1) {
1140 a1[i1].c = i2;
1141 a2[i2].c = i1;
1142 a1[i1].unmap1 = unmap;
1143 a1[i1].unmap2 = unmap;
1144 a2[i2].unmap1 = unmap;
1145 a2[i2].unmap2 = unmap;
1146 a1[i1].dist = Calpha1[i1]->GetDist2 ( Calpha2[i2] );
1147 a2[i2].dist = a1[i1].dist;
1148 i1++;
1149 i2++;
1150 }
1151
1152 while (i1<=D.e1) {
1153 a1[i1].c = i2;
1154 a2[i2].c = i1;
1155 if (i1<D.e1) {
1156 a1[i1].unmap1 = i1+1; // do not unmap
1157 a1[i1].unmap2 = D.e1; // until any of
1158 a2[i2].unmap1 = i2+1; // next atoms
1159 a2[i2].unmap2 = D.e2; // is unmapped
1160 }
1161 a1[i1].dist = Calpha1[i1]->GetDist2 ( Calpha2[i2] );
1162 a2[i2].dist = a1[i1].dist;
1163 i1++;
1164 i2++;
1165 }
1166
1167 SSED1[D.sse1-1].m = D.sse2;
1168 SSED2[D.sse2-1].m = D.sse1;
1169
1170 }
1171
1172 }
1173
isMC(int pos1,int pos2)1174 bool ssm::Superpose::isMC ( int pos1, int pos2 ) {
1175 // Returns true if matching the Calpha pair in the positions
1176 // (pos1,pos2) of the chains would contradict to the already
1177 // aligned pairs and allowMC is set false.
1178 int i;
1179
1180 if (allowMC) return false;
1181
1182 i = pos1 + 1;
1183 while (i<nres1)
1184 if (a1[i].c>=0) break;
1185 else i++;
1186 if (i<nres1) {
1187 if (pos2>=a1[i].c) {
1188 if ((!strcmp(a1[pos1].chID,a1[i].chID)) &&
1189 (!strcmp(a2[pos2].chID,a2[a1[i].c].chID)))
1190 return true;
1191 }
1192 }
1193
1194 i = pos1 - 1;
1195 while (i>=0)
1196 if (a1[i].c>=0) break;
1197 else i--;
1198 if (i>=0) {
1199 if (pos2<=a1[i].c) {
1200 if ((!strcmp(a1[pos1].chID,a1[i].chID)) &&
1201 (!strcmp(a2[pos2].chID,a2[a1[i].c].chID)))
1202 return true;
1203 }
1204 }
1205
1206 return false;
1207
1208 }
1209
CorrespondSSEs(mmdb::ivector F1,int nF1,mmdb::ivector F2,int nF2,mmdb::realtype rmsd_est)1210 void ssm::Superpose::CorrespondSSEs ( mmdb::ivector F1, int nF1,
1211 mmdb::ivector F2, int nF2,
1212 mmdb::realtype rmsd_est ) {
1213 mmdb::realtype rmsd,rmsd2;
1214 int i,j, k1,k2, i1,j1;
1215
1216 rmsd2 = rmsd_est*rmsd_est;
1217
1218 for (i=1;i<=nF1;i++) {
1219 i1 = i-1;
1220 k1 = F1[i]-1;
1221 for (j=1;j<=nF2;j++) {
1222 j1 = j-1;
1223 k2 = F2[j]-1;
1224 if ((SSED1[k1].type==SSED2[k2].type) &&
1225 (SSED1[k1].classID==SSED2[k2].classID) &&
1226 (!isMC(SSED1[k1].pos,SSED2[k2].pos))) {
1227 CalcDistance ( F1[i],F2[j],SDist[i1][j1] );
1228 if ((SDist[i1][j1].na<=0) ||
1229 (SDist[i1][j1].cosine<minCosine) ||
1230 (SDist[i1][j1].rmsd>rmsd2))
1231 SDist[i1][j1].rmsd = -1.0;
1232 } else
1233 SDist[i1][j1].rmsd = -1.0;
1234 }
1235 }
1236
1237 do {
1238 rmsd = mmdb::MaxReal;
1239 i1 = -1;
1240 j1 = -1;
1241 for (i=0;i<nF1;i++)
1242 for (j=0;j<nF2;j++)
1243 if ((SDist[i][j].rmsd>=0.0) && (SDist[i][j].rmsd<rmsd) &&
1244 (!isMC(SDist[i][j].pos1,SDist[i][j].pos2))) {
1245 rmsd = SDist[i][j].rmsd;
1246 i1 = i;
1247 j1 = j;
1248 }
1249 if (i1>=0) {
1250 AlignSSEs ( SDist[i1][j1],UNMAP_YES );
1251 for (j=0;i<nF2;i++)
1252 SDist[i1][j].rmsd = -1.0;
1253 for (i=0;i<nF1;i++)
1254 SDist[i][j1].rmsd = -1.0;
1255 }
1256 } while (i1>=0);
1257
1258 }
1259
1260
ExpandContact(mmdb::RContact c,int & ip,int & im,mmdb::realtype maxDist2)1261 void ssm::Superpose::ExpandContact ( mmdb::RContact c,
1262 int & ip, int & im,
1263 mmdb::realtype maxDist2 ) {
1264 // ip expands contact (c.id1,c.id2) in positive direction of indices
1265 // im expands contact (c.id1,c.id2) in negative direction of indices
1266 // once ip reaches end of the chain, it is assigned -1
1267 // once im reaches begining of the chain, it is assigned -1
1268 mmdb::realtype dist;
1269 int i1,i2;
1270
1271 if (ip>=0) {
1272 i1 = c.id1 + ip;
1273 i2 = c.id2 + ip;
1274 if (!isMC(i1,i2)) {
1275 while ((i1<nres1) && (i2<nres2))
1276 if ((a1[i1].c<0) && (a2[i2].c<0) &&
1277 a1[i1].CompatibleSSE(a2[i2])) {
1278 dist = Calpha1[i1]->GetDist2 ( Calpha2[i2] );
1279 if (dist<maxDist2) {
1280 a1[i1].c = i2;
1281 a2[i2].c = i1;
1282 a1[i1].dist = dist;
1283 a2[i2].dist = dist;
1284 i1++;
1285 i2++;
1286 } else
1287 break;
1288 } else
1289 i1 = nres1;
1290 } else
1291 i1 = nres1;
1292 if ((i1<nres1) && (i2<nres2)) ip = i1 - c.id1;
1293 else ip = -1;
1294 }
1295
1296 if (im>=0) {
1297 i1 = c.id1 - im;
1298 i2 = c.id2 - im;
1299 if (!isMC(i1,i2)) {
1300 while ((i1>=0) && (i2>=0))
1301 if ((a1[i1].c<0) && (a2[i2].c<0) &&
1302 a1[i1].CompatibleSSE(a2[i2])) {
1303 dist = Calpha1[i1]->GetDist2 ( Calpha2[i2] );
1304 if (Calpha1[i1]->GetDist2(Calpha2[i2])<maxDist2) {
1305 a1[i1].c = i2;
1306 a2[i2].c = i1;
1307 a1[i1].dist = dist;
1308 a2[i2].dist = dist;
1309 i1--;
1310 i2--;
1311 } else
1312 break;
1313 } else
1314 i1 = -1;
1315 } else
1316 i1 = -1;
1317 if ((i1>=0) && (i2>=0)) im = c.id1 - i1;
1318 else im = -1;
1319 }
1320
1321 }
1322
1323
1324
CorrespondContacts(mmdb::PManager M1,mmdb::realtype rmsd_est)1325 void ssm::Superpose::CorrespondContacts ( mmdb::PManager M1,
1326 mmdb::realtype rmsd_est ) {
1327 // Find the closest contacts between yet unmapped atoms
1328 // and gradually expand them
1329 mmdb::PContact contact;
1330 mmdb::ivector ip,im;
1331 mmdb::vect3 v1,v2;
1332 mmdb::realtype l1,l2,cosine;
1333 int ncontacts,i,j,m1,m2,n1,n2,i1,i2;
1334
1335 // 1. Find all contacts in the range og 0.0 - Rmsd
1336 contact = NULL;
1337 ncontacts = 0;
1338 M1->SeekContacts ( Calpha1,nres1,Calpha2,nres2,0.0,rmsd_est,0,
1339 contact,ncontacts,0,NULL,0,
1340 mmdb::BRICK_ON_2 | mmdb::BRICK_READY );
1341
1342 // 2. Leave only open contacts (i.e. those between yet unmapped
1343 // atoms), with regard to the SSE type
1344 j = 0;
1345 for (i=0;i<ncontacts;i++) {
1346 i1 = contact[i].id1;
1347 i2 = contact[i].id2;
1348 if ((!a1[i1].excluded) && (a1[i1].c<0) &&
1349 (!a2[i2].excluded) && (a2[i2].c<0) &&
1350 a1[i1].CompatibleSSE(a2[i2])) {
1351 if (j<i) contact[j].Copy ( contact[i] );
1352 contact[j].dist = contact[j].dist*contact[j].dist;
1353 j++;
1354 }
1355 }
1356 ncontacts = j;
1357
1358
1359 // 3. Leave only unique shortest contacts, that is, if a1[i]-a2[j]
1360 // is the shortest contact for atom a1[i], it has also to be
1361 // the shortest contact for atom a2[j].
1362
1363 if (ncontacts>0) {
1364
1365 SortContacts ( contact,ncontacts,mmdb::CNSORT_DINC );
1366
1367 // get memory for contact expansion shifts
1368 mmdb::GetVectorMemory ( ip,ncontacts,0 );
1369 mmdb::GetVectorMemory ( im,ncontacts,0 );
1370
1371 j = 0;
1372 for (i=0;i<ncontacts;i++) {
1373 i1 = contact[i].id1;
1374 i2 = contact[i].id2;
1375 if ((a1[i1].c<0) && (a2[i2].c<0) && (!isMC(i1,i2))) {
1376 // check that chains contact at sufficiently small angle
1377 if ((i1>0) && (i2>0)) {
1378 m1 = i1-1;
1379 m2 = i2-1;
1380 } else {
1381 m1 = i1;
1382 m2 = i2;
1383 }
1384 if ((i1<nres1-1) && (i2<nres2-1)) {
1385 n1 = i1+1;
1386 n2 = i2+1;
1387 } else {
1388 n1 = i1;
1389 n2 = i2;
1390 }
1391 if (m1!=n1) {
1392 v1[0] = Calpha1[n1]->x - Calpha1[m1]->x;
1393 v1[1] = Calpha1[n1]->y - Calpha1[m1]->y;
1394 v1[2] = Calpha1[n1]->z - Calpha1[m1]->z;
1395 v2[0] = Calpha2[n2]->x - Calpha2[m2]->x;
1396 v2[1] = Calpha2[n2]->y - Calpha2[m2]->y;
1397 v2[2] = Calpha2[n2]->z - Calpha2[m2]->z;
1398 l1 = v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2];
1399 l2 = v2[0]*v2[0] + v2[1]*v2[1] + v2[2]*v2[2];
1400 if ((l1>mmdb::MachEps) && (l2>mmdb::MachEps)) {
1401 cosine = (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])/sqrt(l1*l2);
1402 if (cosine>=minCosine) {
1403 // The contact is valid. Because the contacts are sorted
1404 // by decreasing the contact distance, we simply accept
1405 // all firstly encountered contacts and by this
1406 // invalidate all other contacts of the same atoms
1407 if (j<i) contact[j].Copy ( contact[i] );
1408 // close contact
1409 a1[i1].c = i2;
1410 a2[i2].c = i1;
1411 a1[i1].dist = contact[j].dist;
1412 a2[i2].dist = contact[j].dist;
1413 ip[j] = 1;
1414 im[j] = 1;
1415 j++;
1416 }
1417 }
1418 }
1419 }
1420 }
1421 ncontacts = j;
1422
1423 for (i=0;i<ncontacts-1;i++) {
1424 for (j=0;j<i;j++)
1425 if ((ip[j]>=0) || (im[j]>0))
1426 ExpandContact ( contact[j],ip[j],im[j],contact[i].dist );
1427 ExpandContact ( contact[i],ip[i],im[i],contact[i+1].dist );
1428 }
1429
1430 for (i=0;i<ncontacts;i++)
1431 if (((ip[i]>=0) || (im[i]>0)) && (contact[i].dist>=0.0))
1432 ExpandContact ( contact[i],ip[i],im[i],rmsd_est*rmsd_est );
1433
1434 mmdb::FreeVectorMemory ( ip,0 );
1435 mmdb::FreeVectorMemory ( im,0 );
1436
1437 }
1438
1439 if (contact) delete[] contact;
1440
1441 }
1442
1443
RecoverGaps(mmdb::PPAtom Ca1,PSpAtom at1,int nat1,mmdb::PPAtom Ca2,PSpAtom at2,int nat2,mmdb::realtype thresh)1444 void ssm::Superpose::RecoverGaps (
1445 mmdb::PPAtom Ca1, PSpAtom at1, int nat1,
1446 mmdb::PPAtom Ca2, PSpAtom at2, int nat2,
1447 mmdb::realtype thresh ) {
1448 // Checks if additional mappings may be made when moving
1449 // from sides of gaps to their centers.
1450 mmdb::realtype thr2,d,d2;
1451 int i, k1,k2,n1,n2,sw;
1452 bool B1,B2;
1453
1454 d = 0.0; // to supress warnings
1455 d2 = 0.0; // to supress warnings
1456
1457 thr2 = thresh*thresh;
1458 i = 0;
1459 while (i<nat1) {
1460
1461 // skip all closed contacts
1462 while (i<nat1)
1463 if (at1[i].c>=0) i++;
1464 else break;
1465
1466 if (i<nat1) {
1467
1468 k1 = i; // begining of a gap in 1st structure
1469 if (i>0) n1 = at1[i-1].c+1; // begining of the gap in
1470 else n1 = -1; // 2nd structure
1471
1472 // skip all open contacts
1473 while (i<nat1)
1474 if (at1[i].c<0) i++;
1475 else break;
1476 k2 = i-1; // end of the gap in 1st structure
1477 if (i<nat1) n2 = at1[i].c-1; // end of the gap in 2nd structure
1478 else n2 = -1;
1479
1480 // try to cover the gap
1481 if ((n1<0) && (n2>=0)) {
1482 // The gap is in the begining of 1st structure.
1483 // Recover it from the end to the begining.
1484 while ((k2>=0) && (n2>=0))
1485 if ((at2[n2].c<0) && at2[n2].CompatibleSSE(at1[k2])) {
1486 d = Ca1[k2]->GetDist2 ( Ca2[n2] );
1487 if (d<=thr2) {
1488 at1[k2].c = n2;
1489 at2[n2].c = k2;
1490 at1[k2].dist = d;
1491 at2[n2].dist = d;
1492 n2--;
1493 k2--;
1494 } else
1495 break;
1496 } else
1497 break;
1498
1499 } else if ((n1>=0) && (n2<0)) {
1500 // The gap is in the end of 1st structure.
1501 // Recover it from begining to the end.
1502 while ((k1<nat1) && (n1<nat2))
1503 if (at2[n1].c<0 && at2[n1].CompatibleSSE(at1[k1])) {
1504 d = Ca1[k1]->GetDist2 ( Ca2[n1] );
1505 if (d<=thr2) {
1506 at1[k1].c = n1;
1507 at2[n1].c = k1;
1508 at1[k1].dist = d;
1509 at2[n1].dist = d;
1510 n1++;
1511 k1++;
1512 } else
1513 break;
1514 } else
1515 break;
1516
1517 } else if ((n1>=0) && (n2>=0)) {
1518 // The gap is surrounded by closed contacts in the 1st
1519 // structure. Recover it gradually from begining and the end.
1520 B1 = true;
1521 B2 = true;
1522 while ((k1<=k2) && (n1<=n2) && (B1 || B2)) {
1523 if (B1) {
1524 if ((at2[n1].c>=0) ||
1525 (!at2[n1].CompatibleSSE(at1[k1]))) B1 = false;
1526 else {
1527 d = Ca1[k1]->GetDist2 ( Ca2[n1] );
1528 if (d>thr2) B1 = false;
1529 }
1530 }
1531 if (B2) {
1532 if ((at2[n2].c>=0) ||
1533 (!at2[n2].CompatibleSSE(at1[k2]))) B2 = false;
1534 else {
1535 d2 = Ca1[k2]->GetDist2 ( Ca2[n2] );
1536 if (d2>thr2) B2 = false;
1537 }
1538 }
1539 if (B1 && B2) {
1540 if (d>=d2) sw = 2;
1541 else sw = 1;
1542 } else if (B1) {
1543 sw = 1;
1544 } else if (B2) {
1545 sw = 2;
1546 } else
1547 sw = 0;
1548 if (sw==1) {
1549 at1[k1].c = n1;
1550 at2[n1].c = k1;
1551 at1[k1].dist = d;
1552 at2[n1].dist = d;
1553 n1++;
1554 k1++;
1555 } else if (sw==2) {
1556 at1[k2].c = n2;
1557 at2[n2].c = k2;
1558 at1[k2].dist = d2;
1559 at2[n2].dist = d2;
1560 n2--;
1561 k2--;
1562 }
1563
1564 }
1565
1566 }
1567
1568 }
1569
1570 }
1571
1572 }
1573
1574
1575
CleanShortSections(PSpAtom at1,int nat1,PSpAtom at2)1576 void ssm::Superpose::CleanShortSections ( PSpAtom at1, int nat1,
1577 PSpAtom at2 ) {
1578 // Checks if continuous mapped sections of sequences c1 and c2
1579 // follow the same direction (along protein section) and are
1580 // long enough to be significant. If sections are too short
1581 // then they are unmapped.
1582 int i,j,j2,k,shortSect;
1583
1584 if (nmisdr<2) shortSect = shortSect1;
1585 else shortSect = shortSect2;
1586 nmisdr = 0;
1587
1588 k = -1;
1589 for (i=0;i<nat1;i++)
1590 if (at1[i].c>=0) {
1591 if (k<0) k = i; // begining of a continuous section
1592 else if (at1[i].c<=at1[i-1].c) {
1593 // a misdirection has been found
1594 nmisdr++;
1595 if (i-k<=shortSect) {
1596 // the previous section is too short to be significant;
1597 // waste it:
1598 for (j=k;j<i;j++) {
1599 j2 = at1[j].c;
1600 if ((at1[j ].unmap1!=UNMAP_NO) &&
1601 (at2[j2].unmap1!=UNMAP_NO)) {
1602 at2[j2].c = -1;
1603 at1[j ].c = -1;
1604 }
1605 }
1606 }
1607 k = i; // mark begining of the new continuous section
1608 }
1609 } else if (k>=0) {
1610 // end of the continuous section
1611 if (i-k<=shortSect) {
1612 // the section is too short - dispose of it:
1613 for (j=k;j<i;j++) {
1614 j2 = at1[j].c;
1615 if ((at1[j ].unmap1!=UNMAP_NO) &&
1616 (at2[j2].unmap1!=UNMAP_NO)) {
1617 at2[j2].c = -1;
1618 at1[j ].c = -1;
1619 }
1620 }
1621 }
1622 k = -1; // mark end of section, or begining of a gap
1623 }
1624
1625 if (k>=0) {
1626 // end of last continuous section
1627 if (nat1-k<=shortSect) {
1628 // the section is too short - dispose of it:
1629 for (j=k;j<nat1;j++) {
1630 j2 = at1[j].c;
1631 if ((at1[j ].unmap1!=UNMAP_NO) &&
1632 (at2[j2].unmap1!=UNMAP_NO)) {
1633 at2[j2].c = -1;
1634 at1[j ].c = -1;
1635 }
1636 }
1637 }
1638 }
1639
1640 }
1641
1642
CalcNGaps(PSpAtom a,int nres,int & Ng,int & Nm)1643 void ssm::Superpose::CalcNGaps ( PSpAtom a, int nres,
1644 int & Ng, int & Nm ) {
1645 int i,k,m;
1646
1647 Ng = 0;
1648 Nm = 0;
1649 m = -1;
1650 k = -1;
1651 for (i=0;i<nres;i++)
1652 if (a[i].c0>=0) {
1653 if (k<0) k = i; // begining of a continuous section
1654 else if (a[i].c0<=a[i-1].c0) {
1655 // a misdirection has been found - end of section
1656 Ng++;
1657 Nm++;
1658 k = i; // mark begining of the new continuous section
1659 }
1660 } else if (k>=0) {
1661 // end of the continuous section
1662 Ng++;
1663 if ((m>=0) && (a[k].c0<=m)) Nm++;
1664 m = a[k].c0;
1665 k = -1; // mark end of section, or begining of a gap
1666 }
1667
1668 if (a[nres-1].c0<0) Ng--;
1669
1670 }
1671
CalcNCombs(PGraph G,PSSEDesc SSED,int nSSEs,PSpAtom a,int nres)1672 mmdb::realtype ssm::Superpose::CalcNCombs ( PGraph G, PSSEDesc SSED,
1673 int nSSEs, PSpAtom a, int nres ) {
1674 UNUSED_ARGUMENT(nres);
1675 mmdb::ivector F;
1676 mmdb::realtype nc;
1677 int i,j,k,n;
1678 bool matched;
1679 //Boolean keepIt;
1680
1681 mmdb::GetVectorMemory ( F,nSSEs,1 );
1682
1683 // find matched SSEs
1684 k = 0;
1685
1686 for (i=0;i<nSSEs;i++)
1687 if (SSED[i].match>0) {
1688 matched = false;
1689 for (j=SSED[i].pos;(j<=SSED[i].pend) && (!matched);j++)
1690 matched = (a[j].c0>=0);
1691 if (matched)
1692 F[++k] = i+1;
1693 }
1694
1695 for (i=1;i<k;i++)
1696 for (j=i+1;j<=k;j++)
1697 if (F[j]<F[i]) {
1698 n = F[i];
1699 F[i] = F[j];
1700 F[j] = n;
1701 }
1702
1703 // get number of combinations
1704 nc = mmdb::RMax ( 1.0,G->CalcCombinations(F,k) );
1705
1706 mmdb::FreeVectorMemory ( F,1 );
1707
1708 return mmdb::RMax(1.0,nc);
1709
1710 }
1711
1712
1713 // -----------------------------------------------------------------
1714
Compare(int i,int j)1715 int ssm::SortDist::Compare ( int i, int j ) {
1716 // sort by decreasing: rd[i+1]<=rd[i]
1717
1718 if ((sd[i].unmap1<=sd[j].index) &&
1719 (sd[j].index<=sd[i].unmap2)) return -1;
1720 if ((sd[j].unmap1<=sd[i].index) &&
1721 (sd[i].index<=sd[j].unmap2)) return +1;
1722
1723 if (sd[j].dist>sd[i].dist) return +1;
1724 if (sd[j].dist<sd[i].dist) return -1;
1725
1726 return 0;
1727
1728 }
1729
Swap(int i,int j)1730 void ssm::SortDist::Swap ( int i, int j ) {
1731 mmdb::realtype r;
1732 int k;
1733 r = sd[i].dist; sd[i].dist = sd[j].dist; sd[j].dist = r;
1734 k = sd[i].index; sd[i].index = sd[j].index; sd[j].index = k;
1735 k = sd[i].unmap1; sd[i].unmap1 = sd[j].unmap1; sd[j].unmap1 = k;
1736 k = sd[i].unmap2; sd[i].unmap2 = sd[j].unmap2; sd[j].unmap2 = k;
1737 }
1738
1739
Sort(PSortDistData sdata,int len)1740 void ssm::SortDist::Sort ( PSortDistData sdata, int len ) {
1741 sd = sdata;
1742 QuickSort::Sort ( sdata,len );
1743 }
1744
1745 // -----------------------------------------------------------------
1746
GetSSESpseCenters(RSSEDesc Q1,RSSEDesc Q2,RSSEDesc T1,RSSEDesc T2,mmdb::realtype & qc1,mmdb::realtype & qc2,mmdb::realtype & tc1,mmdb::realtype & tc2)1747 void ssm::Superpose::GetSSESpseCenters ( RSSEDesc Q1, RSSEDesc Q2,
1748 RSSEDesc T1, RSSEDesc T2,
1749 mmdb::realtype & qc1,
1750 mmdb::realtype & qc2,
1751 mmdb::realtype & tc1,
1752 mmdb::realtype & tc2 ) {
1753 mmdb::realtype d2, dq,dt, quality,quality0;
1754 int iq1,pq1,eq1,lq1, it1,pt1,et1, l1;
1755 int iq2,pq2,eq2,lq2, it2,pt2,et2, l2;
1756 int i,j;
1757
1758 l1 = mmdb::IMax(hx_min_len,sd_min_len)/2;
1759
1760 pq1 = mmdb::IMin(Q1.pos+l1,Q1.pend);
1761 pq2 = mmdb::IMin(Q2.pos+l1,Q2.pend);
1762 pt1 = mmdb::IMin(T1.pos+l1,T1.pend);
1763 pt2 = mmdb::IMin(T2.pos+l1,T2.pend);
1764
1765 eq1 = mmdb::IMax(pq1,Q1.pend-l1);
1766 eq2 = mmdb::IMax(pq2,Q2.pend-l1);
1767 et1 = mmdb::IMax(pt1,T1.pend-l1);
1768 et2 = mmdb::IMax(pt2,T2.pend-l1);
1769
1770 quality0 = 0.0;
1771 qc1 = mmdb::realtype(Q1.pos+Q1.pend)/2.0;
1772 qc2 = mmdb::realtype(Q2.pos+Q2.pend)/2.0;
1773 tc1 = mmdb::realtype(T1.pos+T1.pend)/2.0;
1774 tc2 = mmdb::realtype(T2.pos+T2.pend)/2.0;
1775 if ((Q1.len<=0) || (Q2.len<=0) || (T1.len<=0) || (T2.len<=0))
1776 return;
1777
1778 for (iq1=pq1;iq1<=eq1;iq1++) {
1779 lq1 = mmdb::IMin(iq1-Q1.pos,Q1.pend-iq1);
1780 for (it1=pt1;it1<=et1;it1++) {
1781 l1 = mmdb::IMin ( lq1,mmdb::IMin(it1-T1.pos,T1.pend-it1) );
1782
1783 for (iq2=pq2;iq2<=eq2;iq2++) {
1784 lq2 = mmdb::IMin(iq2-Q2.pos,Q2.pend-iq2);
1785 for (it2=pt2;it2<=et2;it2++) {
1786 l2 = mmdb::IMin ( lq2,mmdb::IMin(it2-T2.pos,T2.pend-it2) );
1787
1788 d2 = 0.0;
1789 for (i=-l1;i<=l1;i++)
1790 for (j=-l2;j<=l2;j++) {
1791 dq = Calpha1[iq1+i]->GetDist2 ( Calpha1[iq2+j] );
1792 dt = Calpha2[it1+i]->GetDist2 ( Calpha2[it2+j] );
1793 d2 += dq + dt - 2.0*sqrt(dq*dt);
1794 }
1795
1796 dq = (2*l1+1)*(2*l2+1);
1797 quality = sqrt(dq)/(1.0+d2/(dq*Rmsd0*Rmsd0));
1798 if (quality>quality0) {
1799 quality0 = quality;
1800 qc1 = iq1;
1801 qc2 = iq2;
1802 tc1 = it1;
1803 tc2 = it2;
1804 }
1805
1806 }
1807 }
1808
1809 }
1810 }
1811
1812 }
1813
1814
FirstGuess(mmdb::ivector F1,mmdb::ivector F2,int mlen)1815 int ssm::Superpose::FirstGuess ( mmdb::ivector F1, mmdb::ivector F2,
1816 int mlen ) {
1817 mmdb::rvector c1,c2;
1818 mmdb::realtype xc1,yc1,zc1, xc2,yc2,zc2, r1,r2;
1819 mmdb::vect3 vc1,vc2;
1820 int i,j, i1,i2, j1,j2, l, na;
1821
1822
1823 // 1. Mark all atoms as non-corresponding
1824
1825 for (i=0;i<nres1;i++) a1[i].c = -1;
1826 for (i=0;i<nres2;i++) a2[i].c = -1;
1827
1828
1829 // 2. Correspond C-alphas of matched SSEs and calculate
1830 // the mass centers
1831
1832 mmdb::GetVectorMemory ( c1,mlen,1 );
1833 mmdb::GetVectorMemory ( c2,mlen,1 );
1834 if (mlen>1) {
1835 for (i=1;i<=mlen;i++) {
1836 c1[i] = 0.0;
1837 c2[i] = 0.0;
1838 }
1839 l = 0;
1840 for (i=1;i<mlen;i++)
1841 for (j=i+1;j<=mlen;j++) {
1842 GetSSESpseCenters ( SSED1[F1[i]-1],SSED1[F1[j]-1],
1843 SSED2[F2[i]-1],SSED2[F2[j]-1],
1844 xc1,yc1, xc2,yc2 );
1845 if ((xc1>=0.0) && (yc1>=0.0) && (xc2>=0.0) && (yc2>=0.0)) {
1846 c1[i] += xc1;
1847 c1[j] += yc1;
1848 c2[i] += xc2;
1849 c2[j] += yc2;
1850 } else
1851 l++;
1852 }
1853 l = mlen - 1 - l;
1854 for (i=1;i<=mlen;i++) {
1855 c1[i] /= l;
1856 c2[i] /= l;
1857 }
1858 } else {
1859 i = F1[1] - 1;
1860 c1[1] = mmdb::realtype(SSED1[i].pos+SSED1[i].pend)/2.0;
1861 i = F2[1] - 1;
1862 c2[1] = mmdb::realtype(SSED2[i].pos+SSED2[i].pend)/2.0;
1863 }
1864
1865 na = 0;
1866 xc1 = 0.0;
1867 yc1 = 0.0;
1868 zc1 = 0.0;
1869 xc2 = 0.0;
1870 yc2 = 0.0;
1871 zc2 = 0.0;
1872 for (i=1;i<=mlen;i++) {
1873 i1 = F1[i] - 1;
1874 i2 = F2[i] - 1;
1875 SSED1[i1].m = F2[i];
1876 SSED2[i2].m = F1[i];
1877 if ((SSED1[i1].len>0) && (SSED2[i2].len>0)) {
1878 r1 = mmdb::RMin ( c1[i]-SSED1[i1].pos,SSED1[i1].pend-c1[i] );
1879 r2 = mmdb::RMin ( c2[i]-SSED2[i2].pos,SSED2[i2].pend-c2[i] );
1880 r1 = mmdb::RMin ( r1,r2 );
1881 j1 = mmdb::mround ( c1[i]-r1 );
1882 j2 = mmdb::mround ( c2[i]-r1 );
1883 l = mmdb::mround ( 2.0*r1 );
1884 for (j=0;j<=l;j++) {
1885 a1[j1].c = j2;
1886 a2[j2].c = j1;
1887 xc1 += Calpha1[j1]->x;
1888 yc1 += Calpha1[j1]->y;
1889 zc1 += Calpha1[j1]->z;
1890 xc2 += Calpha2[j2]->x;
1891 yc2 += Calpha2[j2]->y;
1892 zc2 += Calpha2[j2]->z;
1893 j1++;
1894 j2++;
1895 na++;
1896 }
1897 }
1898 }
1899
1900 mmdb::FreeVectorMemory ( c1,1 );
1901 mmdb::FreeVectorMemory ( c2,1 );
1902
1903 if (na<=0) return 1;
1904
1905 xc1 /= na;
1906 yc1 /= na;
1907 zc1 /= na;
1908 xc2 /= na;
1909 yc2 /= na;
1910 zc2 /= na;
1911
1912 // 3. Calculate the correlation matrix
1913
1914 for (i=1;i<=3;i++)
1915 for (j=1;j<=3;j++)
1916 A[i][j] = 0.0;
1917
1918 for (i1=0;i1<nres1;i1++) {
1919 i2 = a1[i1].c;
1920 if (i2>=0) {
1921 vc1[0] = Calpha1[i1]->x - xc1;
1922 vc1[1] = Calpha1[i1]->y - yc1;
1923 vc1[2] = Calpha1[i1]->z - zc1;
1924 vc2[0] = Calpha2[i2]->x - xc2;
1925 vc2[1] = Calpha2[i2]->y - yc2;
1926 vc2[2] = Calpha2[i2]->z - zc2;
1927 for (i=1;i<=3;i++)
1928 for (j=1;j<=3;j++)
1929 A[i][j] += vc1[j-1]*vc2[i-1];
1930 }
1931 }
1932
1933
1934 // 4. Calculate the rotation matrix and add translation
1935
1936 if (CalculateTMatrix()) return 2;
1937
1938 TMatrix[0][3] = xc2 - TMatrix[0][0]*xc1 -
1939 TMatrix[0][1]*yc1 - TMatrix[0][2]*zc1;
1940 TMatrix[1][3] = yc2 - TMatrix[1][0]*xc1 -
1941 TMatrix[1][1]*yc1 - TMatrix[1][2]*zc1;
1942 TMatrix[2][3] = zc2 - TMatrix[2][0]*xc1 -
1943 TMatrix[2][1]*yc1 - TMatrix[2][2]*zc1;
1944
1945 return 0;
1946
1947 }
1948
1949
1950
ChooseFirstRotation(int rotSSE1,int rotSSE2)1951 void ssm::Superpose::ChooseFirstRotation ( int rotSSE1, int rotSSE2 ) {
1952 // In case when only one SSE is matched, the initial rotation
1953 // matrix is defined only up to rotation about that SSE.
1954 // Given that matched SSE from 1st structure, this procedure
1955 // tries to choose a more appropriate rotation such that other
1956 // SSEs are most closely approached by each other.
1957 mmdb::mat44 m,tm,mm;
1958 mmdb::PContact contact;
1959 mmdb::realtype alpha, vx,vy,vz, x0,y0,z0, B, score,score0;
1960 int ia1,ia2,nrot,i,j,k,l,ncontacts,na,na0, i1,i2;
1961
1962 ia1 = rotSSE1 - 1;
1963 ia2 = rotSSE2 - 1;
1964 ncontacts = 0;
1965 for (j=0;j<nSSEs1;j++)
1966 if (j!=ia1)
1967 for (k=0;k<nSSEs2;k++)
1968 if ((k!=ia2) && (SSED1[j].type==SSED2[k].type))
1969 ncontacts++;
1970
1971 if (ncontacts>0) {
1972
1973 contact = new mmdb::Contact[nSSEs1*nSSEs2];
1974
1975 nrot = 36; // number of rotations to try
1976
1977 x0 = SSED1[ia1].xs1;
1978 y0 = SSED1[ia1].ys1;
1979 z0 = SSED1[ia1].zs1;
1980 vx = SSED1[ia1].xs2 - x0;
1981 vy = SSED1[ia1].ys2 - y0;
1982 vz = SSED1[ia1].zs2 - z0;
1983
1984 na0 = -1;
1985 score0 = mmdb::MaxReal;
1986
1987 for (i=0;i<nrot;i++) {
1988
1989 // add rotation to the transformation matrix
1990 alpha = i*2.0*mmdb::Pi/nrot;
1991 mmdb::GetVecTMatrix ( m,alpha,vx,vy,vz,x0,y0,z0 );
1992 for (j=0;j<4;j++)
1993 for (k=0;k<4;k++) {
1994 B = 0.0;
1995 for (l=0;l<4;l++)
1996 B += m[j][l]*TMatrix[l][k];
1997 tm[j][k] = B;
1998 }
1999
2000 // find contacts between SSEs
2001 ncontacts = 0;
2002 for (j=0;j<nSSEs1;j++)
2003 if (j!=ia1) {
2004 SSED1[j].Transform ( tm );
2005 for (k=0;k<nSSEs2;k++)
2006 if ((k!=ia2) && (SSED1[j].type==SSED2[k].type) &&
2007 (SSED1[j].Cosine(SSED2[k])>=minCosine)) {
2008 SSED1[j].CalcScore ( SSED2[k] );
2009 if (SSED1[j].score<=maxContact/4.0) {
2010 contact[ncontacts].id1 = j;
2011 contact[ncontacts].id2 = k;
2012 contact[ncontacts].dist = SSED1[j].score;
2013 ncontacts++;
2014 }
2015 }
2016 }
2017
2018 // count potential atom contacts in mapped SSEs
2019
2020 for (j=0;j<nSSEs1;j++) SSED1[j].m = -1;
2021 for (j=0;j<nSSEs2;j++) SSED2[j].m = -1;
2022 SSED1[ia1].m = ia2;
2023 SSED2[ia2].m = ia1;
2024
2025 SortContacts ( contact,ncontacts,mmdb::CNSORT_DINC );
2026 na = 0;
2027 score = 0.0;
2028 for (j=0;j<ncontacts;j++) {
2029 i1 = contact[j].id1;
2030 i2 = contact[j].id2;
2031 if ((SSED1[i1].m<0) && (SSED2[i2].m<0)) {
2032 l = 0;
2033 if (!allowMC) {
2034 // check for misfolding
2035 k = i1 + 1;
2036 while (k<nSSEs1)
2037 if (SSED1[k].m>=0) break;
2038 else k++;
2039 if (k<nSSEs1) {
2040 if (i2>=SSED1[k].m) l = 1;
2041 }
2042 k = i1 - 1;
2043 while (k>=0)
2044 if (SSED1[k].m>=0) break;
2045 else k--;
2046 if (k>=0) {
2047 if (i2<=SSED1[k].m) l = 1;
2048 }
2049 }
2050 if (!l) {
2051 // count contact
2052 na += mmdb::IMin(SSED1[i1].len,SSED2[i2].len);
2053 // bookkeep the mapping; this invalidate other contacts
2054 // of mapped SSEs
2055 SSED1[i1].m = i2;
2056 SSED2[i2].m = i1;
2057 score += contact[j].dist;
2058 }
2059 }
2060 }
2061
2062 // keep rotation with maximal number of potential mappings
2063 if ((na>na0) || ((na==na0) && (score<score0))) {
2064 na0 = na;
2065 score0 = score;
2066 for (j=0;j<4;j++)
2067 for (k=0;k<4;k++)
2068 mm[j][k] = tm[j][k];
2069 }
2070
2071 }
2072
2073 delete[] contact;
2074
2075 if (na0>0)
2076 for (j=0;j<4;j++)
2077 for (k=0;k<4;k++)
2078 TMatrix[j][k] = mm[j][k];
2079
2080 }
2081
2082 }
2083
2084
MatchQuality(int Nalign,mmdb::realtype Rmsd)2085 mmdb::realtype ssm::Superpose::MatchQuality ( int Nalign, mmdb::realtype Rmsd ) {
2086 if (Nalign==0) return 0.0;
2087 return MatchQuality2 ( Nalign,Rmsd*Rmsd*Nalign );
2088 }
2089
MatchQuality2(int Nalign,mmdb::realtype dist2)2090 mmdb::realtype ssm::Superpose::MatchQuality2 ( int Nalign, mmdb::realtype dist2 ) {
2091 mmdb::realtype NormN,Na2,NormR;
2092
2093 if (Nalign<=0) return 0.0;
2094
2095 NormN = nres1*nres2;
2096 if (NormN<=0.0) return 0.0;
2097 Na2 = Nalign*Nalign;
2098 NormR = dist2/(Nalign*Rmsd0*Rmsd0);
2099
2100 return Na2/((1.0+NormR)*NormN);
2101
2102 }
2103
2104
UnmapExcluded(PSpAtom a1,PSpAtom a2,int nres1)2105 void ssm::Superpose::UnmapExcluded ( PSpAtom a1, PSpAtom a2,
2106 int nres1 ) {
2107 int i;
2108 for (i=0;i<nres1;i++)
2109 if (a1[i].excluded && (a1[i].c>=0)) {
2110 a2[a1[i].c].c = -1;
2111 a1[i].c = -1;
2112 }
2113 }
2114
2115
OptimizeNalign()2116 int ssm::Superpose::OptimizeNalign() {
2117 // Finds maximal Nalign such that RMSD does not exceed the given Rmsd
2118 SectionDist D;
2119 mmdb::vect3 vc1,vc2;
2120 mmdb::realtype xc1,yc1,zc1, xc2,yc2,zc2, B;
2121 mmdb::realtype dist2,Q,Q1,maxRMSD2, contDist;
2122 int na,nl;
2123 int i,j,i1,i2,iter,nobetter_cnt,rc;
2124 bool Done;
2125
2126 // Make iterative atom-to-atom superposition
2127
2128 Q_achieved = -1.0;
2129 iter = 0;
2130 nobetter_cnt = 0;
2131 rc = SPOSE_Ok;
2132 maxRMSD2 = maxRMSD*maxRMSD;
2133
2134 do {
2135
2136 B = mmdb::IMax(iter,iterMin);
2137 B = iter/B;
2138 contDist = minContact + (maxContact-minContact)*B;
2139
2140 // 1. Find corresponding atoms
2141
2142 // 1.1 Transform 1st structure to the current superposition
2143 for (i=0;i<nres1;i++)
2144 Calpha1[i]->Transform ( TMatrix );
2145
2146 for (i=0;i<nSSEs1;i++)
2147 SSED1[i].Transform ( TMatrix );
2148
2149 // 1.2 Mark all atoms as non-corresponding
2150 for (i=0;i<nres1;i++) {
2151 a1[i].c = -1;
2152 a1[i].unmap1 = UNMAP_YES;
2153 a1[i].unmap2 = UNMAP_YES;
2154 }
2155 for (i=0;i<nres2;i++) {
2156 a2[i].c = -1;
2157 a2[i].unmap1 = UNMAP_YES;
2158 a2[i].unmap2 = UNMAP_YES;
2159 }
2160 for (i=0;i<nSSEs1;i++) SSED1[i].m = -1;
2161 for (i=0;i<nSSEs2;i++) SSED2[i].m = -1;
2162
2163 // 1.3 Correspond atoms of matched SSEs first
2164 for (i=1;i<=copyFlen;i++) {
2165 CalcDistance ( copyF1[i],copyF2[i],D );
2166 AlignSSEs ( D,UNMAP_NO );
2167 }
2168
2169 // 1.4 Try to correspond other SSEs
2170 CorrespondSSEs ( FH1,nFH1,FH2,nFH2,contDist );
2171 CorrespondSSEs ( FS1,nFS1,FS2,nFS2,contDist );
2172
2173 // 1.5 Correspond all other atoms
2174
2175 CorrespondContacts ( MMDB1,contDist );
2176
2177 RecoverGaps ( Calpha1,a1,nres1,Calpha2,a2,nres2,2.0*contDist );
2178 RecoverGaps ( Calpha2,a2,nres2,Calpha1,a1,nres1,2.0*contDist );
2179
2180 if (selInclHnd1>0) UnmapExcluded ( a1,a2,nres1 );
2181 if (selInclHnd2>0) UnmapExcluded ( a2,a1,nres2 );
2182
2183 // 1.6 Find atoms, which may be unmapped for optimizing
2184 // the quality function
2185 dist2 = 0.0; // square distance between the structures at
2186 // current rotation
2187 nl = 0; // number of atoms that may be unmapped
2188 na = 0; // total number of corresponding atoms
2189 for (i1=0;i1<nres1;i1++) {
2190 i2 = a1[i1].c;
2191 if (i2>=0) {
2192 na++;
2193 dist2 += a1[i1].dist;
2194 if ((a1[i1].unmap1!=UNMAP_NO) && (a2[i2].unmap1!=UNMAP_NO)) {
2195 sdata[nl].dist = a1[i1].dist;
2196 sdata[nl].index = i1;
2197 sdata[nl].unmap1 = a1[i1].unmap1;
2198 sdata[nl].unmap2 = a1[i1].unmap2;
2199 nl++;
2200 }
2201 }
2202 }
2203
2204 // 1.7 Unmap atoms for optimizing the quality function
2205 if ((nl>0) && (na>3)) {
2206 sortDist.Sort ( sdata,nl );
2207 i = 0;
2208 if (dist2<=na*maxRMSD2) {
2209 Q1 = MatchQuality2 ( na,dist2 );
2210 j = 0;
2211 } else {
2212 Q1 = -1.0;
2213 j = nl;
2214 }
2215 while ((i<nl) && (na>3)) {
2216 dist2 -= sdata[i].dist;
2217 na--;
2218 i++;
2219 if (dist2<=na*maxRMSD2) { // rmsd must be within the limits
2220 Q = MatchQuality2 ( na,dist2 );
2221 if (Q>Q1) {
2222 Q1 = Q;
2223 j = i;
2224 }
2225 }
2226 }
2227 for (i=0;i<j;i++) {
2228 i1 = sdata[i].index;
2229 i2 = a1[i1].c;
2230 a1[i1].c = -1;
2231 a2[i2].c = -1;
2232 }
2233 }
2234
2235 // 1.8 Clean up short sections
2236 CleanShortSections ( a1,nres1,a2 );
2237 CleanShortSections ( a2,nres2,a1 );
2238
2239
2240 // 2. Make rotation with the new correspondence
2241
2242 // 2.1 Calculate centers of mass and rmsd
2243 xc1 = 0.0; // mass
2244 yc1 = 0.0; // center of
2245 zc1 = 0.0; // 1st structure
2246 xc2 = 0.0; // mass
2247 yc2 = 0.0; // center of
2248 zc2 = 0.0; // 1st structure
2249 na = 0; // total number of corresponding atoms
2250 dist2 = 0.0; // square distance between the structures
2251 // at current rotation
2252 for (i1=0;i1<nres1;i1++) {
2253 i2 = a1[i1].c;
2254 if (i2>=0) {
2255 xc1 += cax0[i1];
2256 yc1 += cay0[i1];
2257 zc1 += caz0[i1];
2258 xc2 += Calpha2[i2]->x;
2259 yc2 += Calpha2[i2]->y;
2260 zc2 += Calpha2[i2]->z;
2261 dist2 += a1[i1].dist;
2262 na++;
2263 }
2264 }
2265
2266 // 2.2 Reset rotating coordinates for the next iteration
2267 for (i=0;i<nres1;i++) {
2268 Calpha1[i]->x = cax0[i];
2269 Calpha1[i]->y = cay0[i];
2270 Calpha1[i]->z = caz0[i];
2271 }
2272
2273 Q = MatchQuality2 ( na,dist2 );
2274
2275 nobetter_cnt++;
2276 if (Q>Q_achieved) {
2277 // 2.3 Retain the current approximation
2278 if (Q-Q_achieved>minQStep) nobetter_cnt = 0;
2279 Q_achieved = Q;
2280 for (i=0;i<4;i++)
2281 for (j=0;j<4;j++)
2282 TMx[i][j] = TMatrix[i][j];
2283 for (i=0;i<nres1;i++) {
2284 a1[i].c0 = a1[i].c;
2285 a1[i].dist0 = a1[i].dist;
2286 }
2287 for (i=0;i<nres2;i++) {
2288 a2[i].c0 = a2[i].c;
2289 a2[i].dist0 = a2[i].dist;
2290 }
2291 for (i=0;i<nSSEs1;i++) SSED1[i].match = SSED1[i].m;
2292 for (i=0;i<nSSEs2;i++) SSED2[i].match = SSED2[i].m;
2293 if (na>0) rmsd_achieved = sqrt(dist2/na);
2294 else rmsd_achieved = 0.0;
2295 nalgn = na;
2296 }
2297
2298 Done = (na<=0) || (iter>iterMax) ||
2299 ((iter>iterMin) && (nobetter_cnt>maxHollowIt));
2300
2301 if (!Done) {
2302
2303 // 2.4 Calculate the correlation matrix for the next iteration
2304 xc1 /= na;
2305 yc1 /= na;
2306 zc1 /= na;
2307 xc2 /= na;
2308 yc2 /= na;
2309 zc2 /= na;
2310
2311 for (i=1;i<=3;i++)
2312 for (j=1;j<=3;j++)
2313 A[i][j] = 0.0;
2314
2315 for (i1=0;i1<nres1;i1++) {
2316 i2 = a1[i1].c;
2317 if (i2>=0) {
2318 vc1[0] = Calpha1[i1]->x - xc1;
2319 vc1[1] = Calpha1[i1]->y - yc1;
2320 vc1[2] = Calpha1[i1]->z - zc1;
2321 vc2[0] = Calpha2[i2]->x - xc2;
2322 vc2[1] = Calpha2[i2]->y - yc2;
2323 vc2[2] = Calpha2[i2]->z - zc2;
2324 B = 1.0;
2325 for (i=1;i<=3;i++)
2326 for (j=1;j<=3;j++)
2327 A[i][j] += B*vc1[j-1]*vc2[i-1];
2328 }
2329 }
2330
2331 if (CalculateTMatrix()) rc = SPOSE_SVDFail;
2332 else {
2333
2334 // 5.3.e add translation
2335 TMatrix[0][3] = xc2 - TMatrix[0][0]*xc1 -
2336 TMatrix[0][1]*yc1 - TMatrix[0][2]*zc1;
2337 TMatrix[1][3] = yc2 - TMatrix[1][0]*xc1 -
2338 TMatrix[1][1]*yc1 - TMatrix[1][2]*zc1;
2339 TMatrix[2][3] = zc2 - TMatrix[2][0]*xc1 -
2340 TMatrix[2][1]*yc1 - TMatrix[2][2]*zc1;
2341
2342 }
2343
2344 }
2345
2346 iter++;
2347
2348 } while ((rc==SPOSE_Ok) && (!Done));
2349
2350 if (na<=0) rc = SPOSE_RemoteStruct;
2351 else
2352 for (i=0;i<4;i++)
2353 for (j=0;j<4;j++)
2354 TMatrix[i][j] = TMx[i][j];
2355
2356 return rc;
2357
2358 }
2359
2360
CalcQScore(RSSEDesc SSE1)2361 void ssm::Superpose::CalcQScore ( RSSEDesc SSE1 ) {
2362 // Calculates Q-score of overlaping SSE1 with matched SSE
2363 mmdb::realtype NormN,dist2;
2364 int p1,p2,na,i;
2365
2366 SSE1.Qscore = 0.0;
2367 SSE1.Rscore = 0.0;
2368 SSE1.Xscore = 0.0;
2369 if (SSE1.match<=0) return;
2370
2371 SSED2[SSE1.match-1].Qscore = 0.0;
2372 SSED2[SSE1.match-1].Rscore = 0.0;
2373 SSED2[SSE1.match-1].Xscore = 0.0;
2374 NormN = SSE1.len*SSED2[SSE1.match-1].len;
2375 if (NormN<=0.0) return;
2376
2377 p1 = SSED2[SSE1.match-1].pos;
2378 p2 = SSED2[SSE1.match-1].pend;
2379 if ((SSE1.pos>=0) && (SSE1.pend>=SSE1.pos) &&
2380 (p1>=0) && (p2>=p1)) {
2381 na = 0;
2382 dist2 = 0.0;
2383 for (i=SSE1.pos;i<=mmdb::IMin(nres1-1,SSE1.pend);i++)
2384 if ((p1<=a1[i].c0) && (a1[i].c0<=p2)) {
2385 dist2 += a1[i].dist0;
2386 na++;
2387 }
2388 if (na>0) {
2389 dist2 /= na*Rmsd0*Rmsd0;
2390 SSE1.Rscore = 1.0/(1.0+dist2);
2391 SSE1.Xscore = mmdb::Exp ( -dist2 );
2392 SSE1.Qscore = na*na/((1.0+dist2)*NormN);
2393 } else {
2394 SSE1.Rscore = 0.0;
2395 SSE1.Xscore = 0.0;
2396 SSE1.Qscore = 0.0;
2397 }
2398 }
2399
2400 SSED2[SSE1.match-1].Qscore = SSE1.Qscore;
2401 SSED2[SSE1.match-1].Rscore = SSE1.Rscore;
2402 SSED2[SSE1.match-1].Xscore = SSE1.Xscore;
2403
2404 }
2405
2406
SuperposeCalphas(PGraph G1,PGraph G2,mmdb::ivector F1,mmdb::ivector F2,int mlen,mmdb::PManager M1,mmdb::PManager M2,int selHndIncl1,int selHndIncl2)2407 int ssm::Superpose::SuperposeCalphas (
2408 PGraph G1, // SSE graph of 1st structure
2409 PGraph G2, // SSE graph of 2nd structure
2410 mmdb::ivector F1, // matched vertices of G1 [1..mlen]
2411 mmdb::ivector F2, // matched vertices of G2 [1..mlen]
2412 int mlen, // length of match (F1,F2)
2413 mmdb::PManager M1, // 1st structure
2414 mmdb::PManager M2, // 2nd structure
2415 int selHndIncl1, // sel handle to include atoms from M1
2416 int selHndIncl2 // sel handle to include atoms from M2
2417 ) {
2418 int i,j, rc;
2419 int selHnd1,selHnd2;
2420
2421 driverID = 1;
2422
2423 rc = 0; // termination on errors
2424
2425 rmsd_achieved = 0.0;
2426 nres1 = 0;
2427 nres2 = 0;
2428 nalgn = 0;
2429 ngaps = 0;
2430 for (i=0;i<4;i++) {
2431 for (j=0;j<4;j++)
2432 TMatrix[i][j] = 0.0;
2433 TMatrix[i][i] = 1.0;
2434 }
2435
2436 FreeMemory();
2437
2438 if ((!G1) || (!G2) ||
2439 (!F1) || (!F2) || (mlen<=0) ||
2440 (!M1) || (!M2)) return SPOSE_BadData;
2441
2442 copyF1 = F1;
2443 copyF2 = F2;
2444 copyFlen = mlen;
2445 MMDB1 = M1;
2446 MMDB2 = M2;
2447
2448 // 1. Select Calphas from both structures
2449
2450 selInclHnd1 = selHndIncl1;
2451 selInclHnd2 = selHndIncl2;
2452
2453 // SelectCalphas obtains Calpha1,2 from MMDB manager and
2454 // allocates a1,2
2455 selHnd1 = 0;
2456 selHnd2 = 0;
2457 SelectCalphas ( MMDB1,G1,Calpha1,a1,nres1,selHnd1,
2458 selInclHnd1,selString1 );
2459 SelectCalphas ( MMDB2,G2,Calpha2,a2,nres2,selHnd2,
2460 selInclHnd2,selString2 );
2461 if ((nres1<=0) || (nres2<=0)) {
2462 if (nres1<=0) rc = SPOSE_NoCalphas1;
2463 else rc = SPOSE_NoCalphas2;
2464 nres1 = 0;
2465 nres2 = 0;
2466 if (a1) delete[] a1;
2467 if (a2) delete[] a2;
2468 a1 = NULL;
2469 a2 = NULL;
2470 MMDB1->DeleteSelection ( selHnd1 );
2471 MMDB2->DeleteSelection ( selHnd2 );
2472 return rc;
2473 }
2474
2475 // 2. Map vertices on the selected Calphas
2476
2477 MapSSEs ( Calpha1,a1,nres1,G1,SSED1,nSSEs1 );
2478 MapSSEs ( Calpha2,a2,nres2,G2,SSED2,nSSEs2 );
2479
2480
2481 // 3. Common part for different Superpose drivers
2482
2483 _superpose ( G1,G2,rc );
2484
2485 // 4. Deselect C-alphas
2486
2487 MMDB1->DeleteSelection ( selHnd1 );
2488 MMDB2->DeleteSelection ( selHnd2 );
2489
2490 return rc;
2491
2492 }
2493
2494
2495 // ----------------------- ssm::SuperposeData -----------------------
2496
Init()2497 void ssm::SuperposeData::Init() {
2498 G = NULL;
2499 M = NULL;
2500 a = NULL;
2501 Calpha = NULL;
2502 SSED = NULL;
2503 selstring = NULL;
2504 selHnd = 0;
2505 selHndIncl = 0;
2506 nres = 0;
2507 nSSEs = 0;
2508 }
2509
Dispose()2510 void ssm::SuperposeData::Dispose() {
2511 // This DOES NOT dispose the graph and MMDB instance, however
2512 // selected C-alphas are deselected
2513
2514 if (a) delete[] a;
2515 if (SSED) delete[] SSED;
2516 if (selstring) delete[] selstring;
2517
2518 a = NULL;
2519 nres = 0;
2520 SSED = NULL;
2521 selstring = NULL;
2522 nSSEs = 0;
2523
2524 DeselectCalphas();
2525
2526 }
2527
DeselectCalphas()2528 void ssm::SuperposeData::DeselectCalphas() {
2529 if (M && (selHnd>0)) M->DeleteSelection ( selHnd );
2530 selHnd = 0;
2531 }
2532
SelectCalphas()2533 void ssm::SuperposeData::SelectCalphas() {
2534 if (G && M && (selHnd<=0)) {
2535 G->SelectCalphas ( M,selHnd,selstring );
2536 M->GetSelIndex ( selHnd,Calpha,nres );
2537 }
2538 }
2539
2540 // ----------------------------------------------------------------
2541
SuperposeCalphas(PSuperposeData SD1,PSuperposeData SD2,mmdb::ivector F1,mmdb::ivector F2,int mlen)2542 int ssm::Superpose::SuperposeCalphas (
2543 PSuperposeData SD1, // superposition data of 1st structure
2544 PSuperposeData SD2, // superposition data of 2nd structure
2545 mmdb::ivector F1, // matched vertices of SD1.G [1..mlen]
2546 mmdb::ivector F2, // matched vertices of SD2.G [1..mlen]
2547 int mlen // length of match (F1,F2)
2548 ) {
2549 int i,j, rc;
2550
2551 driverID = 2;
2552
2553 rc = 0; // termination on errors
2554
2555 rmsd_achieved = 0.0;
2556 nres1 = 0;
2557 nres2 = 0;
2558 nalgn = 0;
2559 ngaps = 0;
2560 for (i=0;i<4;i++) {
2561 for (j=0;j<4;j++)
2562 TMatrix[i][j] = 0.0;
2563 TMatrix[i][i] = 1.0;
2564 }
2565
2566 FreeMemory();
2567
2568 if ((!SD1->G) || (!SD2->G) ||
2569 (!SD1->M) || (!SD2->M) ||
2570 (!F1) || (!F2) || (mlen<=0))
2571 return SPOSE_BadData;
2572
2573 copyF1 = F1;
2574 copyF2 = F2;
2575 copyFlen = mlen;
2576 MMDB1 = SD1->M;
2577 MMDB2 = SD2->M;
2578
2579 // 1. Select Calphas from both structures
2580
2581 selInclHnd1 = SD1->selHndIncl;
2582 selInclHnd2 = SD2->selHndIncl;
2583
2584 // SelectCalphas obtains Calpha1,2 from MMDB manager and
2585 // allocates a1,2
2586 SelectCalphas ( MMDB1,SD1->G,SD1->Calpha,SD1->a,SD1->nres,
2587 SD1->selHnd,SD1->selHndIncl,SD1->selstring );
2588 SelectCalphas ( MMDB2,SD2->G,SD2->Calpha,SD2->a,SD2->nres,
2589 SD2->selHnd,SD2->selHndIncl,SD2->selstring );
2590
2591 if ((SD1->nres<=0) || (SD2->nres<=0)) {
2592 if (SD1->nres<=0) rc = SPOSE_NoCalphas1;
2593 else rc = SPOSE_NoCalphas2;
2594 SD1->Dispose();
2595 SD2->Dispose();
2596 return rc;
2597 }
2598
2599 Calpha1 = SD1->Calpha;
2600 Calpha2 = SD2->Calpha;
2601 a1 = SD1->a;
2602 a2 = SD2->a;
2603 nres1 = SD1->nres;
2604 nres2 = SD2->nres;
2605
2606 // 2. Map vertices on the selected Calphas
2607
2608 MapSSEs ( Calpha1,a1,nres1,SD1->G,SD1->SSED,SD1->nSSEs );
2609 MapSSEs ( Calpha2,a2,nres2,SD2->G,SD2->SSED,SD2->nSSEs );
2610
2611 SSED1 = SD1->SSED;
2612 SSED2 = SD2->SSED;
2613 nSSEs1 = SD1->nSSEs;
2614 nSSEs2 = SD2->nSSEs;
2615
2616 // 3. Common part for different Superpose drivers
2617
2618 _superpose ( SD1->G,SD2->G,rc );
2619
2620
2621 // 4. Prevent external data from being lost
2622 // NOTE that "ssm::Superpose::GetSuperposition()" and
2623 // "ssm::Superpose::GetSSEDescX();" do not work with this driver.
2624 // The superposition data are to be read directly from
2625 // returned SD's
2626
2627 a1 = NULL;
2628 a2 = NULL;
2629 SSED1 = NULL;
2630 SSED2 = NULL;
2631 Calpha1 = NULL;
2632 Calpha2 = NULL;
2633 nres1 = 0;
2634 nres2 = 0;
2635
2636 return rc;
2637
2638 }
2639
2640
_superpose(PGraph G1,PGraph G2,int & rc)2641 void ssm::Superpose::_superpose ( PGraph G1, PGraph G2,
2642 int & rc ) {
2643 mmdb::mat44 TMx0;
2644 int i,j,i1,i2,AD_alloc;
2645
2646 IdentifyUnmatchedSSEs ( FH1,nFH1,FS1,nFS1,copyF1,copyFlen,G1 );
2647 IdentifyUnmatchedSSEs ( FH2,nFH2,FS2,nFS2,copyF2,copyFlen,G2 );
2648
2649 // 3. Allocate memory for corresponding atoms and shortest contacts
2650
2651 sdata = new SortDistData[nres1];
2652
2653 mmdb::GetVectorMemory ( cax0,nres1,0 );
2654 mmdb::GetVectorMemory ( cay0,nres1,0 );
2655 mmdb::GetVectorMemory ( caz0,nres1,0 );
2656 for (i=0;i<nres1;i++) {
2657 cax0[i] = Calpha1[i]->x;
2658 cay0[i] = Calpha1[i]->y;
2659 caz0[i] = Calpha1[i]->z;
2660 }
2661
2662 AD_alloc = 0;
2663 for (i=0;i<nSSEs1;i++)
2664 if (SSED1[i].len>AD_alloc) AD_alloc = SSED1[i].len;
2665 for (i=0;i<nSSEs2;i++)
2666 if (SSED2[i].len>AD_alloc) AD_alloc = SSED2[i].len;
2667
2668 // matrix AD is used as a temporary storage at calculating
2669 // an optimal superposition of SSE pairs
2670 mmdb::GetMatrixMemory ( AD,AD_alloc,AD_alloc,0,0 );
2671
2672 SDistAlloc = mmdb::IMax ( nFH1,nFS1 );
2673 j = mmdb::IMax ( nFH2,nFS2 );
2674 if ((SDistAlloc>0) && (j>0)) {
2675 SDist = new PSectionDist[SDistAlloc];
2676 for (i=0;i<SDistAlloc;i++)
2677 SDist[i] = new SectionDist[j];
2678 }
2679
2680
2681 // 4. Superpose secondary structures; TMx will be the initial guess
2682
2683 if (copyFlen==1) {
2684 FirstGuess ( copyF1,copyF2,copyFlen );
2685 ChooseFirstRotation ( copyF1[1],copyF2[1] );
2686 } else
2687 SuperposeSSGraphs ( G1,copyF1,G2,copyF2,copyFlen );
2688 for (i=0;i<4;i++)
2689 for (j=0;j<4;j++) {
2690 TMx0[i][j] = TMatrix[i][j];
2691 TMx [i][j] = TMatrix[i][j];
2692 }
2693
2694
2695 // 5. Make iterative atom-to-atom superposition
2696
2697 MMDB1->MakeBricks ( Calpha2,nres2,1.25*maxContact );
2698 rc = OptimizeNalign();
2699
2700
2701 // 6. Finalize the things
2702
2703 if (nalgn<=0) {
2704 for (i=0;i<4;i++)
2705 for (j=0;j<4;j++)
2706 TMatrix[i][j] = TMx0[i][j];
2707 for (i=0;i<nres1;i++) a1[i].c0 = -1;
2708 for (i=0;i<nres2;i++) a2[i].c0 = -1;
2709 rmsd_achieved = -1.0;
2710 nalgn = 0;
2711 ngaps = 0;
2712 seqIdent = 0.0;
2713 ncombs = 1.0;
2714 } else {
2715 CalcNGaps ( a1,nres1,i1,i );
2716 CalcNGaps ( a2,nres2,i2,j );
2717 ngaps = mmdb::IMax ( i1,i2 );
2718 nmd = mmdb::IMax ( i ,j );
2719 seqIdent = 0.0;
2720 for (i=0;i<nres1;i++)
2721 if (a1[i].c0>=0) {
2722 if (!strcasecmp(Calpha1[i]->GetResName(),
2723 Calpha2[a1[i].c0]->GetResName()))
2724 seqIdent += 1.0;
2725 }
2726 seqIdent = seqIdent/nalgn;
2727 ncombs = CalcNCombs ( G1,SSED1,nSSEs1,a1,nres1 ) *
2728 CalcNCombs ( G2,SSED2,nSSEs2,a2,nres2 );
2729 // get scores of matched SSEs:
2730 for (i=0;i<nSSEs1;i++)
2731 if (SSED1[i].match>0) {
2732 SSED1[i].Transform ( TMatrix );
2733 SSED1[i].CalcScore ( SSED2[SSED1[i].match-1] );
2734 CalcQScore ( SSED1[i] );
2735 }
2736 }
2737
2738
2739 // 6.1 Release memory
2740
2741 if (SDist) {
2742 for (i=0;i<SDistAlloc;i++)
2743 if (SDist[i]) delete[] SDist[i];
2744 delete[] SDist;
2745 SDist = NULL;
2746 }
2747 SDistAlloc = 0;
2748
2749 mmdb::FreeVectorMemory ( FH1,1 );
2750 mmdb::FreeVectorMemory ( FS1,1 );
2751 mmdb::FreeVectorMemory ( FH2,1 );
2752 mmdb::FreeVectorMemory ( FS2,1 );
2753 nFH1 = 0;
2754 nFS1 = 0;
2755 nFH2 = 0;
2756 nFS2 = 0;
2757
2758 mmdb::FreeMatrixMemory ( AD,AD_alloc,0,0 );
2759
2760 mmdb::FreeVectorMemory ( cax0,0 );
2761 mmdb::FreeVectorMemory ( cay0,0 );
2762 mmdb::FreeVectorMemory ( caz0,0 );
2763
2764 if (sdata) delete[] sdata;
2765
2766 }
2767
2768