1 // $Id: ssm_align.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_Align <implementation>
31 // ~~~~~~~~~
32 // **** Project : Structure alignment in 3D
33 // ~~~~~~~~~
34 // **** Classes : ssm::Align ( Secondary Structure Matching )
35 // ~~~~~~~~~ ssm::XAlign ( Output alignment )
36 // ssm::XTAlign ( Text output alignment )
37 //
38 // E. Krissinel, 2002-2014
39 //
40 // =================================================================
41 //
42
43 #include <string.h>
44 #include <ctype.h>
45
46 #include "mmdb2/mmdb_tables.h"
47 #include "ssm_align.h"
48
49 namespace ssm {
50 const char * SSM_DATE = "18-09-2013";
51 }
52
53 // --------------------------- ssm::Align ------------------------
54
Align()55 ssm::Align::Align() : mmdb::io::Stream() {
56 InitAlign();
57 }
58
Align(mmdb::io::RPStream Object)59 ssm::Align::Align ( mmdb::io::RPStream Object )
60 : mmdb::io::Stream ( Object ) {
61 InitAlign();
62 }
63
~Align()64 ssm::Align::~Align() {
65 FreeMemory();
66 }
67
FreeMemory()68 void ssm::Align::FreeMemory() {
69 mmdb::FreeVectorMemory ( Ca1 ,0 );
70 mmdb::FreeVectorMemory ( dist1 ,0 );
71 mmdb::FreeVectorMemory ( Ca2 ,0 );
72 mmdb::FreeVectorMemory ( pqvalues,0 );
73 if (G1) delete G1;
74 if (G2) delete G2;
75 G1 = NULL;
76 G2 = NULL;
77 nMatches = 0;
78 }
79
InitAlign()80 void ssm::Align::InitAlign() {
81
82 mmdb::Mat4Init ( TMatrix ); // transformation matrix
83
84 cnCheck = CONNECT_Flexible;
85 rmsd = 0.0; // core rmsd achieved
86 Qscore = 0.0; // core Q achieved
87 nres1 = 0; // number of residues in query structure
88 nres2 = 0; // number of residues in target structure
89 nsel1 = 0; // number of residues in query selection
90 nsel2 = 0; // number of residues in target selection
91 nalgn = 0; // number of aligned residues
92 ngaps = 0; // number of gaps
93 nmd = 0; // number of misdirections
94 seqIdentity = 0; // sequence identity
95 ncombs = 1.0; // number of SSE combinations
96
97 selHndCa1 = 0; // selection handle to used C-alphas in
98 // query structure
99 Ca1 = NULL; // C-alpha correspondence vector for
100 // query structure
101 dist1 = NULL; // optimized distances between the query
102 // and target C-alphas
103 selHndCa2 = 0; // selection handle to used C-alphas in
104 // target structure
105 Ca2 = NULL; // C-alpha correspondence vector for
106 // target structure
107
108 G1 = NULL;
109 G2 = NULL; // retained SSE graphs
110
111 pqvalues = NULL;
112 nMatches = 0;
113
114 }
115
116
MapSelections(int & selHndCa,mmdb::PManager M,PGraph G,int selHnd,mmdb::ivector & newID)117 void ssm::Align::MapSelections ( int & selHndCa, mmdb::PManager M,
118 PGraph G, int selHnd,
119 mmdb::ivector & newID ) {
120 mmdb::PPAtom a;
121 int nr,i,k;
122 G->SelectCalphas ( M,selHndCa,NULL );
123 if (selHnd>0) {
124 M->GetSelIndex ( selHndCa,a,nr );
125 mmdb::GetVectorMemory ( newID,nr,0 );
126 k = 0;
127 for (i=0;i<nr;i++)
128 if (a[i]->isInSelection(selHnd)) newID[i] = k++;
129 else newID[i] = -1;
130 M->Select ( selHndCa,mmdb::STYPE_ATOM,selHnd,mmdb::SKEY_AND );
131 } else
132 newID = NULL;
133 }
134
135
MakeSelections(mmdb::PManager M1,int selHnd1,mmdb::PManager M2,int selHnd2)136 void ssm::Align::MakeSelections ( mmdb::PManager M1, int selHnd1,
137 mmdb::PManager M2, int selHnd2 ) {
138 mmdb::ivector newID1,newID2;
139 int i,k;
140
141 MapSelections ( selHndCa1,M1,G1,selHnd1,newID1 );
142 MapSelections ( selHndCa2,M2,G2,selHnd2,newID2 );
143
144 if (newID2) {
145 k = 0;
146 if (newID1) {
147 for (i=0;i<nres2;i++)
148 if (newID2[i]>=0) {
149 if (Ca2[i]>=0) Ca2[k++] = newID1[Ca2[i]];
150 else Ca2[k++] = -1;
151 }
152 } else {
153 for (i=0;i<nres2;i++)
154 if (newID2[i]>=0) Ca2[k++] = Ca2[i];
155 }
156 }
157
158 if (newID1) {
159 k = 0;
160 if (newID2) {
161 for (i=0;i<nres1;i++)
162 if (newID1[i]>=0) {
163 if (Ca1[i]>=0) Ca1[k] = newID2[Ca1[i]];
164 else Ca1[k] = -1;
165 dist1[k] = dist1[i];
166 k++;
167 }
168 } else {
169 for (i=0;i<nres1;i++)
170 if (newID1[i]>=0) {
171 Ca1[k] = Ca1[i];
172 dist1[k] = dist1[i];
173 k++;
174 }
175 }
176 }
177
178 nsel1 = M1->GetSelLength ( selHndCa1 );
179 nsel2 = M2->GetSelLength ( selHndCa2 );
180
181 mmdb::FreeVectorMemory ( newID1,0 );
182 mmdb::FreeVectorMemory ( newID2,0 );
183
184 }
185
186
align(mmdb::PManager M1,mmdb::PManager M2,PRECISION precision,CONNECTIVITY connectivity,int selHnd1,int selHnd2)187 int ssm::Align::align ( mmdb::PManager M1, mmdb::PManager M2,
188 PRECISION precision,
189 CONNECTIVITY connectivity,
190 int selHnd1, int selHnd2 ) {
191 PPMatch Match;
192 mmdb::ivector F1,F2;
193 mmdb::realtype Q1;
194 int i,nm;
195
196 FreeMemory();
197
198 SetMatchPrecision ( precision );
199 SetConnectivityCheck ( connectivity );
200 cnCheck = connectivity;
201
202 U.SetUniqueMatch ( true );
203 U.SetBestMatch ( true );
204
205 G1 = GetSSGraph ( M1,selHnd1,i );
206 if (!G1) return i;
207
208 G2 = GetSSGraph ( M2,selHnd2,i );
209 if (!G2) return i+2;
210
211 U.MatchGraphs ( G1,G2,1 );
212
213 U.GetMatches ( Match,nMatches );
214 if (nMatches<=0) return RC_NoHits;
215
216 mmdb::GetVectorMemory ( pqvalues,nMatches,0 );
217 for (i=0;i<nMatches;i++)
218 pqvalues[i] = -1.0;
219
220 Qscore = -0.5;
221 for (i=0;i<nMatches;i++)
222 if (Match[i]) {
223 Match[i]->GetMatch ( F1,F2,nm );
224 superpose.SuperposeCalphas ( G1,G2,F1,F2,nm,M1,M2,
225 selHnd1,selHnd2 );
226 Q1 = superpose.GetCalphaQ();
227 if ((Q1>0.0) && (Q1>Qscore)) {
228 Qscore = Q1;
229 superpose.GetSuperposition ( Ca1,dist1,nres1,Ca2,nres2,
230 TMatrix,rmsd,nalgn,ngaps,
231 seqIdentity, nmd,ncombs );
232 }
233 pqvalues[i] = Q1;
234 }
235
236 if (Qscore>0.0) {
237 MakeSelections ( M1,selHnd1, M2,selHnd2 );
238 return RC_Ok;
239 }
240
241 return RC_NoSuperposition;
242
243 }
244
245
AlignSelectedMatch(mmdb::PManager M1,mmdb::PManager M2,PRECISION precision,CONNECTIVITY connectivity,int selHnd1,int selHnd2,int nselect)246 int ssm::Align::AlignSelectedMatch ( mmdb::PManager M1,
247 mmdb::PManager M2,
248 PRECISION precision,
249 CONNECTIVITY connectivity,
250 int selHnd1,
251 int selHnd2,
252 int nselect ) {
253 PPMatch Match;
254 mmdb::ivector F1,F2;
255 mmdb::realtype Q1;
256 int i,nGMatches,nm;
257
258 FreeMemory();
259
260 SetMatchPrecision ( precision );
261 SetConnectivityCheck ( connectivity );
262 cnCheck = connectivity;
263
264 U.SetUniqueMatch ( true );
265 U.SetBestMatch ( true );
266
267 G1 = GetSSGraph ( M1,selHnd1,i );
268 if (!G1) return i;
269
270 G2 = GetSSGraph ( M2,selHnd2,i );
271 if (!G2) return i+2;
272
273 U.MatchGraphs ( G1,G2,1 );
274
275 U.GetMatches ( Match,nGMatches );
276 if (nGMatches<=0) return RC_NoHits;
277
278 if (nselect>=nGMatches) {
279 // printf(" There are only %d matches for this alignment", nGMatches);
280 return RC_TooFewMatches;
281 }
282
283 Qscore = -0.5;
284
285 if (Match[nselect]) {
286
287 Match[nselect]->GetMatch ( F1,F2,nm );
288 superpose.SuperposeCalphas ( G1,G2,F1,F2,nm,M1,M2,
289 selHnd1,selHnd2 );
290
291 Q1 = superpose.GetCalphaQ();
292 if (Q1>0.0) {
293 superpose.GetSuperposition ( Ca1,dist1,nres1,Ca2,nres2,TMatrix,
294 rmsd,nalgn,ngaps,seqIdentity,
295 nmd,ncombs );
296 MakeSelections ( M1,selHnd1, M2,selHnd2 );
297 return RC_Ok;
298 }
299
300 }
301
302 return RC_NoSuperposition;
303
304 }
305
306
write(mmdb::io::RFile f)307 void ssm::Align::write ( mmdb::io::RFile f ) {
308 int i,j;
309
310 for (i=0;i<4;i++)
311 for (j=0;j<4;j++)
312 f.WriteReal ( &(TMatrix[i][j]) );
313
314 f.WriteInt ( &cnCheck );
315 f.WriteReal ( &rmsd );
316 f.WriteInt ( &nres1 );
317 f.WriteInt ( &nres2 );
318 f.WriteInt ( &nsel1 );
319 f.WriteInt ( &nsel2 );
320 f.WriteInt ( &nalgn );
321 f.WriteInt ( &ngaps );
322 f.WriteInt ( &nmd );
323 f.WriteReal ( &seqIdentity );
324 f.WriteReal ( &ncombs );
325
326 if (Ca1 && (nsel1>0))
327 for (i=0;i<nsel1;i++) {
328 f.WriteInt ( &(Ca1 [i]) );
329 f.WriteReal ( &(dist1[i]) );
330 }
331 if (Ca2 && (nsel2>0))
332 for (i=0;i<nsel2;i++)
333 f.WriteInt ( &(Ca2[i]) );
334
335 StreamWrite ( f,G1 );
336 StreamWrite ( f,G2 );
337
338 }
339
read(mmdb::io::RFile f)340 void ssm::Align::read ( mmdb::io::RFile f ) {
341 int i,j;
342
343 FreeMemory();
344
345 for (i=0;i<4;i++)
346 for (j=0;j<4;j++)
347 f.ReadReal ( &(TMatrix[i][j]) );
348
349 f.ReadInt ( &cnCheck );
350 f.ReadReal ( &rmsd );
351 f.ReadInt ( &nres1 );
352 f.ReadInt ( &nres2 );
353 f.ReadInt ( &nsel1 );
354 f.ReadInt ( &nsel2 );
355 f.ReadInt ( &nalgn );
356 f.ReadInt ( &ngaps );
357 f.ReadInt ( &nmd );
358 f.ReadReal ( &seqIdentity );
359 f.ReadReal ( &ncombs );
360
361 if (nsel1>0) {
362 mmdb::GetVectorMemory ( Ca1 ,nsel1,0 );
363 mmdb::GetVectorMemory ( dist1,nsel1,0 );
364 for (i=0;i<nsel1;i++) {
365 f.ReadInt ( &(Ca1 [i]) );
366 f.ReadReal ( &(dist1[i]) );
367 }
368 }
369 if (nsel2>0) {
370 mmdb::GetVectorMemory ( Ca2,nsel2,0 );
371 for (i=0;i<nsel2;i++)
372 f.ReadInt ( &(Ca2[i]) );
373 }
374
375 StreamRead ( f,G1 );
376 StreamRead ( f,G2 );
377
378 }
379
380 namespace ssm {
381 MakeStreamFunctions(Align)
382 }
383
384
385 // ----------------------------- ssm::XAlign --------------------------
386
XAlign()387 ssm::XAlign::XAlign() {
388 XBlock1 = NULL;
389 nBlock1 = 0;
390 XBlock2 = NULL;
391 nBlock2 = 0;
392 algnLen = 0;
393 }
394
~XAlign()395 ssm::XAlign::~XAlign() {
396 FreeMemory();
397 }
398
FreeMemory()399 void ssm::XAlign::FreeMemory() {
400 if (XBlock1) delete[] XBlock1;
401 if (XBlock2) delete[] XBlock2;
402 XBlock1 = NULL;
403 nBlock1 = 0;
404 XBlock2 = NULL;
405 nBlock2 = 0;
406 algnLen = 0;
407 }
408
409
customInit()410 void ssm::XAlign::customInit() {}
411
align(PGraph g1,mmdb::PPAtom Calpha1,mmdb::ivector Ca1,int nat1,PGraph g2,mmdb::PPAtom Calpha2,mmdb::ivector Ca2,int nat2,mmdb::rvector dist1,int & nr)412 void ssm::XAlign::align (
413 PGraph g1, mmdb::PPAtom Calpha1, mmdb::ivector Ca1, int nat1,
414 PGraph g2, mmdb::PPAtom Calpha2, mmdb::ivector Ca2, int nat2,
415 mmdb::rvector dist1, int & nr ) {
416 int i,j;
417
418 FreeMemory();
419
420 a1 = Ca1;
421 a2 = Ca2;
422 alpha1 = Calpha1;
423 alpha2 = Calpha2;
424 sg1 = g1;
425 sg2 = g2;
426 d1 = dist1;
427 na1 = nat1;
428 na2 = nat2;
429
430 nCols1 = makeXBlocks ( Ca1,nat1,XBlock1,nBlock1 );
431 nCols2 = makeXBlocks ( Ca2,nat2,XBlock2,nBlock2 );
432 nRows = nat1 + nat2 + 2;
433
434 maxdist = 0.0;
435 for (i=0;i<nat1;i++)
436 if (Ca1[i]>=0) {
437 if (dist1[i]>maxdist) maxdist = dist1[i];
438 }
439 if (maxdist<=1.0e-2) maxdist = 1.0;
440
441 customInit();
442 nr = 0;
443 for (i=0;i<nBlock1;i++)
444 for (j=0;j<nBlock2;j++)
445 alignXBlocks ( XBlock1[i],XBlock2[j],nr );
446
447 algnLen = nr;
448
449 }
450
451
makeXBlocks(mmdb::ivector Ca,int nat,RPXBlock xBlock,int & nBlocks)452 int ssm::XAlign::makeXBlocks ( mmdb::ivector Ca, int nat, RPXBlock xBlock,
453 int & nBlocks ) {
454 // Ca is considered as blocks of non-negative,
455 // increasing-by-one numbers Ca[i]>=0, and negative
456 // Ca[i]<0 surrounding them. Block boundaries are drawn
457 // at the middle of negative-Ca[i] gaps.
458 // nBlocks returns the number of such blocks, each block is
459 // identified by the initial and final indices i1 and i2, and
460 // by "index mass center" mc used for sorting.
461 // Returns the number of fold-columns.
462 PXBlock XB1;
463 int nAlloc,i,j,i1,i2,ip1,ip2,iv,k,icol;
464 mmdb::realtype mc;
465
466 if (xBlock) delete[] xBlock;
467 xBlock = NULL;
468 nBlocks = 0;
469 nAlloc = 0;
470
471 i = 0;
472 i1 = 0; // begining of a block
473 // begining of first block, check for leading negatives
474 while (i<nat)
475 if (Ca[i]<0) i++;
476 else break;
477 do {
478 if (i<nat) {
479 // check for increasing-by-one positives
480 ip1 = i;
481 iv = Ca[i++];
482 mc = iv;
483 k = 1;
484 while (i<nat)
485 if (Ca[i]==iv+1) {
486 iv = Ca[i++];
487 mc += iv;
488 k++;
489 } else
490 break;
491 mc /= k;
492 ip2 = i-1; // increasing-by-one has stopped, check for negatives
493 while (i<nat)
494 if (Ca[i]<0) i++;
495 else break;
496 // get i2 as end of a block
497 if (i>=nat) i2 = nat-1; // the last block
498 else if (i-ip2>1) i2 = (ip2+i)/2; // take the medium
499 else i2 = ip2;
500 } else {
501 i2 = nat-1; // the only empty (all-negatives) block
502 ip1 = -1;
503 ip2 = -1;
504 mc = 0.0;
505 }
506 // create new block
507 if (nBlocks>=nAlloc) {
508 nAlloc += 20;
509 XB1 = new XBlock[nAlloc];
510 for (j=0;j<nBlocks;j++) {
511 XB1[j].i1 = xBlock[j].i1;
512 XB1[j].i2 = xBlock[j].i2;
513 XB1[j].ip1 = xBlock[j].ip1;
514 XB1[j].ip2 = xBlock[j].ip2;
515 XB1[j].mc = xBlock[j].mc;
516 XB1[j].icol = xBlock[j].icol;
517 }
518 delete[] xBlock;
519 xBlock = XB1;
520 }
521 xBlock[nBlocks].i1 = i1;
522 xBlock[nBlocks].i2 = i2;
523 xBlock[nBlocks].ip1 = ip1;
524 xBlock[nBlocks].ip2 = ip2;
525 xBlock[nBlocks].mc = mc;
526 xBlock[nBlocks].icol = 0;
527 nBlocks++;
528 i1 = i2+1;
529 } while (i<nat);
530
531 // assign fold-columns to the blocks
532 icol = 0;
533 do {
534 icol++;
535 i = 0;
536 iv = 0;
537 while (i<nBlocks) {
538 mc = mmdb::MaxReal;
539 k = -1;
540 for (j=i;j<nBlocks;j++)
541 if ((xBlock[j].icol==0) && (xBlock[j].mc<mc)) {
542 mc = xBlock[j].mc;
543 k = j;
544 }
545 if (k>=0) {
546 xBlock[k].icol = icol;
547 i = k+1;
548 iv = 1;
549 } else
550 i = nBlocks;
551 }
552 } while (iv);
553
554 return icol-1;
555
556 }
557
558
alignXBlocks(RXBlock B1,RXBlock B2,int & nr)559 void ssm::XAlign::alignXBlocks ( RXBlock B1, RXBlock B2, int & nr ) {
560 int l1,l2, i1,i2, sseType1,sseType2, icol;
561
562 if (((a1[B1.ip1]>=B2.ip1) && (a1[B1.ip1]<=B2.ip2)) ||
563 ((a1[B1.ip2]>=B2.ip1) && (a1[B1.ip2]<=B2.ip2)) ||
564 ((a2[B2.ip1]>=B1.ip1) && (a2[B2.ip1]<=B1.ip2)) ||
565 ((a2[B2.ip2]>=B1.ip1) && (a2[B2.ip2]<=B1.ip2))) {
566
567 if (a1[B1.ip1]<B2.ip1) {
568 l1 = 0;
569 l2 = B2.ip1 - B2.i1;
570 i1 = a2[B2.ip1];
571 i2 = B2.i1;
572 } else if (a1[B1.ip1]==B2.ip1) {
573 l1 = B1.ip1 - B1.i1;
574 l2 = B2.ip1 - B2.i1;
575 i1 = B1.i1;
576 i2 = B2.i1;
577 } else {
578 l1 = B1.ip1 - B1.i1; // number of leading unmappings
579 l2 = 0;
580 i1 = B1.i1;
581 i2 = a1[B1.ip1];
582 }
583
584 icol = B1.icol;
585
586 while (l1>l2) {
587 if (alpha1[i1])
588 sseType1 = sg1->GetSSEType ( alpha1[i1]->GetChainID(),i1 );
589 else sseType1 = V_UNKNOWN;
590 makeRow ( alpha1[i1],sseType1,NULL,V_UNKNOWN,
591 d1[i1],nr++,icol,false );
592 i1++;
593 l1--;
594 }
595 while (l2>l1) {
596 if (alpha2[i2])
597 sseType2 = sg2->GetSSEType ( alpha2[i2]->GetChainID(),i2 );
598 else sseType2 = V_UNKNOWN;
599 makeRow ( NULL,V_UNKNOWN,alpha2[i2++],sseType2,
600 -1.0,nr++,icol,false );
601 l2--;
602 }
603 while (l2>0) {
604 if (alpha1[i1])
605 sseType1 = sg1->GetSSEType ( alpha1[i1]->GetChainID(),i1 );
606 else sseType1 = V_UNKNOWN;
607 if (alpha2[i2])
608 sseType2 = sg2->GetSSEType ( alpha2[i2]->GetChainID(),i2 );
609 else sseType2 = V_UNKNOWN;
610 makeRow ( alpha1[i1],sseType1,alpha2[i2++],sseType2,
611 d1[i1],nr++,icol,false );
612 i1++;
613 l2--;
614 }
615
616 l1 = mmdb::IMin ( B1.ip2-i1, B2.ip2-i2 ) + 1;
617 while (l1>0) {
618 if (alpha1[i1])
619 sseType1 = sg1->GetSSEType ( alpha1[i1]->GetChainID(),i1 );
620 else sseType1 = V_UNKNOWN;
621 if (alpha2[i2])
622 sseType2 = sg2->GetSSEType ( alpha2[i2]->GetChainID(),i2 );
623 else sseType2 = V_UNKNOWN;
624 makeRow ( alpha1[i1],sseType1,alpha2[i2++],sseType2,
625 d1[i1],nr++,icol,true );
626 i1++;
627 l1--;
628 }
629
630 if (i1<=B1.ip2) {
631 l1 = 0;
632 l2 = B2.i2 - i2 + 1;
633 } else if (i2<=B2.ip2) {
634 l1 = B1.i2 - i1 + 1;
635 l2 = 0;
636 } else {
637 l1 = B1.i2 - i1 + 1;
638 l2 = B2.i2 - i2 + 1;
639 }
640 while ((l1>0) && (l2>0)) {
641 if (alpha1[i1])
642 sseType1 = sg1->GetSSEType ( alpha1[i1]->GetChainID(),i1 );
643 else sseType1 = V_UNKNOWN;
644 if (alpha2[i2])
645 sseType2 = sg2->GetSSEType ( alpha2[i2]->GetChainID(),i2 );
646 else sseType2 = V_UNKNOWN;
647 makeRow ( alpha1[i1],sseType1,alpha2[i2++],sseType2,
648 d1[i1],nr++,icol,false );
649 i1++;
650 l1--;
651 l2--;
652 }
653 while (l1>0) {
654 if (alpha1[i1])
655 sseType1 = sg1->GetSSEType ( alpha1[i1]->GetChainID(),i1 );
656 else sseType1 = V_UNKNOWN;
657 makeRow ( alpha1[i1],sseType1,NULL,V_UNKNOWN,
658 d1[i1],nr++,icol,false );
659 i1++;
660 l1--;
661 }
662 while (l2>0) {
663 if (alpha2[i2])
664 sseType2 = sg2->GetSSEType ( alpha2[i2]->GetChainID(),i2 );
665 else sseType2 = V_UNKNOWN;
666 makeRow ( NULL,V_UNKNOWN,alpha2[i2++],sseType2,
667 -1.0,nr++,icol,false );
668 l2--;
669 }
670
671 }
672
673 }
674
675
makeRow(mmdb::PAtom A1,int sseType1,mmdb::PAtom A2,int sseType2,mmdb::realtype dist,int rowNo,int icol,bool aligned)676 void ssm::XAlign::makeRow ( mmdb::PAtom A1, int sseType1,
677 mmdb::PAtom A2, int sseType2,
678 mmdb::realtype dist, int rowNo,
679 int icol, bool aligned ) {
680 UNUSED_ARGUMENT(A1);
681 UNUSED_ARGUMENT(sseType1);
682 UNUSED_ARGUMENT(A2);
683 UNUSED_ARGUMENT(sseType2);
684 UNUSED_ARGUMENT(dist);
685 UNUSED_ARGUMENT(rowNo);
686 UNUSED_ARGUMENT(icol);
687 UNUSED_ARGUMENT(aligned);
688 }
689
690
691
692 // ---------------------------- XTAlign --------------------------
693
694 namespace ssm {
695
PrintAtom(mmdb::io::RFile f,int sseType,mmdb::realtype hydropathy,mmdb::ChainID chID,mmdb::ResName resName,int seqNum,mmdb::InsCode insCode)696 void PrintAtom ( mmdb::io::RFile f, int sseType,
697 mmdb::realtype hydropathy,
698 mmdb::ChainID chID, mmdb::ResName resName,
699 int seqNum, mmdb::InsCode insCode ) {
700 char sse[2],hp[2],ch[3],S[200];
701
702 if (sseType==V_HELIX) sse[0] = 'H';
703 else if (sseType==V_STRAND) sse[0] = 'S';
704 else sse[0] = ' ';
705 sse[1] = char(0);
706
707 if ((-5.0<hydropathy) && (hydropathy<5.0)) {
708 if (hydropathy>=-0.5) hp[0] = '-';
709 else if (hydropathy<=-1.5) hp[0] = '+';
710 else hp[0] = '.';
711 } else
712 hp[0] = ' ';
713 hp[1] = char(0);
714
715 if ((!chID[0]) || (chID[0]==' ')) {
716 ch[0] = ' '; ch[1] = ' ';
717 } else {
718 ch[0] = chID[0]; ch[1] = ':';
719 }
720 ch[2] = char(0);
721
722 sprintf ( S," |%1s%1s %2s%3s%4i%1s|",
723 sse,hp,ch,resName,seqNum,insCode );
724 f.Write ( S );
725
726 }
727
728 }
729
Print(mmdb::io::RFile f)730 void ssm::XTAlign::Print ( mmdb::io::RFile f ) {
731 char S[100],SI[10];
732 int i;
733
734 if (alignKey<4) {
735
736 if (alignKey!=2)
737 PrintAtom ( f,sseType1,hydropathy1,chID1,
738 resName1,seqNum1,insCode1 );
739 else
740 f.Write ( " | |" );
741
742 if (alignKey==0) {
743 switch (simindex) {
744 case 5 : strcpy ( SI,"*****" ); break;
745 case 4 : strcpy ( SI,"+++++" ); break;
746 case 3 : strcpy ( SI,"=====" ); break;
747 case 2 : strcpy ( SI,"-----" ); break;
748 case 1 : strcpy ( SI,":::::" ); break;
749 default:
750 case 0 : strcpy ( SI,"....." ); break;
751 }
752 SI[1] = char(0);
753 sprintf ( S," %1s%5.2f%1s%1s",SI,dist,SI,SI );
754 if (S[3]==' ') S[3] = SI[0];
755 f.Write ( S );
756 SI[1] = SI[0];
757 for (i=1;i<loopNo;i++) f.Write ( SI );
758 // f.Write ( " " );
759 } else {
760 f.Write ( " " );
761 for (i=1;i<loopNo;i++) f.Write ( " " );
762 f.Write ( " " );
763 }
764
765 if (alignKey!=3)
766 PrintAtom ( f,sseType2,hydropathy2,chID2,
767 resName2,seqNum2,insCode2 );
768 else
769 f.Write ( " | |" );
770
771 f.LF();
772
773 }
774
775 }
776
777
XAlignText()778 ssm::XAlignText::XAlignText() : ssm::XAlign() {
779 R = NULL;
780 }
781
~XAlignText()782 ssm::XAlignText::~XAlignText() {
783 customFree();
784 }
785
customFree()786 void ssm::XAlignText::customFree() {
787 if (R) delete[] R;
788 R = NULL;
789 }
790
customInit()791 void ssm::XAlignText::customInit() {
792 int i;
793 customFree();
794 R = new ssm::XTAlign[nRows];
795 for (i=0;i<nRows;i++)
796 R[i].alignKey = 5;
797 }
798
WipeTextRows()799 void ssm::XAlignText::WipeTextRows() {
800 R = NULL;
801 }
802
makeRow(mmdb::PAtom A1,int sseType1,mmdb::PAtom A2,int sseType2,mmdb::realtype dist,int rowNo,int icol,bool aligned)803 void ssm::XAlignText::makeRow ( mmdb::PAtom A1, int sseType1,
804 mmdb::PAtom A2, int sseType2,
805 mmdb::realtype dist, int rowNo, int icol,
806 bool aligned ) {
807
808 if (aligned) R[rowNo].alignKey = 0;
809 else R[rowNo].alignKey = 1;
810
811 if (A1) {
812 R[rowNo].sseType1 = sseType1;
813 R[rowNo].hydropathy1 = A1->GetAAHydropathy();
814 R[rowNo].seqNum1 = A1->GetSeqNum ();
815 strcpy ( R[rowNo].chID1 ,A1->GetChainID() );
816 strcpy ( R[rowNo].resName1,A1->GetResName() );
817 strcpy ( R[rowNo].insCode1,A1->GetInsCode() );
818 } else
819 R[rowNo].alignKey = 2;
820
821 if (A2) {
822 R[rowNo].sseType2 = sseType2;
823 R[rowNo].hydropathy2 = A2->GetAAHydropathy();
824 R[rowNo].seqNum2 = A2->GetSeqNum ();
825 strcpy ( R[rowNo].chID2 ,A2->GetChainID() );
826 strcpy ( R[rowNo].resName2,A2->GetResName() );
827 strcpy ( R[rowNo].insCode2,A2->GetInsCode() );
828 } else
829 R[rowNo].alignKey = 3;
830
831 if ((!A1) && (!A2)) R[rowNo].alignKey = 4;
832
833 R[rowNo].simindex = -5;
834 R[rowNo].dist = -1.0;
835 if (aligned) {
836 if (A1 && A2) R[rowNo].simindex = A1->GetAASimilarity ( A2 );
837 else R[rowNo].simindex = -5;
838 R[rowNo].dist = dist;
839 }
840
841 R[rowNo].loopNo = icol;
842
843 }
844
GetAlignments(mmdb::pstr & algn1,mmdb::pstr & algn2)845 void ssm::XAlignText::GetAlignments ( mmdb::pstr & algn1,
846 mmdb::pstr & algn2 ) {
847 char rn1[10];
848 char rn2[10];
849 int i;
850 if (algn1) delete[] algn1;
851 if (algn2) delete[] algn2;
852 if (algnLen>0) {
853 algn1 = new char[algnLen+1];
854 algn2 = new char[algnLen+1];
855 for (i=0;i<algnLen;i++) {
856 if (R[i].alignKey<=3) {
857 if (R[i].alignKey!=2)
858 mmdb::Get1LetterCode ( R[i].resName1,rn1 );
859 else strcpy ( rn1,"-" );
860 if (R[i].alignKey!=3)
861 mmdb::Get1LetterCode ( R[i].resName2,rn2 );
862 else strcpy ( rn2,"-" );
863 if (R[i].alignKey==0) {
864 rn1[0] = char(toupper(int(rn1[0])));
865 rn2[0] = char(toupper(int(rn2[0])));
866 } else {
867 rn1[0] = char(tolower(int(rn1[0])));
868 rn2[0] = char(tolower(int(rn2[0])));
869 }
870 } else {
871 strcpy ( rn1,"-" );
872 strcpy ( rn2,"-" );
873 }
874 algn1[i] = rn1[0];
875 algn2[i] = rn2[0];
876 }
877 algn1[algnLen] = char(0);
878 algn2[algnLen] = char(0);
879 } else {
880 algn1 = NULL;
881 algn2 = NULL;
882 }
883 }
884
885
886 namespace ssm {
887
PrintAlignTable(mmdb::io::RFile f,mmdb::PManager M1,mmdb::PManager M2,PAlign SSMAlign)888 void PrintAlignTable ( mmdb::io::RFile f,
889 mmdb::PManager M1, mmdb::PManager M2,
890 PAlign SSMAlign ) {
891 XAlignText CXA;
892 PXTAlign XTA;
893 mmdb::PPAtom Calpha1,Calpha2;
894 int nat1,nat2,nr,j;
895
896 M1->GetSelIndex ( SSMAlign->selHndCa1,Calpha1,nat1 );
897 M2->GetSelIndex ( SSMAlign->selHndCa2,Calpha2,nat2 );
898
899 CXA.align ( SSMAlign->G1,Calpha1,SSMAlign->Ca1,nat1,
900 SSMAlign->G2,Calpha2,SSMAlign->Ca2,nat2,
901 SSMAlign->dist1,nr );
902 f.LF();
903
904 if (SSMAlign->cnCheck!=CONNECT_None) {
905 f.WriteLine ( " .-------------.----------.-------------." );
906 f.WriteLine ( " | Query | Dist.(A) | Target |" );
907 f.WriteLine ( " |-------------+----------+-------------|" );
908 } else {
909 f.WriteLine (
910 " .-------------.----------.-----------------------------------");
911 f.WriteLine (
912 " | Query | Dist.(A) | Target" );
913 f.WriteLine (
914 " |-------------+----------+-----------------------------------");
915 }
916
917 XTA = CXA.GetTextRows();
918 for (j=0;j<nr;j++)
919 XTA[j].Print ( f );
920
921 if (SSMAlign->cnCheck!=CONNECT_None)
922 f.WriteLine ( " `-------------'----------'-------------'" );
923 else
924 f.WriteLine (
925 " `-------------'----------'-----------------------------------");
926 f.LF();
927
928 f.WriteLine ( " Notations:" );
929 f.WriteLine ( " S/H residue belongs to a strand/helix" );
930 f.WriteLine ( " +/-/. hydrophylic/hydrophobic/neutral residue" );
931 f.WriteLine ( " ** identical residues matched: similarity 5" );
932 f.WriteLine ( " ++ similarity 4" );
933 f.WriteLine ( " == similarity 3" );
934 f.WriteLine ( " -- similarity 2" );
935 f.WriteLine ( " :: similarity 1" );
936 f.WriteLine ( " .. dissimilar residues: similarity 0" );
937
938 }
939
940 }
941