2 Last changed Time-stamp: <2010-06-24 15:12:09 ivo>
3 c Christoph Flamm and Ivo L Hofacker
4 {xtof,ivo}@tbi.univie.ac.at
6 $Id: baum.c,v 1.9 2008/05/21 10:15:45 ivo Exp $
13 #include "fold_vars.h"
21 #define SAME_STRAND(I,J) (((I)>=cut_point)||((J)<cut_point))
22 #define ORDER(x,y) if ((x)->nummer>(y)->nummer) {tempb=x; x=y; y=tempb;}
24 /* item of structure ringlist */
25 typedef struct _baum {
26 int nummer; /* number of base in sequence */
27 char typ; /* 'r' virtualroot, 'p' or 'q' paired, 'u' unpaired */
28 unsigned short base; /* 0<->unknown, 1<->A, 2<->C, 3<->G, 4<->U */
36 static char UNUSED rcsid[]="$Id: baum.c,v 1.9 2008/05/21 10:15:45 ivo Exp $";
37 static short *pairList = NULL;
38 static short *typeList = NULL;
39 static short *aliasList = NULL;
40 static baum *rl = NULL; /* ringlist */
41 static baum *wurzl = NULL; /* virtualroot of ringlist-tree */
42 static char **ptype = NULL;
44 static int comp_struc(const void *A, const void *B);
45 /* PUBLIC FUNCTIONES */
46 void ini_or_reset_rl (void);
48 void update_tree (int i, int j);
49 void clean_up_rl (void);
51 /* PRIVATE FUNCTIONES */
52 static void ini_ringlist(void);
53 static void reset_ringlist(void);
54 static void struc2tree (char *struc);
55 static void close_bp_en (baum *i, baum *j);
56 static void close_bp (baum *i, baum *j);
57 static void open_bp (baum *i);
58 static void open_bp_en (baum *i);
59 static void inb (baum *root);
60 static void inb_nolp (baum *root);
61 static void dnb (baum *rli);
62 static void dnb_nolp (baum *rli);
63 static void fnb (baum *rli);
64 static void make_ptypes(const short *S);
65 /* debugging tool(s) */
67 static void rl_status(void);
70 /* convert structure in bracked-dot-notation to a ringlist-tree */
71 static void struc2tree(char *struc) {
73 int ipos, jpos, balance = 0;
76 struc_copy = (char *)calloc(GSV.len+1, sizeof(char));
78 strcpy(struc_copy,struc);
80 for (ipos = 0; ipos < GSV.len; ipos++) {
81 if (struc_copy[ipos] == ')') {
83 struc_copy[ipos] = '.';
85 while (struc_copy[--ipos] != '(');
86 struc_copy[ipos] = '.';
96 "struc2tree(): start structure is not balanced !\n%s\n%s\n",
101 GSV.currE = GSV.startE =
102 (float )energy_of_structure_pt(GAV.farbe,
103 pairList, typeList, aliasList, 0) / 100.0;
106 for(i = 0; i < GSV.len; i++) {
107 if (pairList[i+1]>i+1)
108 rl[i].loop_energy = loop_energy(pairList, typeList, aliasList,i+1);
110 wurzl->loop_energy = loop_energy(pairList, typeList, aliasList,0);
117 static void ini_ringlist(void) {
120 /* needed by function energy_of_struct_pt() from Vienna-RNA-1.4 */
121 pairList = (short *)calloc(GSV.len + 2, sizeof(short));
122 assert(pairList != NULL);
123 typeList = (short *)calloc(GSV.len + 2, sizeof(short));
124 assert(typeList != NULL);
125 aliasList = (short *)calloc(GSV.len + 2, sizeof(short));
126 assert(aliasList != NULL);
127 pairList[0] = typeList[0] = aliasList[0] = GSV.len;
128 ptype = (char **)calloc(GSV.len + 2, sizeof(char *));
129 assert(ptype != NULL);
130 for (i=0; i<=GSV.len; i++) {
131 ptype[i] = (char*)calloc(GSV.len + 2, sizeof(char));
132 assert(ptype[i] != NULL);
135 /* allocate virtual root */
136 wurzl = (baum *)calloc(1, sizeof(baum));
137 assert(wurzl != NULL);
138 /* allocate ringList */
139 rl = (baum *)calloc(GSV.len+1, sizeof(baum));
141 /* allocate PostOrderList */
143 /* initialize virtualroot */
146 /* connect virtualroot to ringlist-tree in down direction */
147 wurzl->down = &rl[GSV.len];
148 /* initialize post-order list */
152 /* initialize rest of ringlist-tree */
153 for(i = 0; i < GSV.len; i++) {
155 GAV.currform[i] = '.';
156 GAV.prevform[i] = 'x';
159 /* decode base to numeric value */
160 c = encode_char(GAV.farbe[i]);
161 rl[i].base = typeList[i+1] = c;
162 aliasList[i+1] = alias[typeList[i+1]];
163 /* astablish links for node of the ringlist-tree */
165 rl[i].next = &rl[i+1];
166 rl[i].prev = ((i == 0) ? &rl[GSV.len] : &rl[i-1]);
167 rl[i].up = rl[i].down = NULL;
169 GAV.currform[GSV.len] = GAV.prevform[GSV.len] = '\0';
170 make_ptypes(aliasList);
174 /* make ringlist circular in next, prev direction */
176 rl[i].prev = &rl[i-1];
177 /* make virtual basepair for virtualroot */
184 void ini_or_reset_rl(void) {
186 /* if there is no ringList-tree make a new one */
190 /* start structure */
191 struc2tree(GAV.startform);
192 GSV.currE = GSV.startE = energy_of_structure(GAV.farbe, GAV.startform, 0);
195 /* stop structure(s) */
199 qsort(GAV.stopform, GSV.maxS, sizeof(char *), comp_struc);
200 for (i = 0; i< GSV.maxS; i++)
201 GAV.sE[i] = energy_of_structure(GAV.farbe_full, GAV.stopform[i], 0);
206 initialize_cofold(GSV.len);
207 /* fold sequence to get Minimum free energy structure (Mfe) */
208 GAV.sE[0] = cofold(GAV.farbe_full, GAV.stopform[0]);
210 /* revaluate energy of Mfe (maye differ if --logML=logarthmic */
211 GAV.sE[0] = energy_of_structure(GAV.farbe_full, GAV.stopform[0], 0);
213 GSV.stopE = GAV.sE[0];
214 ini_nbList(strlen(GAV.farbe_full)*strlen(GAV.farbe_full));
217 /* reset ringlist-tree to start conditions */
219 if(GTV.start) struc2tree(GAV.startform);
221 GSV.currE = GSV.startE;
227 static void reset_ringlist(void) {
230 for(i = 0; i < GSV.len; i++) {
231 GAV.currform[i] = '.';
232 GAV.prevform[i] = 'x';
235 rl[i].next = &rl[i + 1];
236 rl[i].prev = ((i == 0) ? &rl[GSV.len] : &rl[i - 1]);
237 rl[i].up = rl[i].down = NULL;
240 rl[i].prev = &rl[i-1];
244 /* update ringlist-tree */
245 void update_tree(int i, int j) {
247 baum *rli, *rlj, *tempb;
249 if ( abs(i) < GSV.len) { /* >> single basepair move */
250 if ((i > 0) && (j > 0)) { /* insert */
253 close_bp_en(rli, rlj);
255 else if ((i < 0)&&(j < 0)) { /* delete */
261 if (i > 0) { /* i remains the same, j shifts */
267 close_bp_en(rli, rlj);
269 else { /* j remains the same, i shifts */
277 close_bp_en(rli, rlj);
280 } /* << single basepair move */
281 else { /* >> double basepair move */
282 if ((i > 0) && (j > 0)) { /* insert */
283 rli = &rl[i-GSV.len-2];
284 rlj = &rl[j-GSV.len-2];
285 close_bp_en(rli->next, rlj->prev);
286 close_bp_en(rli, rlj);
288 else if ((i < 0)&&(j < 0)) { /* delete */
290 rli = &rl[i-GSV.len-2];
292 open_bp_en(rli->next);
294 } /* << double basepair move */
298 /* open a particular base pair */
299 void open_bp(baum *i) {
301 baum *in; /* points to i->next */
303 /* change string representation */
304 GAV.currform[i->nummer] = '.';
305 GAV.currform[i->down->nummer] = '.';
307 /* change pairtable representation */
308 pairList[1 + i->nummer] = 0;
309 pairList[1 + i->down->nummer] = 0;
311 /* change tree representation */
315 i->next = i->down->next;
319 i->down = in->prev->up = NULL;
322 /* close a particular base pair */
323 void close_bp (baum *i, baum *j) {
325 baum *jn; /* points to j->next */
327 /* change string representation */
328 GAV.currform[i->nummer] = '(';
329 GAV.currform[j->nummer] = ')';
331 /* change pairtable representation */
332 pairList[1 + i->nummer] = 1+ j->nummer;
333 pairList[1 + j->nummer] = 1 + i->nummer;
335 /* change tree representation */
348 /* for a given tree, generate postorder-list */
349 static void make_poList (baum *root) {
353 if (!root) root = wurzl;
356 /* foreach base in ringlist ... */
357 for (rli = stop->next; rli != stop; rli = rli->next) {
358 /* ... explore subtee if bp found */
359 if (rli->typ == 'p') {
360 /* fprintf(stderr, "%d >%d<\n", poListop, rli->nummer); */
361 poList[poListop++] = rli;
362 if ( poListop > GSV.len+1 ) {
363 fprintf(stderr, "Something went wrong in make_poList()\n");
373 /* for a given ringlist, generate all structures
374 with one additional basepair */
375 static void inb(baum *root) {
378 int E_old, E_new_in, E_new_out;
379 baum *stop,*rli,*rlj;
381 E_old = root->loop_energy;
383 /* loop ringlist over all possible i positions */
384 for(rli=stop->next;rli!=stop;rli=rli->next){
385 /* potential i-position is already paired */
386 if(rli->typ=='p') continue;
387 /* loop ringlist over all possible j positions */
388 for(rlj=rli->next;rlj!=stop;rlj=rlj->next){
389 /* base pair must enclose at least 3 bases */
390 if(rlj->nummer - rli->nummer < MYTURN) continue;
391 /* potential j-position is already paired */
392 if(rlj->typ=='p') continue;
393 /* if i-j can form a base pair ... */
394 if(ptype[rli->nummer][rlj->nummer]){
395 /* close the base bair and ... */
397 E_new_in = loop_energy(pairList, typeList, aliasList,rli->nummer+1);
398 E_new_out = loop_energy(pairList, typeList, aliasList,root->nummer+1);
399 /* ... evaluate energy of the structure */
400 EoT = (int) (GSV.currE*100 + ((GSV.currE<0)?-0.4:0.4)) + E_new_in + E_new_out - E_old ;
401 /* assert(EoT == energy_of_struct_pt(GAV.farbe, pairList, typeList, aliasList)); */
402 /* open the base pair again... */
404 /* ... and put the move and the enegy
405 of the structure into the neighbour list */
406 update_nbList(1 + rli->nummer, 1 + rlj->nummer, EoT);
412 /* for a given ringlist, generate all structures (canonical)
413 with one additional base pair (BUT WITHOUT ISOLATED BASE PAIRS) */
414 static void inb_nolp(baum *root) {
417 baum *stop, *rli, *rlj;
420 /* loop ringlist over all possible i positions */
421 for (rli=stop->next;rli!=stop;rli=rli->next) {
422 /* potential i-position is already paired */
423 if (rli->typ=='p') continue;
424 /* loop ringlist over all possible j positions */
425 for (rlj=rli->next;rlj!=stop;rlj=rlj->next) {
426 /* base pair must enclose at least 3 bases */
427 if (rlj->nummer - rli->nummer < MYTURN) continue;
428 /* potential j-position is already paired */
429 if (rlj->typ=='p') continue;
430 /* if i-j can form a base pair ... */
431 if (ptype[rli->nummer][rlj->nummer]) {
432 /* ... and extends a helix ... */
433 if (((rli->prev==stop && rlj->next==stop) && stop->typ != 'x') ||
434 (rli->next == rlj->prev)) {
435 /* ... close the base bair and ... */
437 /* ... evaluate energy of the structure */
438 EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
439 /* open the base pair again... */
441 /* ... and put the move and the enegy
442 of the structure into the neighbour list */
443 update_nbList(1 + rli->nummer, 1 + rlj->nummer, EoT);
445 /* if double insertion is possible ... */
446 else if ((rlj->nummer - rli->nummer >= MYTURN+2)&&
447 (rli->next->typ != 'p' && rlj->prev->typ != 'p') &&
448 (rli->next->next != rlj->prev->prev) &&
449 (ptype[rli->next->nummer][rlj->prev->nummer])) {
450 /* close the two base bair and ... */
451 close_bp(rli->next, rlj->prev);
453 /* ... evaluate energy of the structure */
454 EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
455 /* open the two base pair again ... */
458 /* ... and put the move and the enegy
459 of the structure into the neighbour list */
460 update_nbList(1+rli->nummer+GSV.len+1, 1+rlj->nummer+GSV.len+1, EoT);
467 /* for a given ringlist, generate all structures
468 with one less base pair */
469 static void dnb(baum *rli){
471 int EoT, E_old_in, E_old_out, E_new;
477 /* ... evaluate energy of the structure */
479 for (r=rli->next; r->up==NULL; r=r->next);
480 E_old_in = rli->loop_energy;
481 E_old_out = r->up->loop_energy;
482 E_new = loop_energy(pairList,typeList,aliasList,r->up->nummer+1);
483 EoT = (int) (GSV.currE*100 + ((GSV.currE<0)?-0.4:0.4)) -
484 E_old_in - E_old_out + E_new;
486 /* assert(EoT== energy_of_struct_pt(GAV.farbe, pairList, typeList, aliasList));*/
488 update_nbList(-(1 + rli->nummer), -(1 + rlj->nummer), EoT);
491 /* for a given ringlist, generate all structures (canonical)
492 with one less base pair (BUT WITHOUT ISOLATED BASE PAIRS) */
493 static void dnb_nolp(baum *rli) {
497 baum *rlin = NULL; /* pointers to following pair in helix, if any */
499 baum *rlip = NULL; /* pointers to preceding pair in helix, if any */
504 /* immediate interior base pair ? */
505 if (rlj->next == rlj->prev) {
510 /* immediate exterior base pair and not virtualroot ? */
511 if (rli->prev == rli->next && rli->next->typ != 'x') {
512 rlip = rli->next->up;
516 /* double delete ? */
517 if (rlip==NULL && rlin && rljn->next != rljn->prev ) {
518 /* open the two base pairs ... */
521 /* ... evaluate energy of the structure ... */
522 EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
523 /* ... and put the move and the enegy
524 of the structure into the neighbour list ... */
525 update_nbList(-(1+rli->nummer+GSV.len+1),-(1+rlj->nummer+GSV.len+1), EoT);
526 /* ... and close the two base pairs again */
527 close_bp(rlin, rljn);
529 } else { /* single delete */
530 /* the following will work only if boolean expr are shortcicuited */
531 if (rlip==NULL || (rlip->prev == rlip->next && rlip->prev->typ != 'x'))
532 if (rlin ==NULL || (rljn->next == rljn->prev)) {
533 /* open the base pair ... */
535 /* ... evaluate energy of the structure ... */
536 EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
537 /* ... and put the move and the enegy
538 of the structure into the neighbour list ... */
539 update_nbList(-(1 + rli->nummer),-(1 + rlj->nummer), EoT);
540 /* and close the base pair again */
546 /* for a given ringlist, generate all structures
547 with one shifted base pair */
548 static void fnb(baum *rli) {
551 baum *rlj, *stop, *help_rli, *help_rlj;
555 /* examin interior loop of bp(ij); (.......)
557 for (rlj = stop->next; rlj != stop; rlj = rlj->next) {
558 /* prevent shifting to paired position */
559 if ((rlj->typ=='p')||(rlj->typ=='q')) continue;
560 /* j-position of base pair shifts to k position (ij)->(ik) i<k<j */
561 if ( (rlj->nummer-rli->nummer >= MYTURN)
562 && (ptype[rli->nummer][rlj->nummer]) ) {
563 /* open original basepair */
565 /* close shifted version of original basepair */
567 /* evaluate energy of the structure */
568 EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
569 /* put the move and the enegy of the structure into the neighbour list */
570 update_nbList(1+rli->nummer, -(1+rlj->nummer), EoT);
571 /* open shifted basepair */
573 /* restore original basepair */
576 /* i-position of base pair shifts to position k (ij)->(kj) i<k<j */
577 if ( (stop->nummer-rlj->nummer >= MYTURN)
578 && (ptype[stop->nummer][rlj->nummer]) ) {
579 /* open original basepair */
581 /* close shifted version of original basepair */
583 /* evaluate energy of the structure */
584 EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
585 /* put the move and the enegy of the structure into the neighbour list */
586 update_nbList(-(1 + rlj->nummer), 1 + stop->nummer, EoT);
587 /* open shifted basepair */
589 /* restore original basepair */
593 /* examin exterior loop of bp(ij); (.......)
594 i or j moves <- -> */
595 for (rlj=rli->next;rlj!=rli;rlj=rlj->next) {
596 if ((rlj->typ=='p') || (rlj->typ=='q' ) || (rlj->typ=='x')) continue;
597 x=rlj->nummer-rli->nummer;
599 /* j-position of base pair shifts to position k */
600 if ((x >= MYTURN) && (ptype[rli->nummer][rlj->nummer])) {
601 if (rli->nummer<rlj->nummer) {
609 /* open original basepair */
611 /* close shifted version of original basepair */
612 close_bp(help_rli,help_rlj);
613 /* evaluate energy of the structure */
614 EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
615 /* put the move and the enegy of the structure into the neighbour list */
616 update_nbList(1 + rli->nummer, -(1 + rlj->nummer), EoT);
617 /* open shifted base pair */
619 /* restore original basepair */
622 x = rlj->nummer-stop->nummer;
624 /* i-position of base pair shifts to position k */
625 if ((x >= MYTURN) && (ptype[stop->nummer][rlj->nummer])) {
626 if (stop->nummer < rlj->nummer) {
634 /* open original basepair */
636 /* close shifted version of original basepair */
637 close_bp(help_rli, help_rlj);
638 /* evaluate energy of the structure */
639 EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
640 /* put the move and the enegy of the structure into the neighbour list */
641 update_nbList(-(1 + rlj->nummer), 1 + stop->nummer, EoT);
642 /* open shifted basepair */
644 /* restore original basepair */
650 /* for a given tree (structure),
651 generate all neighbours according to moveset */
652 void move_it (void) {
656 energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0)/100.;
658 if ( GTV.noLP ) { /* canonical neighbours only */
660 for (i = 0; i < GSV.len; i++) {
662 if (pairList[i+1]>i+1) {
663 inb_nolp(rl+i); /* insert pair neighbours */
664 dnb_nolp(rl+i); /* delete pair neighbour */
668 else { /* all neighbours */
670 for (i = 0; i < GSV.len; i++) {
672 if (pairList[i+1]>i+1) {
673 inb(rl+i); /* insert pair neighbours */
674 dnb(rl+i); /* delete pair neighbour */
675 if ( GTV.noShift == 0 ) fnb(rl+i);
683 void clean_up_rl(void) {
685 free(pairList); pairList=NULL;
686 free(typeList); typeList = NULL;
687 free(aliasList); aliasList = NULL;
689 free(wurzl); wurzl=NULL;
690 for (i=0; i<=GSV.len; i++)
697 static int comp_struc(const void *A, const void *B) {
699 aE = (int)(100 * energy_of_structure(GAV.farbe_full, ((char **)A)[0], 0));
700 bE = (int)(100 * energy_of_structure(GAV.farbe_full, ((char **)B)[0], 0));
706 static void rl_status(void) {
710 printf("\n%s\n%s\n", GAV.farbe, GAV.currform);
711 for (i=0; i <= GSV.len; i++) {
712 printf("%2d %c %c %2d %2d %2d %2d\n",
714 i == GSV.len ? 'X': GAV.farbe[i],
716 rl[i].up==NULL?0:(rl[i].up)->nummer,
717 rl[i].down==NULL?0:(rl[i].down)->nummer,
718 (rl[i].prev)->nummer,
719 (rl[i].next)->nummer);
726 static void make_ptypes(const short *S) {
730 for (l=1; l<=2; l++) {
731 int type,ntype=0,otype=0;
733 if (!SAME_STRAND(i,j)) j=cut_point;
735 type = pair[S[i]][S[j]];
736 while ((i>=1)&&(j<=n)) {
737 if ((i>1)&&(j<n)) ntype = pair[S[i-1]][S[j+1]];
738 if (noLonelyPairs && (!otype) && (!ntype))
739 type = 0; /* i.j can only form isolated pairs */
740 ptype[i-1][j-1] = ptype[j-1][i-1] = (char) type;
748 static void close_bp_en (baum *i, baum *j) {
749 /* close bp and update energy */
752 i->loop_energy = loop_energy(pairList,typeList,aliasList,i->nummer+1);
753 for (r=i->next; r->up==NULL; r=r->next);
754 r->up->loop_energy = loop_energy(pairList,typeList,aliasList,r->up->nummer+1);
757 static void open_bp_en (baum *i) {
758 /* open bp and update energy */
762 for (r=i->next; r->up==NULL; r=r->next);
763 r->up->loop_energy = loop_energy(pairList,typeList,aliasList,r->up->nummer+1);