1#ifndef RNAPUZZLER_CONFIG_H 2#define RNAPUZZLER_CONFIG_H 3 4/* 5 * RNApuzzler template config 6 * 7 * c Daniel Wiegreffe, Daniel Alexander, Dirk Zeckzer 8 * ViennaRNA package 9 */ 10#include <stdlib.h> 11#include <math.h> 12 13#include "ViennaRNA/utils/basic.h" 14 15#include "definitions.inc" 16#include "vector_math.inc" 17#include "../headers/configtree_struct.h" 18#include "../headers/tBaseInformation_struct.h" 19#include "ViennaRNA/plotting/RNApuzzler/RNApuzzler.h" 20 21PRIVATE double 22getArcAngle(const config *cfg, 23 const int currentArc); 24 25 26PRIVATE double 27getArcAngleDegree(const config *cfg, 28 const int currentArc); 29 30 31/** 32 * @brief cfgPrintConfig 33 * - prints the given config to the terminal 34 * @param config 35 * - config to print 36 */ 37PRIVATE void 38cfgPrintConfig(config *config); 39 40 41/** 42 * @brief cfgGenerateConfig 43 * - generates configurations for each loop 44 * @param pair_table 45 * - the RNA's pairing information 46 * @param baseInformation 47 * - array of tBaseInformation annotations (for saving configs) 48 * @param unpaired 49 * - default length of backbone lines 50 * @param paired 51 * - default distance between paired bases 52 */ 53PRIVATE void 54cfgGenerateConfig(const short *const pair_table, 55 tBaseInformation *baseInformation, 56 const double unpaired, 57 const double paired); 58 59 60/** 61 * @brief cfgCloneConfig 62 * - prints the given config to the terminal 63 * @param cfg 64 * - config to clone 65 * @return clone of cfg 66 */ 67PRIVATE config * 68cfgCloneConfig(const config *cfg); 69 70 71/** 72 * @brief cfgFreeConfig 73 * - prints the given config to the terminal 74 * @param cfg 75 * - config to free 76 */ 77PRIVATE void 78cfgFreeConfig(config *cfg); 79 80 81/** 82 * Function to apply a set of config changes. 83 * 84 * By passing -1.0 to radiusNew you will apply the minimum possible radius for this loop while not allowing to shrink the loop. 85 * In case it would actually shrink a default relative increase is applied to the radius. 86 * 87 * By passing 0.0 to radiusNew you set the minimum possible radius no matter if there would be some shrinking or not. 88 * @param loopName the node to change 89 * @param deltaCfg array of angles to change (in degree) 90 * @param radiusNew desired radius to set while applying deltas 91 * @param puzzler 92 */ 93PRIVATE double 94cfgApplyChanges(config *cfg, 95 const char loopName, 96 const double *deltaCfg, 97 const double radiusNew, 98 const vrna_plot_options_puzzler_t *puzzler); 99 100 101/** 102 * @brief cfgIsValid 103 * - check if config is valid 104 * @param config 105 * - the config to check 106 * @param deltaCfg 107 * - array of config changes. 108 * contains diff-values for each config angle. 109 * in degree format 110 * @return true iff config is valid 111 */ 112PRIVATE short 113cfgIsValid(config *config, 114 const double *deltaCfg); 115 116 117/** 118 * @brief intToMotiv 119 * - converts a given number into a readable format 120 * @param _int 121 * - motiv name as number (int) 122 * @return 123 * - motiv name as char 124 */ 125PRIVATE char 126intToMotiv(const int _int); 127 128 129/* 130 *-------------------------------------------------------------------------- 131 *--- create, copy, and free config ------------------------------------ 132 *-------------------------------------------------------------------------- 133 */ 134 135/** 136 * @brief cfgCreateConfig 137 * - constructor-like method for creating a config 138 * @param radius 139 * - radius used for drawing that loop 140 * @return 141 * - an initialized config struct 142 */ 143PRIVATE config * 144cfgCreateConfig(const double radius) 145{ 146 config *cfg = (config *)vrna_alloc(1 * sizeof(config)); 147 148 cfg->radius = radius; 149 cfg->minRadius = radius; 150 cfg->defaultRadius = radius; 151 152 cfg->cfgArcs = NULL; 153 cfg->numberOfArcs = 0; 154 155 return cfg; 156} 157 158 159/** 160 * @brief cfgCreateConfigArc 161 * - constructor-like method for adding a new config entry to a given config 162 * @param angle 163 * - angle (radiant) that can be found between stems 'from' and 'to' later on 164 * @param numberOfArcSegments 165 * - number of arc segments between stems 'from' and 'to' 166 * @return 167 * - an initialized configArc struct 168 */ 169PRIVATE configArc 170cfgCreateConfigArc(const double angle, 171 const int numberOfArcSegments) 172{ 173 configArc newConfigArc; 174 175 newConfigArc.numberOfArcSegments = numberOfArcSegments; 176 177 newConfigArc.arcAngle = angle; 178 179 return newConfigArc; 180} 181 182 183PRIVATE config * 184cfgCloneConfig(const config *cfg) 185{ 186 config *clonedCfg = (config *)vrna_alloc(sizeof(config)); 187 188 clonedCfg->radius = cfg->radius; 189 clonedCfg->minRadius = cfg->minRadius; 190 clonedCfg->defaultRadius = cfg->defaultRadius; 191 clonedCfg->numberOfArcs = cfg->numberOfArcs; 192 193 const int numberOfArcs = cfg->numberOfArcs; 194 clonedCfg->cfgArcs = (configArc *)vrna_alloc(numberOfArcs * sizeof(configArc)); 195 196 for (int currentArc = 0; currentArc < numberOfArcs; ++currentArc) { 197 clonedCfg->cfgArcs[currentArc].numberOfArcSegments = 198 cfg->cfgArcs[currentArc].numberOfArcSegments; 199 clonedCfg->cfgArcs[currentArc].arcAngle = cfg->cfgArcs[currentArc].arcAngle; 200 } 201 202 return clonedCfg; 203} 204 205 206PRIVATE void 207cfgFreeConfig(config *cfg) 208{ 209 free(cfg->cfgArcs); 210 free(cfg); 211} 212 213 214/* documentation at header file */ 215PRIVATE void 216cfgPrintConfig(config *cfg) 217{ 218 short useDegree = 1; 219 220 for (int currentArc = 0; currentArc < cfg->numberOfArcs; ++currentArc) { 221 double angle = getArcAngle(cfg, currentArc); 222 223 if (useDegree) 224 angle = toDegree(angle); 225 } 226} 227 228 229/* 230 *-------------------------------------------------------------------------- 231 *--- access to config elements ---------------------------------------- 232 *-------------------------------------------------------------------------- 233 */ 234PRIVATE double 235getArcAngle(const config *cfg, 236 const int currentArc) 237{ 238 return cfg->cfgArcs[currentArc].arcAngle; 239} 240 241 242PRIVATE double 243getArcAngleDegree(const config *cfg, 244 const int currentArc) 245{ 246 return toDegree(getArcAngle(cfg, currentArc)); 247} 248 249 250/* 251 *-------------------------------------------------------------------------- 252 *--- radius computation ----------------------------------------------- 253 *-------------------------------------------------------------------------- 254 */ 255 256/** 257 * Approximate the radius of a circle required to draw m base pairs 258 * with a distance of 'a' to each other, and n (unpaired) consecutive 259 * nucleotides with a distance of b over a specified angle. 260 * 261 * Uses a Newton iteration to approximate solution 262 * 263 * @param a paired 264 * @param b unpaired 265 * @param m #stems 266 * @param n #backbones 267 * @param angle angle 268 */ 269PRIVATE double 270approximateConfigArcRadius(const double a, 271 const double b, 272 const short m, 273 const short n, 274 double angle) 275{ 276 const short MAX_ITERATIONS = 1000; 277 278 /* 279 * calculation: 280 * 281 * be s the length of a line at the circle (paired or unpaired / a or b) 282 * the angle over such a single line be alpha 283 * alpha = angle / ( m + n ) 284 * 285 * for such a single line the following equation holds (where r is the radius of the circle) 286 * sin( alpha / 2 ) = ( s / 2 ) / r 287 * r = ( s / 2 ) / sin( alpha / 2 ) 288 * r = ( s / 2 ) / sin( ( angle / ( m + n ) ) / 2 ) 289 * 290 * now we replace s with a or b to get the upper or lower bound for the radius interval 291 */ 292 double lowerBound = (b / 2) / sin((angle / (m + n)) / 2); 293 double upperBound = (a / 2) / sin((angle / (m + n)) / 2); 294 /* printf("new: lower %f upper %f\n", lowerBound, upperBound); */ 295 double rtn = 0.5 * (lowerBound + upperBound); 296 297 /* 298 * there is a minimum valid radius! 299 * if rtn is smaller than 0.5*a or 0.5*b than the result will become nan 300 */ 301 rtn = fmax(rtn, 0.5 * a); 302 rtn = fmax(rtn, 0.5 * b); 303 304 /* 305 * printf("lower %f upper %f %f\n", lowerBound, upperBound, rtn); 306 * printf("Angle %f\n", angle); 307 */ 308 309 /* short isERROR = 0; */ 310 311 int j = 0; 312 313 for (j = 0; j < MAX_ITERATIONS; j++) { 314 double dx = 2 * (m * asin(a / (2 * rtn)) + n * asin(b / (2 * rtn)) - (angle / 2)) 315 / (-(a * m / (rtn * sqrt(rtn * rtn - a * a / 4)) + b * n / 316 (rtn * sqrt(rtn * rtn - b * b / 4)))); 317 rtn -= dx; 318 319 if (fabs(dx) < EPSILON_3) 320 321 break; 322 } 323 324 if (rtn < lowerBound) 325 rtn = lowerBound; 326 else if (rtn > upperBound) 327 rtn = upperBound; 328 329 return rtn; 330} 331 332 333PRIVATE double 334approximateConfigRadius(const config *cfg, 335 const double unpaired, 336 const double paired) 337{ 338 /* 339 * calculate a fitting radius for each arc without compressing or stretching arc segments 340 * return the maximum of those values as the radius fitting for the loop 341 */ 342 double r = 0; 343 344 for (int currentArc = 0; currentArc < cfg->numberOfArcs; ++currentArc) { 345 int stems = 1; 346 int numberOfArcSegments = (cfg->cfgArcs[currentArc]).numberOfArcSegments; 347 double angle = getArcAngle(cfg, currentArc); 348 349 double tempR = approximateConfigArcRadius(paired, unpaired, stems, numberOfArcSegments, angle); 350 351 if (tempR > r) 352 353 r = tempR; 354 } 355 356 return r; 357} 358 359 360/* -------------------------------------------------------------------------------------------------------------------------- */ 361 362/** 363 * @brief cfgGenerateDefaultConfig 364 * - generates a config that resembles a drawing without any 365 * constraints given by config input for the given loop 366 * @param pair_table 367 * - the RNA's pairing information 368 * @param start 369 * - index of the loop's first base 370 * @param unpaired 371 * - default length of backbone lines 372 * @param paired 373 * - default distance between paired bases 374 * @param radius 375 * - radius for the given loop 376 * @return 377 * - complete config for that loop 378 */ 379PRIVATE config * 380cfgGenerateDefaultConfig(const short *const pair_table, 381 const int start, 382 const int unpaired, 383 const int paired, 384 const double radius) 385{ 386 /* create loop configuration */ 387 config *cfg = cfgCreateConfig(radius); 388 389 /* compute angles for paired and unpaired bases */ 390 double anglePaired = 2 * asin(paired / (2 * radius)); /* angle over paired */ 391 double angleUnpaired = 2 * asin(unpaired / (2 * radius)); /* angle over unpaired */ 392 393 /* initialize values for first arc */ 394 int arcUnpaired = 0; 395 double angleArc; /* alpha + numBackbones * beta */ 396 397 /* start with first base after parent stem */ 398 int i = start + 1; 399 400 while (i <= pair_table[start]) { 401 /* until last base at parent stem */ 402 if (pair_table[i] == 0) { 403 /* arc */ 404 i++; 405 } else { 406 /* increment number of arcs */ 407 ++(cfg->numberOfArcs); 408 409 if (i != pair_table[start]) 410 /* skip subtree at stem */ 411 i = pair_table[i] + 1; 412 else 413 /* parent stem -> finish */ 414 break; 415 } 416 } 417 418 cfg->cfgArcs = (configArc *)vrna_alloc(cfg->numberOfArcs * sizeof(configArc)); 419 420 /* start with first base after parent stem */ 421 i = start + 1; 422 int currentArc = 0; 423 int numberOfArcSegments = 0; 424 425 while (i <= pair_table[start]) { 426 /* until last base at parent stem */ 427 if (pair_table[i] == 0) { 428 /* arc */ 429 arcUnpaired++; 430 i++; 431 } else { 432 /* stem: create arc */ 433 angleArc = anglePaired + (arcUnpaired + 1) * angleUnpaired; 434 numberOfArcSegments = arcUnpaired + 1; 435 cfg->cfgArcs[currentArc] = cfgCreateConfigArc(angleArc, numberOfArcSegments); 436 ++currentArc; 437 438 if (i != pair_table[start]) { 439 /* initialize values for next arc */ 440 arcUnpaired = 0; 441 442 /* skip subtree at stem */ 443 i = pair_table[i] + 1; 444 } else { 445 /* parent stem -> finish */ 446 break; 447 } 448 } 449 } 450 451 return cfg; 452} 453 454 455PRIVATE void 456cfgGenHandleStem(int baseNr, 457 const short *const pair_table, 458 tBaseInformation *baseInformation, 459 const double unpaired, 460 const double paired); 461 462 463/** 464 * @brief cfgGenHandleLoop 465 * - recursively iterates through the RNA and generates default configurations. 466 * Alternates with corresponding handleStem method. 467 * @param baseNr 468 * - index of the loop's first base 469 * @param pair_table 470 * - the RNA's pairing information 471 * @param baseInformation 472 * - array of tBaseInformation annotations (to save config) 473 */ 474PRIVATE void 475cfgGenHandleLoop(int baseNr, 476 const short *const pair_table, 477 tBaseInformation *baseInformation, 478 const double unpaired, 479 const double paired) 480{ 481 int start = baseNr; 482 int end = pair_table[baseNr]; 483 484 int unpairedCount = 0; 485 int stemCount = 1; 486 487 /* count stems and unpaired bases to use for bulge detection */ 488 int i = start + 1; 489 490 while (i < end) { 491 if (pair_table[i] == 0) { 492 /* unpaired base */ 493 unpairedCount++; 494 i++; 495 } else if (pair_table[i] > i) { 496 /* found new stem */ 497 stemCount++; 498 i = pair_table[i]; 499 } else { 500 /* returned from stem */ 501 i++; 502 } 503 } 504 505 short isBulge = (stemCount == 2 && unpairedCount == 1); 506 if (isBulge) { 507 if (pair_table[start + 1] == 0) 508 /* unpaired on left strand */ 509 cfgGenHandleStem(start + 2, pair_table, baseInformation, unpaired, paired); 510 else 511 /* unpaired on the right strand */ 512 cfgGenHandleStem(start + 1, pair_table, baseInformation, unpaired, paired); 513 } else { 514 int m = stemCount; /* compare RNApuzzler.c -> f_handle_loop */ 515 int n = unpairedCount + stemCount; /* compare RNApuzzler.c -> f_handle_loop */ 516 double defaultRadius = approximateConfigArcRadius(paired, unpaired, m, n, MATH_TWO_PI); 517 config *cfgLoop = cfgGenerateDefaultConfig(pair_table, 518 start, 519 unpaired, 520 paired, 521 defaultRadius); 522 baseInformation[start].config = cfgLoop; 523 524 int i = start + 1; 525 while (i < end) { 526 if (pair_table[i] == 0) { 527 /* unpaired base */ 528 i++; 529 } else if (pair_table[i] > i) { 530 /* found new stem */ 531 cfgGenHandleStem(i, pair_table, baseInformation, unpaired, paired); 532 i = pair_table[i]; 533 } else { 534 /* returned from stem */ 535 i++; 536 } 537 } 538 } 539} 540 541 542/** 543 * @brief cfgGenHandleStem 544 * - recursively iterates through the RNA and generates default configurations. 545 * Alternates with corresponding handleLoop method. 546 * @param baseNr 547 * - index of the stem's first base 548 * @param pair_table 549 * - the RNA's pairing information 550 * @param baseInformation 551 * - array of tBaseInformation annotations (to save config) 552 */ 553PRIVATE void 554cfgGenHandleStem(int baseNr, 555 const short *const pair_table, 556 tBaseInformation *baseInformation, 557 const double unpaired, 558 const double paired) 559{ 560 /* iterating over the stem and calling cfgGenHandleLoop as soon as a loop is found */ 561 short continueStem = 1; 562 int i = baseNr; 563 564 while (continueStem) { 565 if (pair_table[i + 1] == pair_table[i] - 1) { 566 i++; 567 } else { 568 /* found unpaired above stem */ 569 cfgGenHandleLoop(i, pair_table, baseInformation, unpaired, paired); 570 continueStem = 0; 571 } 572 } 573} 574 575 576/* documentation at header file */ 577PRIVATE void 578cfgGenerateConfig(const short *const pair_table, 579 tBaseInformation *baseInformation, 580 const double unpaired, 581 const double paired) 582{ 583 /* 584 * iterate over RNA 585 * for each loop generate a default config 586 */ 587 int length = pair_table[0]; 588 int i = 1; 589 590 while (i < length) { 591 if (pair_table[i] == 0) { 592 /* unpaired at exterior loop */ 593 i++; 594 } else if (pair_table[i] > i) { 595 /* found stem */ 596 cfgGenHandleStem(i, pair_table, baseInformation, unpaired, paired); 597 i = pair_table[i]; 598 } else { 599 /* returned from stem */ 600 i++; 601 } 602 } 603} 604 605 606/* 607 *-------------------------------------------------------------------------- 608 *--- set and update config elements ----------------------------------- 609 *-------------------------------------------------------------------------- 610 */ 611 612/* documentation at header file */ 613/** 614 * @brief cfgSetRadius 615 * - changes the value of radius for a config to the given value 616 * @param config 617 * - config that is being altered 618 * @param radius 619 * - new radius 620 */ 621PRIVATE void 622cfgSetRadius(config *cfg, 623 const double radius) 624{ 625 cfg->radius = radius; 626} 627 628 629/** 630 * @brief cfgUpdateMinRadius 631 * - updates the minimum possible radius for the given config 632 * @param config 633 * @param unpaired 634 * @param paired 635 */ 636PRIVATE void 637cfgUpdateMinRadius(config *cfg, 638 const double unpaired, 639 const double paired) 640{ 641 double minRadius = approximateConfigRadius(cfg, unpaired, paired); 642 643 cfg->minRadius = minRadius; 644} 645 646 647PRIVATE double 648cfgApplyChanges(config *cfg, 649 const char loopName, 650 const double *deltaCfg, 651 const double radiusNew, 652 const vrna_plot_options_puzzler_t *puzzler) 653{ 654 /* start with adjusting config angles; if applicable */ 655 if (deltaCfg != NULL) { 656 for (int currentArc = 0; currentArc < cfg->numberOfArcs; currentArc++) 657 658 (cfg->cfgArcs[currentArc]).arcAngle += deltaCfg[currentArc]; 659 } 660 661 /* then, adjust config radius */ 662 double oldRadius = cfg->radius; 663 double newRadius = -1.0; 664 if (radiusNew > 0.0) { 665 /* 666 * in case the input is a positive value 667 * we set the minimum of valid and input as new radius 668 */ 669 cfgUpdateMinRadius(cfg, puzzler->unpaired, puzzler->paired); 670 newRadius = fmax(radiusNew, cfg->minRadius); 671 cfgSetRadius(cfg, newRadius); 672 } else if (radiusNew == 0.0) { 673 /* 674 * set the minRadius as new value 675 * (this allows to shrink a loop) 676 */ 677 cfgUpdateMinRadius(cfg, puzzler->unpaired, puzzler->paired); 678 newRadius = cfg->minRadius; 679 cfgSetRadius(cfg, newRadius); 680 } else if (radiusNew == -1.0) { 681 /* 682 * set the minRadius as new value 683 * (this forbidds to shrink a loop) 684 */ 685 cfgUpdateMinRadius(cfg, puzzler->unpaired, puzzler->paired); 686 687 if (cfg->minRadius - EPSILON_0 > oldRadius) { 688 newRadius = cfg->minRadius; 689 } else { 690 double defaultIncrease = 1.05; 691 newRadius = oldRadius * defaultIncrease; 692 } 693 694 cfgSetRadius(cfg, newRadius); 695 } else { 696 /* all unhandled inputs result in errors */ 697 newRadius = -1.0; 698 } 699 700 return newRadius; 701} 702 703 704PRIVATE short 705cfgIsValid(config *cfg, 706 const double *deltaCfg) 707{ 708 if (deltaCfg == NULL) 709 710 return 0; 711 712 double sumAngles = 0.0; 713 short validSingleAngles = 1; 714 715 for (int currentArc = 0; currentArc < cfg->numberOfArcs; currentArc++) { 716 double angle = getArcAngle(cfg, currentArc) + deltaCfg[currentArc]; 717 sumAngles += angle; 718 719 short validAngle = 0.0 < angle && angle < MATH_TWO_PI; 720 validSingleAngles = validSingleAngles && validAngle; 721 } 722 723 short validSumAngles = (fabs(sumAngles - MATH_TWO_PI) < EPSILON_3); 724 725 return validSingleAngles && validSumAngles; 726} 727 728 729#endif 730