1#ifndef RNAPUZZLER_COORDINATE_TRANSFORM 2#define RNAPUZZLER_COORDINATE_TRANSFORM 3 4/* 5 * RNAturtle algorithm 6 * 7 * c Daniel Wiegreffe, Daniel Alexander, Dirk Zeckzer 8 * ViennaRNA package 9 */ 10 11#include "../headers/tBaseInformation_struct.h" 12#include "vector_math.inc" 13#include "drawingconfig.inc" 14 15/** 16 * Compute affin angles for all loops 17 */ 18PRIVATE void 19computeAffineCoordinates(short const *const pair_table, 20 const double paired, 21 const double unpaired, 22 tBaseInformation *const baseInformation); 23 24 25/** 26 * Calculate the coordinates for the drawing with the given affin angles 27 */ 28PRIVATE void 29affineToCartesianCoordinates(tBaseInformation *const baseInformation, 30 unsigned short const length, 31 double *const x, 32 double *const y); 33 34 35/* -------------------------------------------------------------------------------------------------------------------------- */ 36 37/** 38 * Handle bases of the exterior loop 39 * 40 * @returns first paired base 41 */ 42PRIVATE short 43handleExteriorBases(const short *const pair_table, 44 short currentBase, 45 tBaseInformation *baseInformation, 46 const int direction) 47{ 48 short const length = pair_table[0]; 49 50 if (currentBase > 1) { 51 baseInformation[currentBase].angle = baseInformation[currentBase].angle + direction * 52 MATH_PI_HALF; 53 baseInformation[currentBase].baseType = TYPE_EXTERIOR; 54 } 55 56 while (currentBase < length && pair_table[currentBase] <= 0) { 57 baseInformation[currentBase + 1].angle = 0.0; 58 baseInformation[currentBase].baseType = TYPE_EXTERIOR; 59 currentBase++; 60 } 61 /* 62 * For following stem, since it cannot alter the first two bases 63 * which may contain loop exit angles. 64 */ 65 if ((currentBase + 1) <= (length)) { 66 baseInformation[currentBase + 1].angle = direction * MATH_PI_HALF; 67 baseInformation[currentBase].baseType = TYPE_EXTERIOR; 68 } else { 69 baseInformation[currentBase].baseType = TYPE_EXTERIOR; 70 } 71 72 return currentBase; 73} 74 75 76PRIVATE void 77countLoopPairs(short *const m, 78 short *const n, 79 short i, 80 const short *const pair_table) 81{ 82 short const end = pair_table[i++]; /* end is 1st base of 2nd stem half */ 83 84 *m = *n = 1; /* Base pair of enclosing stem, consecutive pair from last stem base to first loop base */ 85 while (i < end) { 86 if (pair_table[i] <= 0 || pair_table[i] < i) { 87 (*n)++; 88 i++; /* Go to consecutive base */ 89 } else { 90 /* base pair ahead */ 91 (*m)++; 92 i = pair_table[i]; /* Go to base pair partner */ 93 } 94 } 95 96 return; 97} 98 99 100/* -------------------------------------------------------------------------------------------------------------------------- */ 101PRIVATE void 102handleStem(const short *const pair_table, 103 short currentBase, 104 const double paired, 105 const double unpaired, 106 tBaseInformation *const baseInformation, 107 int direction); 108 109 110/* -------------------------------------------------------------------------------------------------------------------------- */ 111 112/*---------------------------------------------------------------------------*/ 113 114/** 115 * Detect bulges (internal loops with only one unpaired base) 116 */ 117PRIVATE int 118detectBulge(short start, 119 const short *const pair_table) 120{ 121 int bulge = 0; 122 int end = pair_table[start]; 123 int iterate = 1; 124 int old = 0; 125 int i = start + 1; 126 int side = 0; 127 128 do { 129 if (pair_table[i] > 0) { 130 /* There was a basepair before */ 131 if (iterate > 0) { 132 /* Stem entry */ 133 if (pair_table[i] == old) { 134 side++; 135 i++; 136 } /* Bulge found */ 137 else { 138 /* It's only a bulge if the paired base is the start base or the bulge is on the other strand */ 139 if (start == pair_table[i] || pair_table[i] == end - 2) { 140 bulge = pair_table[i]; 141 break; 142 } else { 143 break; 144 } 145 } 146 } /* Found first basepair */ 147 else { 148 iterate++; 149 old = i; 150 i = pair_table[i]; 151 } 152 } /* Consecutive Base found */ 153 else { 154 /* There was a basepair before */ 155 if (iterate > 0) { 156 iterate = 0; 157 i++; 158 } /* Just the next consecutive base */ 159 else { 160 i++; 161 } 162 } 163 } while (i > start); 164 165 return bulge; 166} 167 168 169PRIVATE void 170handleLoop(short i, 171 const short *const pair_table, 172 const double paired, 173 const double unpaired, 174 tBaseInformation *baseInformation, 175 int direction) 176{ 177 int start = i; 178 short const end = pair_table[i]; /* end is 1st base of 2nd stem half */ 179 short m, n; /* m - no. of base pairs; n - no. of consecutive pairs */ 180 181 countLoopPairs(&m, &n, i, pair_table); 182 183 184 int bulge = detectBulge(i, pair_table); 185 186 if (bulge > 0 && n - m == 1) { 187 /* detected bulge with a single unpaired base */ 188 int length = (unpaired * (n - m + 1)) / 2; 189 190 double alpha = acos(unpaired / (2 * length)); 191 192 if (pair_table[i + 1] == 0) { 193 baseInformation[i + 1].angle = baseInformation[i + 1].angle + direction * alpha; 194 baseInformation[i].baseType = TYPE_BULGE; 195 baseInformation[pair_table[i]].baseType = TYPE_BULGE; 196 i++; 197 198 baseInformation[i + 1].angle = -direction * alpha * 2; 199 baseInformation[i].baseType = TYPE_BULGE; 200 i++; 201 202 baseInformation[i + 1].angle = direction * alpha; 203 baseInformation[i].baseType = TYPE_BULGE; 204 baseInformation[pair_table[i]].baseType = TYPE_BULGE; 205 206 handleStem(pair_table, i, paired, unpaired, baseInformation, direction); 207 208 i = pair_table[i]; 209 } else { 210 baseInformation[i + 1].angle = baseInformation[i + 1].angle + 0.0; 211 baseInformation[i].baseType = TYPE_BULGE; 212 i++; 213 214 baseInformation[i + 1].angle = baseInformation[i + 1].angle + 0.0; 215 baseInformation[i + 1].baseType = TYPE_BULGE; 216 217 baseInformation[i + 2].angle = baseInformation[i + 2].angle + 0.0; 218 baseInformation[i + 1].baseType = TYPE_BULGE; 219 220 handleStem(pair_table, i, paired, unpaired, baseInformation, direction); 221 222 i = pair_table[i]; 223 224 baseInformation[i + 1].angle = baseInformation[i + 1].angle + direction * alpha; 225 baseInformation[i].baseType = TYPE_BULGE; 226 i++; 227 228 baseInformation[i + 1].angle = -direction * alpha * 2; 229 baseInformation[i].baseType = TYPE_BULGE; 230 i++; 231 232 baseInformation[i + 1].angle = direction * alpha; 233 baseInformation[i].baseType = TYPE_BULGE; 234 } 235 } else { 236 /* loop is not bulge with a single unpaired base */ 237 238 /* 239 * ****************************************************************** 240 * ***************** Loop Drawing Algorithm - Start ***************** 241 * ****************************************************************** 242 */ 243 244 245 /* variables for config drawing */ 246 double angle_over_paired; /* alpha at calculation */ 247 double current_angle; /* phi at calculation */ 248 double current_bb_angle; /* beta at calculation */ 249 double current_distance; /* b at calculation */ 250 double current_delta_ab; /* delta_ab at calculation */ 251 double current_delta_bb; /* delta_bb at calculation */ 252 int current_stem_count; 253 254 config *cfg = baseInformation[start].config; 255 int currentArc = 0; 256 double r = cfg->radius; 257 258 angle_over_paired = 2 * asin(paired / (2 * r)); 259 260 current_angle = getArcAngle(cfg, currentArc); 261 current_bb_angle = (current_angle - angle_over_paired) / 262 (cfg->cfgArcs[currentArc]).numberOfArcSegments; 263 current_distance = sqrt(2 * r * r * (1 - cos(current_bb_angle))); 264 current_delta_ab = 0.5 * (MATH_PI + angle_over_paired + current_bb_angle); 265 current_delta_bb = MATH_PI + current_bb_angle; 266 ++currentArc; 267 268 /* consider the direction that was given before (e.g. because of a bulge) */ 269 baseInformation[i + 1].angle = baseInformation[i + 1].angle + direction * 270 (MATH_PI - current_delta_ab); 271 272 baseInformation[i].distance = current_distance; 273 274 current_stem_count = 0; 275 276 /* Double Loop Check: Loop is followed directly by a loop */ 277 if (baseInformation[i].baseType == TYPE_LOOP1) 278 baseInformation[i].baseType = TYPE_LOOP2; 279 else 280 baseInformation[i].baseType = TYPE_LOOP1; 281 282 i++; 283 284 while (i < end) { 285 if (pair_table[i] <= 0) { 286 /* unpaired */ 287 288 /* handle backbone element */ 289 baseInformation[i + 1].angle = -direction * (current_delta_bb - MATH_PI); 290 baseInformation[i].distance = current_distance; 291 292 baseInformation[i].baseType = TYPE_LOOP1; 293 i++; 294 } else if (pair_table[i] > i) { 295 /* base pair ahead */ 296 297 /* i is the beginning of a stem and therefore the end of the current arc */ 298 baseInformation[i + 1].angle = direction * (MATH_PI - current_delta_ab); 299 current_stem_count++; 300 baseInformation[i].baseType = TYPE_LOOP1; 301 handleStem(pair_table, i, paired, unpaired, baseInformation, direction); /* Recurse multiloop */ 302 i = pair_table[i]; /* Go to base pair partner */ 303 } else { 304 /* base pair behind */ 305 306 307 /* 308 * i is the end of a stem that returned to the loop 309 * recalculation if angles and distances is nessecary because the next arc has been reached 310 */ 311 if (current_stem_count == 1) { 312 current_stem_count = 0; 313 current_angle = getArcAngle(cfg, currentArc); 314 current_bb_angle = (current_angle - angle_over_paired) / 315 (cfg->cfgArcs[currentArc]).numberOfArcSegments; 316 current_distance = sqrt(2 * r * r * (1 - cos(current_bb_angle))); 317 current_delta_ab = 0.5 * (MATH_PI + angle_over_paired + current_bb_angle); 318 current_delta_bb = MATH_PI + current_bb_angle; 319 ++currentArc; 320 } 321 322 /* consider the direction that was given before (e.g. because of a bulge) */ 323 baseInformation[i + 1].angle = baseInformation[i + 1].angle + direction * 324 (MATH_PI - current_delta_ab); 325 baseInformation[i].distance = current_distance; 326 baseInformation[i].baseType = TYPE_LOOP1; 327 328 i++; 329 } 330 } 331 332 if ((i + 1) <= pair_table[0]) 333 baseInformation[i + 1].angle = direction * (MATH_PI - current_delta_ab); 334 335 baseInformation[i].baseType = TYPE_LOOP1; 336 337 338 /* 339 * ****************************************************************** 340 * ****************** Loop Drawing Algorithm - End ****************** 341 * ****************************************************************** 342 */ 343 } 344 345 return; 346} 347 348 349/* -------------------------------------------------------------------------------------------------------------------------- */ 350 351PRIVATE void 352handleStem(const short *const pair_table, 353 short i, 354 const double paired, 355 const double unpaired, 356 tBaseInformation *const baseInformation, 357 int direction) 358{ 359 short end = pair_table[i] + 1; /* first base after stem */ 360 361 /* 362 * Skip first two bases of stem, angle already written by 363 * either loop (+exit angle) or handle_unpaired 364 * i+=2; 365 */ 366 baseInformation[i].baseType = TYPE_STEM; 367 i++; 368 369 while (pair_table[i] > 0 && /* i is base paired */ 370 (pair_table[i] == end - 1 || /* First position of stem, continue! */ 371 pair_table[i] + 1 == pair_table[i - 1] /* Detect bulges on opposite strand */ 372 ) 373 ) { 374 baseInformation[i + 1].angle = 0.0; 375 baseInformation[i].baseType = TYPE_STEM; 376 i++; 377 } 378 if (pair_table[i] == end - 1) { 379 } else { 380 handleLoop(--i, pair_table, paired, unpaired, baseInformation, direction); /* i - last base of first stem half */ 381 } 382 383 i = pair_table[i]; /* set i as base pairing partner of last stem base */ 384 baseInformation[i].baseType = TYPE_STEM; 385 i++; 386 387 while (i < end && i < pair_table[0]) { 388 baseInformation[i].baseType = TYPE_STEM; 389 i++; 390 } 391 392 return; 393} 394 395 396/* ------------------------------------------------------------------------------ */ 397 398/** 399 * Compute angle angles for all loops 400 */ 401PRIVATE void 402computeAffineCoordinates(short const *const pair_table, 403 const double paired, 404 const double unpaired, 405 tBaseInformation *const baseInformation) 406{ 407 /* Initialization */ 408 short const length = pair_table[0]; 409 short currentBase = 1; 410 const int direction = -1; 411 412 baseInformation[0].angle = 0.0; 413 414 /* 415 * For first stem, since handle_stem cannot alter the first two bases 416 * which may contain loop exit angles (in other stems). 417 */ 418 if (2 <= length) { 419 baseInformation[1].angle = baseInformation[0].angle; 420 baseInformation[2].angle = baseInformation[1].angle; 421 } 422 423 /* reset dangle_count */ 424 int dangle_count = 0; 425 while (currentBase < length) { 426 if (pair_table[currentBase] <= 0) { 427 if (currentBase > 1) 428 429 baseInformation[currentBase - 1].baseType = TYPE_EXTERIOR; 430 431 currentBase = handleExteriorBases(pair_table, currentBase, baseInformation, direction); 432 /* returns first paired base currentBase */ 433 dangle_count++; 434 } 435 436 if (currentBase < length) { 437 /* it was not the dangling end */ 438 if (pair_table[currentBase] - pair_table[currentBase - 1] != 1 439 && pair_table[currentBase] != 0 440 && pair_table[currentBase - 1] != 0) { 441 if (currentBase == 1) { 442 if (dangle_count < 1) { 443 baseInformation[0].angle = 444 baseInformation[1].angle = 445 baseInformation[2].angle = -MATH_PI_HALF; 446 baseInformation[currentBase].baseType = TYPE_EXTERIOR; 447 } 448 449 handleStem(pair_table, currentBase, paired, unpaired, baseInformation, direction); 450 currentBase = pair_table[currentBase] + 1; 451 452 /* Check, if there is a minor dangling end at the end */ 453 if (currentBase == length) { 454 baseInformation[currentBase - 1].baseType = TYPE_EXTERIOR; 455 baseInformation[currentBase].baseType = TYPE_EXTERIOR; 456 baseInformation[currentBase].angle = -MATH_PI_HALF; 457 } 458 459 continue; 460 } else { 461 baseInformation[currentBase].angle = baseInformation[currentBase].angle + 462 direction * MATH_PI_HALF; 463 baseInformation[currentBase + 1].distance = unpaired; 464 baseInformation[currentBase - 1].baseType = TYPE_EXTERIOR; 465 baseInformation[currentBase + 1].angle = baseInformation[currentBase + 1].angle + 466 direction * MATH_PI_HALF; 467 baseInformation[currentBase].baseType = TYPE_EXTERIOR; 468 469 dangle_count++; 470 } 471 } 472 473 handleStem(pair_table, currentBase, paired, unpaired, baseInformation, direction); 474 475 currentBase = pair_table[currentBase] + 1; /* currentBase is next base after stem */ 476 477 if (currentBase == length) { 478 baseInformation[currentBase - 1].baseType = TYPE_EXTERIOR; 479 currentBase = handleExteriorBases(pair_table, 480 currentBase, 481 baseInformation, 482 direction); 483 } 484 } 485 } 486 baseInformation[length].baseType = TYPE_EXTERIOR; 487} 488 489 490/* ------------------------------------------------------------------------------ */ 491 492/** 493 * Calculate the coordinates for the drawing with the given angle angles 494 */ 495PRIVATE void 496affineToCartesianCoordinates(tBaseInformation *const baseInformation, 497 unsigned short const length, 498 double *const x, 499 double *const y) 500{ 501 if (length < 1) 502 return; 503 504 double angle = 0.0; 505 x[0] = y[0] = EXTERIOR_Y; 506 for (int i = 1; i < length; i++) { 507 angle = angle - baseInformation[i + 1].angle; 508 509 x[i] = x[i - 1] + baseInformation[i].distance * cos(angle); 510 y[i] = y[i - 1] + baseInformation[i].distance * sin(angle); 511 } 512 return; 513} 514 515 516#endif 517