1 /* @source embpdb *************************************************************
2 **
3 ** Algorithms for handling protein structural data.
4 ** For use with the Atom, Chain and Pdb objects defined in ajpdb.h
5 ** Also for use with Hetent, Het, Vdwres, Vdwall, Cmap and Pdbtosp objects
6 ** (also in ajpdb.h).
7 **
8 ** @author CopyrightCopyright (c) 2004 Jon Ison
9 ** @version $Revision: 1.24 $
10 ** @modified $Date: 2012/07/14 14:52:40 $ by $Author: rice $
11 ** @@
12 **
13 ** This library is free software; you can redistribute it and/or
14 ** modify it under the terms of the GNU Lesser General Public
15 ** License as published by the Free Software Foundation; either
16 ** version 2.1 of the License, or (at your option) any later version.
17 **
18 ** This library is distributed in the hope that it will be useful,
19 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
20 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
21 ** Lesser General Public License for more details.
22 **
23 ** You should have received a copy of the GNU Lesser General Public
24 ** License along with this library; if not, write to the Free Software
25 ** Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
26 ** MA 02110-1301, USA.
27 **
28 ****************************************************************************/
29
30 #include "ajlib.h"
31
32 #include "embpdb.h"
33 #include "ajpdb.h"
34 #include "ajdomain.h"
35 #include "ajlist.h"
36
37 #include <math.h>
38
39
40
41
42 /* ======================================================================= */
43 /* ============================ private data ============================= */
44 /* ======================================================================= */
45
46
47
48
49 /* ======================================================================= */
50 /* ================= Prototypes for private functions ==================== */
51 /* ======================================================================= */
52
53
54
55
56 /* ======================================================================= */
57 /* ========================== private functions ========================== */
58 /* ======================================================================= */
59
60
61
62
63 /* ======================================================================= */
64 /* =========================== constructors ============================== */
65 /* ======================================================================= */
66
67
68
69
70 /* @section Constructors ****************************************************
71 **
72 ** All constructors return a pointer to a new instance. It is the
73 ** responsibility of the user to first destroy any previous instance. The
74 ** target pointer does not need to be initialised to NULL, but it is good
75 ** programming practice to do so anyway.
76 **
77 ****************************************************************************/
78
79
80
81
82 /* ======================================================================= */
83 /* =========================== destructors =============================== */
84 /* ======================================================================= */
85
86
87
88
89 /* @section Structure Destructors *******************************************
90 **
91 ** All destructor functions receive the address of the instance to be
92 ** deleted. The original pointer is set to NULL so is ready for re-use.
93 **
94 ****************************************************************************/
95
96
97
98
99 /* ======================================================================= */
100 /* ============================ Assignments ============================== */
101 /* ======================================================================= */
102
103
104
105
106 /* @section Assignments *****************************************************
107 **
108 ** These functions overwrite the instance provided as the first argument
109 ** A NULL value is always acceptable so these functions are often used to
110 ** create a new instance by assignment.
111 **
112 ****************************************************************************/
113
114
115
116
117 /* ======================================================================= */
118 /* ============================= Modifiers =============================== */
119 /* ======================================================================= */
120
121
122
123
124 /* @section Modifiers *******************************************************
125 **
126 ** These functions use the contents of an instance and update them.
127 **
128 ****************************************************************************/
129
130
131
132
133 /* ======================================================================= */
134 /* ========================== Operators ===================================*/
135 /* ======================================================================= */
136
137
138
139
140 /* @section Operators *******************************************************
141 **
142 ** These functions use the contents of an instance but do not make any
143 ** changes.
144 **
145 ****************************************************************************/
146
147
148
149
150 /* @func embPdbidToSp *********************************************************
151 **
152 ** Read a pdb identifier code and writes the equivalent swissprot identifier
153 ** code. Relies on list of Pdbtosp objects sorted by PDB code, which is
154 ** usually obtained by a call to ajPdbtospReadAllNew.
155 **
156 ** @param [r] pdb [const AjPStr] Pdb identifier code
157 ** @param [w] spr [AjPStr*] Swissprot identifier code
158 ** @param [r] list [const AjPList] Sorted list of Pdbtosp objects
159 **
160 ** @return [AjBool] True if a swissprot identifier code was found
161 ** for the Pdb code.
162 **
163 ** @release 2.9.0
164 ** @@
165 ****************************************************************************/
166
embPdbidToSp(const AjPStr pdb,AjPStr * spr,const AjPList list)167 AjBool embPdbidToSp(const AjPStr pdb, AjPStr *spr, const AjPList list)
168 {
169 AjPPdbtosp *arr = NULL; /* Array derived from list */
170 ajint dim = 0; /* Size of array */
171 ajint idx = 0; /* Index into array for the Pdb code */
172
173
174 if(!pdb || !list)
175 {
176 ajWarn("Bad args passed to embPdbidToSp");
177
178 return ajFalse;
179 }
180
181
182 dim = (ajuint) ajListToarray(list, (void ***) &(arr));
183
184 if(!dim)
185 {
186 ajWarn("Empty list passed to embPdbidToSp");
187
188 return ajFalse;
189 }
190
191
192 if( (idx = ajPdbtospArrFindPdbid(arr, dim, pdb))==-1)
193 return ajFalse;
194
195 ajStrAssignS(spr, arr[idx]->Spr[0]);
196
197 return ajTrue;
198 }
199
200
201
202
203 /* @func embPdbidToAcc ********************************************************
204 **
205 ** Read a pdb identifier code and writes the equivalent accession number.
206 ** Relies on list of Pdbtosp objects sorted by PDB code, which is usually
207 ** obtained by a call to ajPdbtospReadAllNew.
208 **
209 ** @param [r] pdb [const AjPStr] Pdb identifier code
210 ** @param [w] acc [AjPStr*] Accession number
211 ** @param [r] list [const AjPList] Sorted list of Pdbtosp objects
212 **
213 ** @return [AjBool] True if a swissprot identifier code was found for the
214 ** Pdb code.
215 **
216 ** @release 2.9.0
217 ** @@
218 ****************************************************************************/
219
embPdbidToAcc(const AjPStr pdb,AjPStr * acc,const AjPList list)220 AjBool embPdbidToAcc(const AjPStr pdb, AjPStr *acc, const AjPList list)
221 {
222 AjPPdbtosp *arr = NULL; /* Array derived from list */
223 AjPPdbtosp *arrfree = NULL;
224 ajint dim = 0; /* Size of array */
225 ajint idx = 0; /* Index into array for the Pdb code */
226
227
228 if(!pdb || !list)
229 {
230 ajWarn("Bad args passed to embPdbidToAcc");
231
232 return ajFalse;
233 }
234
235
236 dim = (ajuint) ajListToarray(list, (void ***) &(arr));
237
238 if(!dim)
239 {
240 ajWarn("Empty list passed to embPdbidToAcc");
241
242 return ajFalse;
243 }
244
245
246 if( (idx = ajPdbtospArrFindPdbid(arr, dim, pdb))==-1)
247 {
248 arrfree = (AjPPdbtosp*) arr;
249 AJFREE(arrfree);
250
251 return ajFalse;
252 }
253
254 ajStrAssignS(acc, arr[idx]->Acc[0]);
255 arrfree = (AjPPdbtosp*) arr;
256 AJFREE(arrfree);
257
258 return ajTrue;
259 }
260
261
262
263
264 /* @func embPdbidToScop ******************************************************
265 **
266 ** Writes a list of scop identifier codes for the domains that a Pdb object
267 ** contains. The domain data is taken from a list of scop objects.
268 **
269 ** @param [r] pdb [const AjPPdb] Pointer to pdb object
270 ** @param [r] list_allscop [const AjPList] Pointer to SCOP list of SCOP
271 ** classification objects
272 ** @param [w] list_pdbscopids [AjPList*] Pointer to list of scop domain ids
273 ** in the current pdb object
274 **
275 ** @return [AjBool] ajTrue on success
276 **
277 ** @release 2.9.0
278 ** @@
279 ****************************************************************************/
280
embPdbidToScop(const AjPPdb pdb,const AjPList list_allscop,AjPList * list_pdbscopids)281 AjBool embPdbidToScop(const AjPPdb pdb, const AjPList list_allscop,
282 AjPList *list_pdbscopids)
283 {
284
285 AjIList iter = NULL; /* List iterator for SCOP classification list */
286 AjPScop ptr = NULL;
287 AjPStr tmpPdbId = NULL;
288 AjPStr tmpDomId = NULL;
289 ajint found = 0;
290
291 iter=ajListIterNewread(list_allscop);
292
293
294 while((ptr=(AjPScop)ajListIterGet(iter)))
295 {
296 ajStrAssignS(&tmpPdbId, ptr->Pdb);
297 ajStrFmtLower(&tmpPdbId);
298
299 if(ajStrMatchS(pdb->Pdb, tmpPdbId))
300 {
301 ajStrAssignS(&tmpDomId, ptr->Entry);
302 ajStrFmtLower(&tmpDomId);
303 ajListPushAppend(*list_pdbscopids, tmpDomId);
304 tmpDomId = NULL;
305 found = 1;
306 }
307 }
308
309 ajListIterDel(&iter);
310 ajStrDel(&tmpPdbId);
311 ajStrDel(&tmpDomId);
312
313 if(found==1)
314 return ajTrue;
315
316 return ajFalse;
317 }
318
319
320
321
322 /* @func embAtomInContact *****************************************************
323 **
324 ** Determines whether two atoms are in physical contact
325 **
326 ** @param [r] atm1 [const AjPAtom] Atom 1 object
327 ** @param [r] atm2 [const AjPAtom] Atom 1 object
328 ** @param [r] thresh [float] Threshold contact distance
329 ** @param [r] vdw [const AjPVdwall] Vdwall object
330 **
331 ** Contact between two atoms is defined as when the van der Waals surface of
332 ** the first atom comes within the threshold contact distance (thresh) of
333 ** the van der Waals surface of the second atom.
334 **
335 ** @return [AjBool] True if the two atoms form contact
336 **
337 ** @release 2.9.0
338 ** @@
339 **
340 ****************************************************************************/
341
embAtomInContact(const AjPAtom atm1,const AjPAtom atm2,float thresh,const AjPVdwall vdw)342 AjBool embAtomInContact(const AjPAtom atm1, const AjPAtom atm2, float thresh,
343 const AjPVdwall vdw)
344 {
345 float val = 0.0;
346 float val1 = 0.0;
347
348
349
350 /* Check args */
351 if(!atm1 || !atm2 || !vdw)
352 {
353 ajWarn("Bad args passed to embAtomInContact");
354
355 return ajFalse;
356 }
357
358
359 val=((atm1->X - atm2->X) * (atm1->X - atm2->X)) +
360 ((atm1->Y - atm2->Y) * (atm1->Y - atm2->Y)) +
361 ((atm1->Z - atm2->Z) * (atm1->Z - atm2->Z));
362
363
364 /* This calculation uses square root
365 if((sqrt(val) - embVdwRad(atm1, vdw) -
366 embVdwRad(atm2, vdw)) <= thresh)
367 return ajTrue;
368 */
369
370
371 /* Same calculation avoiding square root */
372 val1 = embVdwRad(atm1, vdw) + embVdwRad(atm2, vdw) + thresh;
373
374 if(val <= (val1*val1))
375 return ajTrue;
376
377
378 return ajFalse;
379 }
380
381
382
383
384 /* @func embAtomDistance ******************************************************
385 **
386 ** Returns the distance (Angstroms) between two atoms.
387 **
388 ** @param [r] atm1 [const AjPAtom] Atom 1 object
389 ** @param [r] atm2 [const AjPAtom] Atom 1 object
390 ** @param [r] vdw [const AjPVdwall] Vdwall object
391 **
392 ** Returns the distance (Angstroms) between the van der Waals surface of two
393 ** atoms.
394 **
395 ** @return [float] Distance (Angstroms) between two atoms.
396 **
397 ** @release 2.9.0
398 ** @@
399 **
400 ****************************************************************************/
401
embAtomDistance(const AjPAtom atm1,const AjPAtom atm2,const AjPVdwall vdw)402 float embAtomDistance(const AjPAtom atm1, const AjPAtom atm2,
403 const AjPVdwall vdw)
404 {
405 float val = 0.0;
406 float val1 = 0.0;
407
408
409 val=((atm1->X - atm2->X) * (atm1->X - atm2->X)) +
410 ((atm1->Y - atm2->Y) * (atm1->Y - atm2->Y)) +
411 ((atm1->Z - atm2->Z) * (atm1->Z - atm2->Z));
412
413
414 /* This calculation uses square root */
415 val1= (float) (sqrt(val) - embVdwRad(atm1, vdw) - embVdwRad(atm2, vdw));
416
417 return val1;
418 }
419
420
421
422
423 /* ======================================================================= */
424 /* ============================== Casts ===================================*/
425 /* ======================================================================= */
426
427
428
429
430 /* @section Casts ***********************************************************
431 **
432 ** These functions examine the contents of an instance and return some
433 ** derived information. Some of them provide access to the internal
434 ** components of an instance. They are provided for programming convenience
435 ** but should be used with caution.
436 **
437 ****************************************************************************/
438
439
440
441
442 /* ======================================================================= */
443 /* =========================== Reporters ==================================*/
444 /* ======================================================================= */
445
446
447
448
449 /* @section Reporters *******************************************************
450 **
451 ** These functions return the contents of an instance but do not make any
452 ** changes.
453 **
454 ****************************************************************************/
455
456
457
458
459 /* @func embVdwRad ************************************************************
460 **
461 ** Returns the van der Waals radius of an atom. Returns 1.2 as default.
462 **
463 ** @param [r] atm [const AjPAtom] Atom object
464 ** @param [r] vdw [const AjPVdwall] Vdwall object
465 **
466 ** @return [float] van der Waals radius of the atom
467 **
468 ** @release 2.9.0
469 ** @@
470 **
471 ****************************************************************************/
472
embVdwRad(const AjPAtom atm,const AjPVdwall vdw)473 float embVdwRad(const AjPAtom atm, const AjPVdwall vdw)
474 {
475 ajuint x = 0U;
476 ajuint y = 0U;
477
478 for(x = 0U; x < vdw->N; x++)
479 for(y = 0; y < vdw->Res[x]->N; y++)
480 if(ajStrMatchS(atm->Atm, vdw->Res[x]->Atm[y]))
481 return vdw->Res[x]->Rad[y];
482
483 return (float) 1.2;
484 }
485
486
487
488
489 /* @func embPdbToIdx **********************************************************
490 **
491 ** Reads a Pdb object and writes an integer which gives the index into the
492 ** protein sequence for a residue with a specified pdb residue number and a
493 ** specified chain number.
494 **
495 ** @param [w] idx [ajint*] Residue number (index into sequence)
496 ** @param [r] pdb [const AjPPdb] Pdb object
497 ** @param [r] res [const AjPStr] Residue number (PDB numbering)
498 ** @param [r] chn [ajuint] Chain number
499 **
500 ** @return [AjBool] True on succcess (res was found in pdb object)
501 **
502 ** @release 2.9.0
503 ** @@
504 ****************************************************************************/
embPdbToIdx(ajint * idx,const AjPPdb pdb,const AjPStr res,ajuint chn)505 AjBool embPdbToIdx(ajint *idx, const AjPPdb pdb, const AjPStr res, ajuint chn)
506 {
507 AjIList iter = NULL;
508 AjPResidue residue = NULL;
509
510
511 if(!pdb || !(res) || !(idx))
512 {
513 ajWarn("Bad arg's passed to embPdbToIdx");
514
515 return ajFalse;
516 }
517
518 if((chn > pdb->Nchn) || (!pdb->Chains) || (chn<1))
519 {
520 ajWarn("Bad arg's passed to embPdbToIdx");
521
522 return ajFalse;
523 }
524
525
526 /* Initialise the iterator */
527 iter=ajListIterNewread(pdb->Chains[chn-1]->Residues);
528
529
530 /* Iterate through the list of residues */
531 while((residue = (AjPResidue)ajListIterGet(iter)))
532 {
533 if(residue->Chn!=chn)
534 continue;
535
536 /*
537 ** Hard-coded to work on model 1
538 ** Continue / break if a non-protein residue is found or model no. !=1
539 */
540 if(residue->Mod!=1)
541 break;
542
543 /* if(residue->Type!='P')
544 continue; */
545
546 /* If we have found the residue */
547 if(ajStrMatchS(res, residue->Pdb))
548 {
549 ajListIterDel(&iter);
550 *idx = residue->Idx;
551
552 return ajTrue;
553 }
554 }
555
556 ajWarn("Residue number not found in embPdbToIdx");
557 ajListIterDel(&iter);
558
559 return ajFalse;
560 }
561
562
563
564
565 /* ======================================================================= */
566 /* ========================== Input & Output ============================= */
567 /* ======================================================================= */
568
569
570
571
572 /* @section Input & output **************************************************
573 **
574 ** These functions are used for formatted input and output to file.
575 **
576 ****************************************************************************/
577
578
579
580
581 /* ======================================================================= */
582 /* ======================== Miscellaneous =================================*/
583 /* ======================================================================= */
584
585
586
587
588 /* @section Miscellaneous ***************************************************
589 **
590 ** These functions may have diverse functions that do not fit into the other
591 ** categories.
592 **
593 ****************************************************************************/
594
595
596
597
598 /* @func embPdbListHeterogens *************************************************
599 **
600 ** Function to create a list of arrays of Atom objects for ligands in the
601 ** current Pdb object (a single array for each ligand). An array of int's
602 ** giving the number of Atom objects in each array, is also written. The
603 ** number of ligands is also written.
604 **
605 ** @param [r] pdb [const AjPPdb] Pointer to pdb object
606 ** @param [w] list_heterogens [AjPList*] Pointer to list of heterogen Atom
607 ** arrays
608 ** @param [w] siz_heterogens [AjPInt*] Pointer to integer array of sizes
609 ** (number of Atom objects in each
610 ** array).
611 ** @param [w] nhet [ajint*] Number of arrays in the list that
612 ** was written.
613 ** @param [u] logfile [AjPFile] Log file for error messages
614 **
615 ** @return [AjBool] ajTrue on success
616 **
617 ** @release 2.9.0
618 ** @@
619 ****************************************************************************/
620
embPdbListHeterogens(const AjPPdb pdb,AjPList * list_heterogens,AjPInt * siz_heterogens,ajint * nhet,AjPFile logfile)621 AjBool embPdbListHeterogens(const AjPPdb pdb, AjPList *list_heterogens,
622 AjPInt *siz_heterogens, ajint *nhet,
623 AjPFile logfile)
624 {
625 /*
626 ** NOTE: EVERYTHING IN THE CLEAN PDB FILES IS CURRENTLY CHAIN
627 ** ASSOCIATED! THIS WILL BE CHANGED IN FUTURE
628 */
629 AjIList iter = NULL; /* Iterator for atoms in current pdb object */
630 AjPAtom hetat = NULL; /* Pointer to current Atom object */
631 ajuint i = 0U; /* Counter for chains */
632 ajint prev_gpn = -10000; /* Group number of atom object from
633 previous iteration */
634 AjPList GrpAtmList = NULL; /* List to hold atoms from the current
635 group */
636 AjPAtom *AtmArray = NULL; /* Array of atom objects */
637 ajint n=0; /* number of elements in AtmArray */
638 ajint grp_count = 0; /* No. of groups */
639 ajint arr_count = 0; /* Index for siz_heterogens */
640
641 /* Check args */
642 if((pdb==NULL)||(list_heterogens==NULL)||(siz_heterogens==NULL))
643 {
644 ajWarn("Bad args passed to embPdbListHeterogens\n");
645
646 return ajFalse;
647 }
648
649 if((!(*list_heterogens))||(!(*siz_heterogens)))
650 {
651 ajWarn("Bad args passed to embPdbListHeterogens\n");
652
653 return ajFalse;
654 }
655
656 if(pdb->Ngp>0)
657 ajFmtPrintF(logfile, "\tNGP:%d\n", pdb->Ngp);
658
659 if(pdb->Nchn>0)
660 {
661 for(i = 0U; i < pdb->Nchn; i++)
662 {
663 prev_gpn=-100000; /* Reset prev_gpn for each chain */
664 /* initialise iterator for pdb->Chains[i]->Atoms */
665 iter=ajListIterNewread(pdb->Chains[i]->Atoms);
666
667 /* Iterate through list of Atom objects */
668 while((hetat=(AjPAtom)ajListIterGet(iter)))
669 {
670 /* check for type */
671 if(hetat->Type != 'H')
672 continue;
673
674 /* TEST FOR A NEW GROUP */
675 if(prev_gpn != hetat->Gpn)
676 {
677 grp_count++;
678
679 if(GrpAtmList)
680 {
681 n=(ajuint) (ajListToarray(GrpAtmList, (void ***) &AtmArray));
682 ajListPushAppend(*list_heterogens, AtmArray);
683 /*
684 ** So that ajListToarray doesn't try and free the
685 ** non-NULL pointer
686 */
687 AtmArray=NULL;
688 ajIntPut(siz_heterogens, arr_count, n);
689 (*nhet)++;
690 ajListFree(&GrpAtmList);
691 GrpAtmList=NULL;
692 arr_count++;
693 }
694
695 GrpAtmList=ajListNew();
696 prev_gpn=hetat->Gpn;
697 } /* End of new group loop */
698
699 ajListPushAppend(GrpAtmList, (AjPAtom) hetat);
700 } /* End of list iteration loop */
701
702 /* Free list iterator */
703 ajListIterDel(&iter);
704
705 } /* End of chain for loop */
706
707 if(GrpAtmList)
708 {
709 n=(ajuint) (ajListToarray(GrpAtmList, (void ***) &AtmArray));
710 ajListPushAppend(*list_heterogens, AtmArray);
711 /*
712 ** So that ajListToarray doesn't try and free the non-NULL
713 ** pointer
714 */
715 AtmArray=NULL;
716 ajIntPut(siz_heterogens, arr_count, n);
717 (*nhet)++;
718 ajListFree(&GrpAtmList);
719 GrpAtmList=NULL;
720 }
721
722 GrpAtmList = NULL;
723 prev_gpn = -10000;
724
725 } /* End of chain loop */
726
727
728 return ajTrue;
729 }
730
731
732
733
734 /* @func embPdbResidueIndexI **************************************************
735 **
736 ** Reads a Pdb object and writes an integer array which gives the index into
737 ** the protein sequence for structured residues (residues for which electron
738 ** density was determined) for a given chain. The array length is of course
739 ** equal to the number of structured residues.
740 **
741 ** @param [r] pdb [const AjPPdb] Pdb object
742 ** @param [r] chn [ajuint] Chain number
743 ** @param [w] idx [AjPInt*] Index array
744 **
745 ** @return [AjBool] True on succcess
746 **
747 ** @release 3.0.0
748 ** @@
749 ****************************************************************************/
750
embPdbResidueIndexI(const AjPPdb pdb,ajuint chn,AjPInt * idx)751 AjBool embPdbResidueIndexI(const AjPPdb pdb, ajuint chn, AjPInt *idx)
752 {
753 AjIList iter = NULL;
754 AjPResidue res = NULL;
755 /* ajint this_rn = 0;
756 ajint last_rn = -1000; */
757 ajint resn = 0;
758
759
760 if(!pdb || !(*idx))
761 {
762 ajWarn("Bad arg's passed to embPdbResidueIndexI");
763
764 return ajFalse;
765 }
766
767 if((chn > pdb->Nchn) || (!pdb->Chains))
768 {
769 ajWarn("Bad arg's passed to embPdbResidueIndexI");
770
771 return ajFalse;
772 }
773
774
775 /* Initialise the iterator */
776 iter=ajListIterNewread(pdb->Chains[chn-1]->Residues);
777
778
779 /* Iterate through the list of residues */
780 while((res=(AjPResidue)ajListIterGet(iter)))
781 {
782 if(res->Chn!=chn)
783 continue;
784
785 /* Hard-coded to work on model 1 */
786 /* Continue / break if a non-protein residue is found or
787 model no. !=1 */
788 if(res->Mod!=1)
789 break;
790 /*
791 if(res->Type!='P')
792 continue; */
793
794
795 /* If we are onto a new residue */
796 /*
797 this_rn=res->Idx;
798 if(this_rn!=last_rn)
799 {
800 ajIntPut(&(*idx), resn++, res->Idx);
801 last_rn=this_rn;
802 }*/
803
804 ajIntPut(&(*idx), resn++, res->Idx);
805
806 }
807
808 if(resn==0)
809 {
810 ajWarn("Chain not found in embPdbResidueIndexI");
811 ajListIterDel(&iter);
812 return ajFalse;
813 }
814
815 ajListIterDel(&iter);
816
817 return ajTrue;
818 }
819
820
821
822
823 /* @func embPdbResidueIndexC **************************************************
824 **
825 ** Reads a Pdb object and writes an integer array which gives the index into
826 ** the protein sequence for structured residues (residues for which electron
827 ** density was determined) for a given chain. The array length is of course
828 ** equal to the number of structured residues.
829 **
830 ** @param [r] pdb [const AjPPdb] Pdb object
831 ** @param [r] chn [char] Chain identifier
832 ** @param [w] idx [AjPInt*] Index array
833 **
834 ** @return [AjBool] True on succcess
835 **
836 ** @release 3.0.0
837 ** @@
838 ****************************************************************************/
839
embPdbResidueIndexC(const AjPPdb pdb,char chn,AjPInt * idx)840 AjBool embPdbResidueIndexC(const AjPPdb pdb, char chn, AjPInt *idx)
841 {
842 ajuint chnn = 0U;
843
844 if(!ajPdbChnidToNum(chn, pdb, &chnn))
845 {
846 ajWarn("Chain not found in embPdbResidueIndexC");
847
848 return ajFalse;
849 }
850
851 if(!embPdbResidueIndexI(pdb, chnn, idx))
852 return ajFalse;
853
854 return ajTrue;
855 }
856
857
858
859
860 /* @func embPdbResidueIndexICA ************************************************
861 **
862 ** Reads a Pdb object and writes an integer array which gives the index into
863 ** the protein sequence for structured residues (residues for which electron
864 ** density was determined) for a given chain, EXCLUDING those residues for
865 ** which CA atoms are missing. The array length is of course equal to the
866 ** number of structured residues.
867 **
868 ** @param [r] pdb [const AjPPdb] Pdb object
869 ** @param [r] chn [ajuint] Chain number
870 ** @param [w] idx [AjPUint*] Index array
871 ** @param [w] nres [ajint*] Array length
872 **
873 ** @return [AjBool] True on succcess
874 **
875 ** @release 3.0.0
876 ** @@
877 ****************************************************************************/
878
embPdbResidueIndexICA(const AjPPdb pdb,ajuint chn,AjPUint * idx,ajint * nres)879 AjBool embPdbResidueIndexICA(const AjPPdb pdb,
880 ajuint chn, AjPUint *idx, ajint *nres)
881 {
882 AjIList iter = NULL;
883 AjPAtom atm = NULL;
884 ajint this_rn = 0;
885 ajint last_rn = -1000;
886 ajint resn = 0; /* Sequential count of residues */
887
888 if(!pdb || !(*idx))
889 {
890 ajWarn("Bad arg's passed to embPdbResidueIndexICA");
891
892 return ajFalse;
893 }
894
895 if((chn > pdb->Nchn) || (!pdb->Chains))
896 {
897 ajWarn("Bad arg's passed to embPdbResidueIndexICA");
898
899 return ajFalse;
900 }
901
902
903 /* Initialise the iterator */
904 iter=ajListIterNewread(pdb->Chains[chn-1]->Atoms);
905
906
907 /* Iterate through the list of atoms */
908 while((atm=(AjPAtom)ajListIterGet(iter)))
909 {
910 if(atm->Chn!=chn)
911 continue;
912
913 /*
914 ** Hard-coded to work on model 1
915 ** Continue / break if a non-protein atom is found or model no. !=1
916 */
917 if(atm->Mod!=1)
918 break;
919
920 if(atm->Type!='P')
921 continue;
922
923 /* If we are onto a new residue */
924 this_rn=atm->Idx;
925
926 if(this_rn!=last_rn && ajStrMatchC(atm->Atm, "CA"))
927 {
928 ajUintPut(&(*idx), resn++, atm->Idx);
929 last_rn=this_rn;
930 }
931 }
932
933
934 if(resn==0)
935 {
936 ajWarn("Chain not found in embPdbResidueIndexICA");
937 ajListIterDel(&iter);
938
939 return ajFalse;
940 }
941
942 *nres=resn;
943
944 ajListIterDel(&iter);
945
946 return ajTrue;
947 }
948
949
950
951
952 /* @func embPdbResidueIndexCCA ************************************************
953 **
954 ** Reads a Pdb object and writes an integer array which gives the index into
955 ** the protein sequence for structured residues (residues for which electron
956 ** density was determined) for a given chain, EXCLUDING those residues for
957 ** which CA atoms are missing. The array length is of course equal to the
958 ** number of structured residues.
959 **
960 ** @param [r] pdb [const AjPPdb] Pdb object
961 ** @param [r] chn [char] Chain identifier
962 ** @param [w] idx [AjPUint*] Index array
963 ** @param [w] nres [ajint*] Array length
964 **
965 ** @return [AjBool] True on succcess
966 **
967 ** @release 3.0.0
968 ** @@
969 ****************************************************************************/
970
embPdbResidueIndexCCA(const AjPPdb pdb,char chn,AjPUint * idx,ajint * nres)971 AjBool embPdbResidueIndexCCA(const AjPPdb pdb, char chn,
972 AjPUint *idx, ajint *nres)
973 {
974 ajuint chnn = 0U;
975
976 if(!ajPdbChnidToNum(chn, pdb, &chnn))
977 {
978 ajWarn("Chain not found in embPdbResidueIndexCCA");
979
980 return ajFalse;
981 }
982
983 if(!embPdbResidueIndexICA(pdb, chnn, idx, nres))
984 return ajFalse;
985
986
987 return ajTrue;
988 }
989
990
991
992
993 /* @func embStrideToThree *****************************************************
994 **
995 ** Reads a string that contains an 8-state STRIDE secondary structure
996 ** assignment and writes a string with the corresponding 3-state assignment.
997 ** The 8 states used in STRIDE are 'H' (alpha helix), 'G' (3-10 helix),
998 ** 'I' (Pi-helix), 'E' (extended conformation), 'B' or 'b' (isolated bridge),
999 ** 'T' (turn) or 'C' (coil, i.e. none of the above). The 3 states used
1000 ** are 'H' (STRIDE 'H', 'G' or 'I'), 'E' (STRIDE 'E', 'B' or 'b') and 'C'
1001 ** (STRIDE 'T' or 'C'). The string is allocated if necessary.
1002 **
1003 ** @param [w] to [AjPStr*] String to write
1004 ** @param [r] from [const AjPStr] String to read
1005 **
1006 ** @return [AjBool] True on succcess
1007 **
1008 ** @release 2.9.0
1009 ** @@
1010 ****************************************************************************/
1011
embStrideToThree(AjPStr * to,const AjPStr from)1012 AjBool embStrideToThree(AjPStr *to, const AjPStr from)
1013 {
1014 if(!from)
1015 {
1016 ajWarn("Bad args passed to embStrideToThree");
1017
1018 return ajFalse;
1019 }
1020 else
1021 ajStrAssignS(to, from);
1022
1023 ajStrExchangeKK(to, 'G', 'H');
1024 ajStrExchangeKK(to, 'I', 'H');
1025 ajStrExchangeKK(to, 'B', 'E');
1026 ajStrExchangeKK(to, 'b', 'E');
1027 ajStrExchangeKK(to, 'T', 'C');
1028
1029 return ajTrue;
1030 }
1031