1 #ifndef _RNA_ALGEBRA_H_ 2 #define _RNA_ALGEBRA_H_ 3 4 #include <assert.h> 5 #include <algorithm> 6 #include <climits> 7 8 #include "algebra.h" 9 #include "debug.h" 10 #include "misc.h" 11 #include "rna_alphabet.h" 12 #include "rnaforester_options.h" 13 14 // TODO in und del sind immer gleich, also kann man das in private indel fknen zusammenfassen 15 // TODO mit einer ordentlichen hierarchie koennte man z.b. den score in der oberklasse setzen 16 // man koennte dann auch einfach alle klassen in diesem file beerben, und um die affinen funktionen erweitern 17 // TODO zerlegen in verschiedene files, fuer profile und szalgebra extra 18 19 // Definitions and typedefs 20 21 const double DBL_NEG = -100000000.0; // TODO the values from limits.h caused problems.. 22 const double DBL_POS = 100000000.0; 23 24 // Class Score 25 26 // Class Score reads scoring parameters from RNAforester command line. 27 class Score { 28 friend std::ostream& operator<< (std::ostream &out, const Score &score); 29 private: 30 bool isAffine_; 31 bool isLocal_; 32 bool isDistance_; 33 bool isRIBOSUM_; 34 35 public: 36 double bp_rep_score_; 37 double bp_indel_score_; 38 double bp_indel_open_score_; 39 double b_match_score_; 40 double b_rep_score_; 41 double b_indel_score_; 42 double b_indel_open_score_; 43 44 Score(Options &options); 45 Score(const Score &s); 46 }; 47 48 49 // RNA Algebra Classes 50 51 52 // Similarity algebra for RNA forests 53 54 class DoubleSimiRNA_Algebra : public RNA_Algebra<double,RNA_Alphabet> { 55 private: 56 Score s_; 57 58 public: empty()59 double empty() const { 60 return 0; 61 }; 62 replacepair(RNA_Alphabet la,RNA_Alphabet lb,double down,RNA_Alphabet ra,RNA_Alphabet rb,double over)63 double replacepair(RNA_Alphabet la, RNA_Alphabet lb, double down, RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 64 return s_.bp_rep_score_+down+over; 65 }; 66 replace(RNA_Alphabet a,double down,RNA_Alphabet b,double over)67 double replace(RNA_Alphabet a,double down, RNA_Alphabet b, double over) const { 68 if (a==ALPHA_BASEPAIR && b == ALPHA_BASEPAIR) 69 return s_.bp_rep_score_+down+over; 70 else { 71 if (a==ALPHA_BASEPAIR || b==ALPHA_BASEPAIR) 72 return INT_MIN; 73 else { 74 if (a==b) 75 return s_.b_match_score_+down+over; 76 else 77 return s_.b_rep_score_+down+over; 78 } 79 } 80 }; 81 del(RNA_Alphabet a,double down,double over)82 double del(RNA_Alphabet a,double down, double over) const { 83 if (a==ALPHA_BASEPAIR) 84 return s_.bp_indel_score_+down+over; 85 else 86 return s_.b_indel_score_+down+over; 87 }; 88 insert(double down,RNA_Alphabet b,double over)89 double insert(double down,RNA_Alphabet b,double over) const { 90 if (b==ALPHA_BASEPAIR) 91 return s_.bp_indel_score_+down+over; 92 else 93 return s_.b_indel_score_+down+over; 94 }; 95 choice(double a,double b)96 double choice(double a, double b) const { 97 return std::max(a,b); 98 }; 99 worst_score()100 double worst_score() const { 101 return INT_MIN; 102 }; 103 DoubleSimiRNA_Algebra(const Score & s)104 DoubleSimiRNA_Algebra(const Score &s) 105 : s_(s) {}; 106 107 // Pair indels 108 109 // three/four possibilities for scores, choose via choice function: 110 // delete pairing only, but leave bases 111 // del p + rep l + mdown + rep r + over 112 // delete pairing and both bases 113 // del p + del l + mdown + del r + over 114 // delete pairing and one base 115 // del p + del l + mdown + rep r + over 116 // del p + rep l + mdown + del r + over deletePairOnly(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)117 double deletePairOnly(RNA_Alphabet la, RNA_Alphabet lb, double mdown, RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 118 // return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_rep_score_ + over; 119 double result = s_.bp_indel_score_; 120 if (la == lb) 121 result += s_.b_match_score_; 122 else 123 result += s_.b_rep_score_; 124 result += mdown; 125 if (ra == rb) 126 result += s_.b_match_score_; 127 else 128 result += s_.b_rep_score_; 129 result += over; 130 return result; 131 } deletePairAndBases(RNA_Alphabet la,double mdown,RNA_Alphabet ra,double over)132 double deletePairAndBases(RNA_Alphabet la, double mdown, RNA_Alphabet ra, double over) const { 133 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_indel_score_ + over; 134 double result = s_.bp_indel_score_; 135 result += s_.b_indel_score_; 136 result += mdown; 137 result += s_.b_indel_score_; 138 result += over; 139 return result; 140 } deletePairAndLeftBase(RNA_Alphabet la,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)141 double deletePairAndLeftBase(RNA_Alphabet la, double mdown, RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 142 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_rep_score_ + over; 143 double result = s_.bp_indel_score_; 144 result += s_.b_indel_score_; 145 result += mdown; 146 if (ra == rb) 147 result += s_.b_match_score_; 148 else 149 result += s_.b_rep_score_; 150 result += over; 151 return result; 152 } deletePairAndRightBase(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet ra,double over)153 double deletePairAndRightBase(RNA_Alphabet la, RNA_Alphabet lb, double mdown, RNA_Alphabet ra, double over) const { 154 // return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_indel_score_ + over; 155 double result = s_.bp_indel_score_; 156 if (la == lb) 157 result += s_.b_match_score_; 158 else 159 result += s_.b_rep_score_; 160 result += mdown; 161 result += s_.b_indel_score_; 162 result += over; 163 return result; 164 } 165 166 167 // three/four possibilities for scores, choose via choice function: 168 // insert pairing only, but leave bases 169 // ins p + rep l + mdown + rep r + over 170 // insert pairing and both bases 171 // ins p + ins l + mdown + ins r + over 172 // insert pairing and one base 173 // ins p + ins l + mdown + rep r + over 174 // ins p + rep l + mdown + ins r + over insertPairOnly(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)175 double insertPairOnly(RNA_Alphabet la, RNA_Alphabet lb, double mdown, RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 176 // return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_rep_score_ + over; 177 double result = s_.b_indel_score_; 178 if (la == lb) { 179 result += s_.b_match_score_; 180 } 181 else { 182 result += s_.b_rep_score_; 183 } 184 result += mdown; 185 if (ra == rb) { 186 result += s_.b_match_score_; 187 } 188 else { 189 result += s_.b_rep_score_; 190 } 191 result += over; 192 return result; 193 } insertPairAndBases(RNA_Alphabet lb,double mdown,RNA_Alphabet rb,double over)194 double insertPairAndBases(RNA_Alphabet lb, double mdown, RNA_Alphabet rb, double over) const { 195 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_indel_score_ + over; 196 double result = s_.bp_indel_score_; 197 result += s_.b_indel_score_; 198 result += mdown; 199 result += s_.b_indel_score_; 200 result += over; 201 return result; 202 } insertPairAndLeftBase(RNA_Alphabet lb,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)203 double insertPairAndLeftBase(RNA_Alphabet lb, double mdown, RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 204 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_rep_score_ + over; 205 double result = s_.bp_indel_score_; 206 result += s_.b_indel_score_; 207 result += mdown; 208 if (ra == rb) 209 result += s_.b_match_score_; 210 else 211 result += s_.b_rep_score_; 212 result += over; 213 return result; 214 } insertPairAndRightBase(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet rb,double over)215 double insertPairAndRightBase(RNA_Alphabet la, RNA_Alphabet lb, double mdown, RNA_Alphabet rb, double over) const { 216 // return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_indel_score_ + over; 217 double result = s_.bp_indel_score_; 218 if (la == lb) 219 result += s_.b_match_score_; 220 else 221 result += s_.b_rep_score_; 222 result += mdown; 223 result += s_.b_indel_score_; 224 result += over; 225 return result; 226 } 227 228 }; 229 230 class AffineDoubleSimiRNA_Algebra : public RNA_AlgebraAffine<double,RNA_Alphabet> { 231 232 private: 233 Score s_; 234 public: empty()235 double empty() const { 236 return 0; 237 }; replacepair(RNA_Alphabet la,RNA_Alphabet lb,double down,RNA_Alphabet ra,RNA_Alphabet rb,double over)238 double replacepair(RNA_Alphabet la, RNA_Alphabet lb, double down, RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 239 return s_.bp_rep_score_+down+over; 240 }; 241 replace(RNA_Alphabet a,double down,RNA_Alphabet b,double over)242 double replace(RNA_Alphabet a,double down, RNA_Alphabet b, double over) const { 243 if (a==ALPHA_BASEPAIR && b == ALPHA_BASEPAIR) 244 return s_.bp_rep_score_+down+over; 245 else { 246 if (a==ALPHA_BASEPAIR || b==ALPHA_BASEPAIR) 247 return INT_MIN; 248 else { 249 if (a==b) 250 return s_.b_match_score_+down+over; 251 else 252 return s_.b_rep_score_+down+over; 253 } 254 } 255 }; 256 del(RNA_Alphabet a,double down,double over)257 double del(RNA_Alphabet a,double down, double over) const { 258 if (a==ALPHA_BASEPAIR) 259 return s_.bp_indel_score_+down+over; 260 else 261 return s_.b_indel_score_+down+over; 262 }; 263 insert(double down,RNA_Alphabet b,double over)264 double insert(double down,RNA_Alphabet b,double over) const { 265 if (b==ALPHA_BASEPAIR) 266 return s_.bp_indel_score_+down+over; 267 else 268 return s_.b_indel_score_+down+over; 269 }; 270 choice(double a,double b)271 double choice(double a, double b) const { 272 return std::max(a,b); 273 }; 274 worst_score()275 double worst_score() const { 276 return INT_MIN; 277 }; 278 279 delO(RNA_Alphabet a,double down,double over)280 double delO(RNA_Alphabet a,double down, double over) const { 281 if (a==ALPHA_BASEPAIR) 282 return s_.bp_indel_open_score_+down+over; 283 else 284 return s_.b_indel_open_score_+down+over; 285 }; 286 insertO(double down,RNA_Alphabet b,double over)287 double insertO(double down,RNA_Alphabet b,double over) const { 288 if (b==ALPHA_BASEPAIR) 289 return s_.bp_indel_open_score_+down+over; 290 else 291 return s_.b_indel_open_score_+down+over; 292 }; 293 AffineDoubleSimiRNA_Algebra(const Score & s)294 AffineDoubleSimiRNA_Algebra(const Score &s) 295 : s_(s) {}; 296 297 // Pair indels 298 299 300 // three/four possibilities for scores, choose via choice function: 301 // delete pairing only, but leave bases 302 // del p + rep l + mdown + rep r + over 303 // delete pairing and both bases 304 // del p + del l + mdown + del r + over 305 // delete pairing and one base 306 // del p + del l + mdown + rep r + over 307 // del p + rep l + mdown + del r + over deletePairOnly(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)308 double deletePairOnly(RNA_Alphabet la, RNA_Alphabet lb, double mdown, 309 RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 310 double result = s_.bp_indel_score_; 311 if (la == lb) 312 result += s_.b_match_score_; 313 else 314 result += s_.b_rep_score_; 315 result += mdown; 316 if (ra == rb) 317 result += s_.b_match_score_; 318 else 319 result += s_.b_rep_score_; 320 result += over; 321 return result; 322 // return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_rep_score_ + over; 323 } 324 deletePairAndBases(RNA_Alphabet la,double mdown,RNA_Alphabet ra,double over)325 double deletePairAndBases(RNA_Alphabet la, double mdown, 326 RNA_Alphabet ra, double over) const { 327 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_indel_score_ + over; 328 double result = s_.bp_indel_score_; 329 result += s_.b_indel_score_; 330 result += mdown; 331 result += s_.b_indel_score_; 332 result += over; 333 return result; 334 } 335 deletePairAndLeftBase(RNA_Alphabet la,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)336 double deletePairAndLeftBase(RNA_Alphabet la, double mdown, 337 RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 338 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_rep_score_ + over; 339 double result = s_.bp_indel_score_; 340 result += s_.b_indel_score_; 341 result += mdown; 342 if (ra == rb) 343 result += s_.b_match_score_; 344 else 345 result += s_.b_rep_score_; 346 result += over; 347 return result; 348 349 } 350 deletePairAndRightBase(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet ra,double over)351 double deletePairAndRightBase(RNA_Alphabet la, RNA_Alphabet lb, double mdown, 352 RNA_Alphabet ra, double over) const { 353 // return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_indel_score_ + over; 354 double result = s_.bp_indel_score_; 355 if (la == lb) 356 result += s_.b_match_score_; 357 else 358 result += s_.b_rep_score_; 359 result += mdown; 360 result += s_.b_indel_score_; 361 result += over; 362 return result; 363 } 364 365 366 // three/four possibilities for scores, choose via choice function: 367 // insert pairing only, but leave bases 368 // ins p + rep l + mdown + rep r + over 369 // insert pairing and both bases 370 // ins p + ins l + mdown + ins r + over 371 // insert pairing and one base 372 // ins p + ins l + mdown + rep r + over 373 // ins p + rep l + mdown + ins r + over insertPairOnly(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)374 double insertPairOnly(RNA_Alphabet la, RNA_Alphabet lb, double mdown, 375 RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 376 // return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_rep_score_ + over; 377 double result = s_.b_indel_score_; 378 if (la == lb) { 379 result += s_.b_match_score_; 380 } 381 else { 382 result += s_.b_rep_score_; 383 } 384 result += mdown; 385 if (ra == rb) { 386 result += s_.b_match_score_; 387 } 388 else { 389 result += s_.b_rep_score_; 390 } 391 result += over; 392 return result; 393 } insertPairAndBases(RNA_Alphabet lb,double mdown,RNA_Alphabet rb,double over)394 double insertPairAndBases(RNA_Alphabet lb, double mdown, 395 RNA_Alphabet rb, double over) const { 396 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_indel_score_ + over; 397 double result = s_.bp_indel_score_; 398 result += s_.b_indel_score_; 399 result += mdown; 400 result += s_.b_indel_score_; 401 result += over; 402 return result; 403 } insertPairAndLeftBase(RNA_Alphabet lb,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)404 double insertPairAndLeftBase(RNA_Alphabet lb, double mdown, 405 RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 406 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_rep_score_ + over; 407 double result = s_.bp_indel_score_; 408 result += s_.b_indel_score_; 409 result += mdown; 410 if (ra == rb) 411 result += s_.b_match_score_; 412 else 413 result += s_.b_rep_score_; 414 result += over; 415 return result; 416 } insertPairAndRightBase(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet rb,double over)417 double insertPairAndRightBase(RNA_Alphabet la, RNA_Alphabet lb, double mdown, 418 RNA_Alphabet rb, double over) const { 419 // return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_indel_score_ + over; 420 double result = s_.bp_indel_score_; 421 if (la == lb) 422 result += s_.b_match_score_; 423 else 424 result += s_.b_rep_score_; 425 result += mdown; 426 result += s_.b_indel_score_; 427 result += over; 428 return result; 429 } 430 431 }; 432 433 // Distance algebra for RNA forests 434 class DoubleDistRNA_Algebra : public RNA_Algebra<double,RNA_Alphabet> { 435 private: 436 Score s_; 437 438 public: empty()439 double empty() const { 440 return 0; 441 }; replacepair(RNA_Alphabet la,RNA_Alphabet lb,double down,RNA_Alphabet ra,RNA_Alphabet rb,double over)442 double replacepair(RNA_Alphabet la, RNA_Alphabet lb, double down, RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 443 return s_.bp_rep_score_+down+over; 444 }; 445 replace(RNA_Alphabet a,double down,RNA_Alphabet b,double over)446 double replace(RNA_Alphabet a,double down, RNA_Alphabet b, double over) const { 447 if (a==ALPHA_BASEPAIR && b == ALPHA_BASEPAIR) 448 return s_.bp_rep_score_+down+over; 449 else { 450 if (a==ALPHA_BASEPAIR || b==ALPHA_BASEPAIR) 451 return INT_MAX; 452 else { 453 if (a==b) 454 return s_.b_match_score_+down+over; 455 else 456 return s_.b_rep_score_+down+over; 457 } 458 } 459 }; 460 del(RNA_Alphabet a,double down,double over)461 double del(RNA_Alphabet a,double down,double over) const { 462 if (a==ALPHA_BASEPAIR) 463 return s_.bp_indel_score_+down+over; 464 else 465 return s_.b_indel_score_+down+over; 466 }; 467 insert(double down,RNA_Alphabet b,double over)468 double insert(double down,RNA_Alphabet b,double over) const { 469 if (b==ALPHA_BASEPAIR) 470 return s_.bp_indel_score_+down+over; 471 else 472 return s_.b_indel_score_+down+over; 473 }; 474 choice(double a,double b)475 double choice(double a, double b) const { 476 return std::min(a,b); 477 }; 478 worst_score()479 double worst_score() const { 480 return INT_MAX; 481 }; 482 DoubleDistRNA_Algebra(const Score & s)483 DoubleDistRNA_Algebra(const Score &s) 484 : s_(s) {}; 485 486 // Pair indel 487 // three/four possibilities for scores, choose via choice function: 488 // delete pairing only, but leave bases 489 // del p + rep l + mdown + rep r + over 490 // delete pairing and both bases 491 // del p + del l + mdown + del r + over 492 // delete pairing and one base 493 // del p + del l + mdown + rep r + over 494 // del p + rep l + mdown + del r + over deletePairOnly(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)495 double deletePairOnly(RNA_Alphabet la, RNA_Alphabet lb, double mdown, RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 496 // return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_rep_score_ + over; 497 double result = s_.bp_indel_score_; 498 if (la == lb) 499 result += s_.b_match_score_; 500 else 501 result += s_.b_rep_score_; 502 result += mdown; 503 if (ra == rb) 504 result += s_.b_match_score_; 505 else 506 result += s_.b_rep_score_; 507 result += over; 508 return result; 509 } deletePairAndBases(RNA_Alphabet la,double mdown,RNA_Alphabet ra,double over)510 double deletePairAndBases(RNA_Alphabet la, double mdown, RNA_Alphabet ra, double over) const { 511 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_indel_score_ + over; 512 double result = s_.bp_indel_score_; 513 result += s_.b_indel_score_; 514 result += mdown; 515 result += s_.b_indel_score_; 516 result += over; 517 return result; 518 } deletePairAndLeftBase(RNA_Alphabet la,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)519 double deletePairAndLeftBase(RNA_Alphabet la, double mdown, RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 520 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_rep_score_ + over; 521 double result = s_.bp_indel_score_; 522 result += s_.b_indel_score_; 523 result += mdown; 524 if (ra == rb) 525 result += s_.b_match_score_; 526 else 527 result += s_.b_rep_score_; 528 result += over; 529 return result; 530 } deletePairAndRightBase(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet ra,double over)531 double deletePairAndRightBase(RNA_Alphabet la, RNA_Alphabet lb, double mdown, RNA_Alphabet ra, double over) const { 532 // return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_indel_score_ + over; 533 double result = s_.bp_indel_score_; 534 if (la == lb) 535 result += s_.b_match_score_; 536 else 537 result += s_.b_rep_score_; 538 result += mdown; 539 result += s_.b_indel_score_; 540 result += over; 541 return result; 542 543 } 544 545 // three/four possibilities for scores, choose via choice function: 546 // insert pairing only, but leave bases 547 // ins p + rep l + mdown + rep r + over 548 // insert pairing and both bases 549 // ins p + ins l + mdown + ins r + over 550 // insert pairing and one base 551 // ins p + ins l + mdown + rep r + over 552 // ins p + rep l + mdown + ins r + over insertPairOnly(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)553 double insertPairOnly(RNA_Alphabet la, RNA_Alphabet lb, double mdown, RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 554 return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_rep_score_ + over; 555 double result = s_.b_indel_score_; 556 if (la == lb) { 557 result += s_.b_match_score_; 558 } 559 else { 560 result += s_.b_rep_score_; 561 } 562 result += mdown; 563 if (ra == rb) { 564 result += s_.b_match_score_; 565 } 566 else { 567 result += s_.b_rep_score_; 568 } 569 result += over; 570 return result; 571 } insertPairAndBases(RNA_Alphabet lb,double mdown,RNA_Alphabet rb,double over)572 double insertPairAndBases(RNA_Alphabet lb, double mdown, RNA_Alphabet rb, double over) const { 573 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_indel_score_ + over; 574 double result = s_.bp_indel_score_; 575 result += s_.b_indel_score_; 576 result += mdown; 577 result += s_.b_indel_score_; 578 result += over; 579 return result; 580 } insertPairAndLeftBase(RNA_Alphabet lb,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)581 double insertPairAndLeftBase(RNA_Alphabet lb, double mdown, RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 582 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_rep_score_ + over; 583 double result = s_.bp_indel_score_; 584 result += s_.b_indel_score_; 585 result += mdown; 586 if (ra == rb) 587 result += s_.b_match_score_; 588 else 589 result += s_.b_rep_score_; 590 result += over; 591 return result; 592 } insertPairAndRightBase(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet rb,double over)593 double insertPairAndRightBase(RNA_Alphabet la, RNA_Alphabet lb, double mdown, RNA_Alphabet rb, double over) const { 594 // return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_indel_score_ + over; 595 double result = s_.bp_indel_score_; 596 if (la == lb) 597 result += s_.b_match_score_; 598 else 599 result += s_.b_rep_score_; 600 result += mdown; 601 result += s_.b_indel_score_; 602 result += over; 603 return result; 604 } 605 606 607 }; 608 609 // Distance algebra for RNA forests 610 class AffineDoubleDistRNA_Algebra : public RNA_AlgebraAffine<double,RNA_Alphabet> { 611 private: 612 Score s_; 613 public: empty()614 double empty() const { 615 return 0; 616 }; replacepair(RNA_Alphabet la,RNA_Alphabet lb,double down,RNA_Alphabet ra,RNA_Alphabet rb,double over)617 double replacepair(RNA_Alphabet la, RNA_Alphabet lb, double down, RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 618 return s_.bp_rep_score_+down+over; 619 }; 620 replace(RNA_Alphabet a,double down,RNA_Alphabet b,double over)621 double replace(RNA_Alphabet a,double down, RNA_Alphabet b, double over) const { 622 if (a==ALPHA_BASEPAIR && b == ALPHA_BASEPAIR) 623 return s_.bp_rep_score_+down+over; 624 else { 625 if (a==ALPHA_BASEPAIR || b==ALPHA_BASEPAIR) 626 return INT_MAX; 627 else { 628 if (a==b) 629 return s_.b_match_score_+down+over; 630 else 631 return s_.b_rep_score_+down+over; 632 } 633 } 634 }; 635 del(RNA_Alphabet a,double down,double over)636 double del(RNA_Alphabet a,double down,double over) const { 637 if (a==ALPHA_BASEPAIR) 638 return s_.bp_indel_score_+down+over; 639 else 640 return s_.b_indel_score_+down+over; 641 }; 642 insert(double down,RNA_Alphabet b,double over)643 double insert(double down,RNA_Alphabet b,double over) const { 644 if (b==ALPHA_BASEPAIR) 645 return s_.bp_indel_score_+down+over; 646 else 647 return s_.b_indel_score_+down+over; 648 }; 649 choice(double a,double b)650 double choice(double a, double b) const { 651 return std::min(a,b); 652 }; 653 worst_score()654 double worst_score() const { 655 return INT_MAX; 656 }; 657 658 delO(RNA_Alphabet a,double down,double over)659 double delO(RNA_Alphabet a,double down,double over) const { 660 if (a==ALPHA_BASEPAIR) 661 return s_.bp_indel_open_score_+down+over; 662 else 663 return s_.b_indel_open_score_+down+over; 664 }; 665 insertO(double down,RNA_Alphabet b,double over)666 double insertO(double down,RNA_Alphabet b,double over) const { 667 if (b==ALPHA_BASEPAIR) 668 return s_.bp_indel_open_score_+down+over; 669 else 670 return s_.b_indel_open_score_+down+over; 671 }; 672 AffineDoubleDistRNA_Algebra(const Score & s)673 AffineDoubleDistRNA_Algebra(const Score &s) 674 : s_(s) {}; 675 676 // three/four possibilities for scores, choose via choice function: 677 // delete pairing only, but leave bases 678 // del p + rep l + mdown + rep r + over 679 // delete pairing and both bases 680 // del p + del l + mdown + del r + over 681 // delete pairing and one base 682 // del p + del l + mdown + rep r + over 683 // del p + rep l + mdown + del r + over deletePairOnly(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)684 double deletePairOnly(RNA_Alphabet la, RNA_Alphabet lb, double mdown, 685 RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 686 // return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_rep_score_ + over; 687 double result = s_.bp_indel_score_; 688 if (la == lb) 689 result += s_.b_match_score_; 690 else 691 result += s_.b_rep_score_; 692 result += mdown; 693 if (ra == rb) 694 result += s_.b_match_score_; 695 else 696 result += s_.b_rep_score_; 697 result += over; 698 return result; 699 } deletePairAndBases(RNA_Alphabet la,double mdown,RNA_Alphabet ra,double over)700 double deletePairAndBases(RNA_Alphabet la, double mdown, 701 RNA_Alphabet ra, double over) const { 702 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_indel_score_ + over; 703 double result = s_.bp_indel_score_; 704 result += s_.b_indel_score_; 705 result += mdown; 706 result += s_.b_indel_score_; 707 result += over; 708 return result; 709 } deletePairAndLeftBase(RNA_Alphabet la,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)710 double deletePairAndLeftBase(RNA_Alphabet la, double mdown, 711 RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 712 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_rep_score_ + over; 713 double result = s_.bp_indel_score_; 714 result += s_.b_indel_score_; 715 result += mdown; 716 if (ra == rb) 717 result += s_.b_match_score_; 718 else 719 result += s_.b_rep_score_; 720 result += over; 721 return result; 722 } deletePairAndRightBase(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet ra,double over)723 double deletePairAndRightBase(RNA_Alphabet la, RNA_Alphabet lb, double mdown, 724 RNA_Alphabet ra, double over) const { 725 // return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_indel_score_ + over; 726 double result = s_.bp_indel_score_; 727 if (la == lb) 728 result += s_.b_match_score_; 729 else 730 result += s_.b_rep_score_; 731 result += mdown; 732 result += s_.b_indel_score_; 733 result += over; 734 return result; 735 } 736 737 // three/four possibilities for scores, choose via choice function: 738 // insert pairing only, but leave bases 739 // ins p + rep l + mdown + rep r + over 740 // insert pairing and both bases 741 // ins p + ins l + mdown + ins r + over 742 // insert pairing and one base 743 // ins p + ins l + mdown + rep r + over 744 // ins p + rep l + mdown + ins r + over insertPairOnly(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)745 double insertPairOnly(RNA_Alphabet la, RNA_Alphabet lb, double mdown, 746 RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 747 // return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_rep_score_ + over; 748 double result = s_.b_indel_score_; 749 if (la == lb) { 750 result += s_.b_match_score_; 751 } 752 else { 753 result += s_.b_rep_score_; 754 } 755 result += mdown; 756 if (ra == rb) { 757 result += s_.b_match_score_; 758 } 759 else { 760 result += s_.b_rep_score_; 761 } 762 result += over; 763 return result; 764 } insertPairAndBases(RNA_Alphabet lb,double mdown,RNA_Alphabet rb,double over)765 double insertPairAndBases(RNA_Alphabet lb, double mdown, 766 RNA_Alphabet rb, double over) const { 767 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_indel_score_ + over; 768 double result = s_.bp_indel_score_; 769 result += s_.b_indel_score_; 770 result += mdown; 771 result += s_.b_indel_score_; 772 result += over; 773 return result; 774 } insertPairAndLeftBase(RNA_Alphabet lb,double mdown,RNA_Alphabet ra,RNA_Alphabet rb,double over)775 double insertPairAndLeftBase(RNA_Alphabet lb, double mdown, 776 RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 777 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_rep_score_ + over; 778 double result = s_.bp_indel_score_; 779 result += s_.b_indel_score_; 780 result += mdown; 781 if (ra == rb) 782 result += s_.b_match_score_; 783 else 784 result += s_.b_rep_score_; 785 result += over; 786 return result; 787 } insertPairAndRightBase(RNA_Alphabet la,RNA_Alphabet lb,double mdown,RNA_Alphabet rb,double over)788 double insertPairAndRightBase(RNA_Alphabet la, RNA_Alphabet lb, double mdown, 789 RNA_Alphabet rb, double over) const { 790 // return s_.bp_indel_score_ + s_.b_rep_score_ + mdown + s_.b_indel_score_ + over; 791 double result = s_.bp_indel_score_; 792 if (la == lb) 793 result += s_.b_match_score_; 794 else 795 result += s_.b_rep_score_; 796 result += mdown; 797 result += s_.b_indel_score_; 798 result += over; 799 return result; 800 } 801 802 803 }; 804 805 /** RIBOSUM85-60 matrix published in RSEARCH: Finding homologs of single structured RNA sequences 806 * R. Klein and S. Eddy, BMC Bioinformatics 2003 Vol.4 807 */ 808 class RIBOSUM8560 : public RNA_Algebra<double,RNA_Alphabet> { 809 private: 810 double baseSubstMtrx_[4][4]; 811 double basepairSubstMtrx_[4][4][4][4]; 812 Score s_; 813 public: empty()814 double empty() const { 815 return 0; 816 }; replacepair(RNA_Alphabet la,RNA_Alphabet lb,double down,RNA_Alphabet ra,RNA_Alphabet rb,double over)817 double replacepair(RNA_Alphabet la, RNA_Alphabet lb, double down, RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 818 int i,j,k,l; 819 i = alpha2RNA_Alpha(la); 820 j = alpha2RNA_Alpha(ra); 821 k = alpha2RNA_Alpha(lb); 822 l = alpha2RNA_Alpha(rb); 823 824 return basepairSubstMtrx_[i][j][k][l]+down+over; 825 }; 826 replace(RNA_Alphabet a,double down,RNA_Alphabet b,double over)827 double replace(RNA_Alphabet a,double down, RNA_Alphabet b, double over) const { 828 assert(!(a==ALPHA_BASEPAIR && b==ALPHA_BASEPAIR)); 829 830 if (a==ALPHA_BASEPAIR || b==ALPHA_BASEPAIR) 831 return INT_MIN; 832 else { 833 int i,j; 834 i = alpha2RNA_Alpha(a); 835 j = alpha2RNA_Alpha(b); 836 837 return baseSubstMtrx_[i][j]+down+over; 838 } 839 }; 840 del(RNA_Alphabet a,double down,double over)841 double del(RNA_Alphabet a,double down, double over) const { 842 if (a==ALPHA_BASEPAIR) 843 return s_.bp_indel_score_+down+over; 844 else 845 return s_.b_indel_score_+down+over; 846 }; 847 insert(double down,RNA_Alphabet b,double over)848 double insert(double down,RNA_Alphabet b,double over) const { 849 if (b==ALPHA_BASEPAIR) 850 return s_.bp_indel_score_+down+over; 851 else 852 return s_.b_indel_score_+down+over; 853 }; 854 choice(double a,double b)855 double choice(double a, double b) const { 856 return std::max(a,b); 857 }; 858 worst_score()859 double worst_score() const { 860 return INT_MIN; 861 }; 862 RIBOSUM8560(const Score & s)863 RIBOSUM8560(const Score &s) 864 : s_(s) { 865 int i,j,k,l; 866 867 // set substitution matrices 868 869 // base replacement 870 baseSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_A]=222; 871 baseSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_C]=-186; 872 baseSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_G]=-146; 873 baseSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U]=-139; 874 875 baseSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_C]=116; 876 baseSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G]=-248; 877 baseSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_U]=-105; 878 879 baseSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_G]=103; 880 baseSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U]=-174; 881 882 baseSubstMtrx_[ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_U]=165; 883 884 // copy triangle 885 for (i=0; i<=ALPHA_PRO_BASE_U; i++) 886 for (j=0; j<i; j++) 887 baseSubstMtrx_[i][j]=baseSubstMtrx_[j][i]; 888 889 // basepair replacement 890 891 // set default score. This score should never be used since the scores for canonical basepairs are defined later 892 for (i=0; i<=ALPHA_PRO_BASE_U; i++) 893 for (j=0; j<=ALPHA_PRO_BASE_U; j++) 894 for (k=i; k<=ALPHA_PRO_BASE_U; k++) 895 for (l=j; l<=ALPHA_PRO_BASE_U; l++) 896 basepairSubstMtrx_[i][j][k][l]=-1000; 897 898 basepairSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U]=449; 899 basepairSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G]=167; 900 basepairSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C]=270; 901 basepairSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U]=59; 902 basepairSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A]=161; 903 basepairSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G]=-51; 904 905 basepairSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G]=536; 906 basepairSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C]=211; 907 basepairSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U]=-27; 908 basepairSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A]=275; 909 basepairSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G]=132; 910 911 basepairSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C]=562; 912 basepairSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U]=121; 913 basepairSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A]=167; 914 basepairSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G]=-8; 915 916 basepairSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U]=347; 917 basepairSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A]=-57; 918 basepairSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G]=-209; 919 920 basepairSubstMtrx_[ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A]=497; 921 basepairSubstMtrx_[ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G]=114; 922 923 basepairSubstMtrx_[ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G]=336; 924 925 // copy triangle 926 for (i=0; i<=ALPHA_PRO_BASE_U; i++) 927 for (j=0; j<=ALPHA_PRO_BASE_U; j++) 928 for (k=0; k<=ALPHA_PRO_BASE_U; k++) 929 for (l=0; l<=ALPHA_PRO_BASE_U; l++) 930 if (k<i || (k==i && l<j)) 931 basepairSubstMtrx_[i][j][k][l] = basepairSubstMtrx_[k][l][i][j]; 932 }; 933 934 // Pair indel 935 936 // three/four possibilities for scores, choose via choice function: 937 // delete pairing only, but leave bases 938 // del p + rep l + mdown + rep r + over 939 // delete pairing and both bases 940 // del p + del l + mdown + del r + over 941 // delete pairing and one base 942 // del p + del l + mdown + rep r + over 943 // del p + rep l + mdown + del r + over deletePairOnly(const RNA_Alphabet la,const RNA_Alphabet lb,const double mdown,const RNA_Alphabet ra,const RNA_Alphabet rb,const double over)944 double deletePairOnly(const RNA_Alphabet la, const RNA_Alphabet lb, const double mdown, 945 const RNA_Alphabet ra, const RNA_Alphabet rb, const double over) const { 946 947 int i, j, k, l; 948 i = alpha2RNA_Alpha(la); 949 j = alpha2RNA_Alpha(ra); 950 k = alpha2RNA_Alpha(lb); 951 l = alpha2RNA_Alpha(rb); 952 953 return s_.bp_indel_score_ + basepairSubstMtrx_[i][j][k][l] + mdown + over; 954 } deletePairAndBases(const RNA_Alphabet la,const double mdown,const RNA_Alphabet ra,const double over)955 double deletePairAndBases(const RNA_Alphabet la, const double mdown, 956 const RNA_Alphabet ra, const double over) const { 957 958 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_indel_score_ + over; 959 double result = s_.bp_indel_score_; 960 result += s_.b_indel_score_; 961 result += mdown; 962 result += s_.b_indel_score_; 963 result += over; 964 return result; 965 } deletePairAndLeftBase(const RNA_Alphabet la,const double mdown,const RNA_Alphabet ra,const RNA_Alphabet rb,const double over)966 double deletePairAndLeftBase(const RNA_Alphabet la, const double mdown, 967 const RNA_Alphabet ra, const RNA_Alphabet rb, const double over) const { 968 969 int k, l; 970 k = alpha2RNA_Alpha(ra); 971 l = alpha2RNA_Alpha(rb); 972 973 return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + baseSubstMtrx_[k][l] + over; 974 } deletePairAndRightBase(const RNA_Alphabet la,const RNA_Alphabet lb,const double mdown,const RNA_Alphabet ra,const double over)975 double deletePairAndRightBase(const RNA_Alphabet la, const RNA_Alphabet lb, const double mdown, 976 const RNA_Alphabet ra, const double over) const { 977 978 int i, j; 979 i = alpha2RNA_Alpha(la); 980 j = alpha2RNA_Alpha(ra); 981 982 return s_.bp_indel_score_ + baseSubstMtrx_[i][j] + mdown + s_.b_indel_score_ + over; 983 } 984 985 // three/four possibilities for scores, choose via choice function: 986 // insert pairing only, but leave bases 987 // ins p + rep l + mdown + rep r + over 988 // insert pairing and both bases 989 // ins p + ins l + mdown + ins r + over 990 // insert pairing and one base 991 // ins p + ins l + mdown + rep r + over 992 // ins p + rep l + mdown + ins r + over insertPairOnly(const RNA_Alphabet la,const RNA_Alphabet lb,const double mdown,const RNA_Alphabet ra,const RNA_Alphabet rb,const double over)993 double insertPairOnly(const RNA_Alphabet la, const RNA_Alphabet lb, const double mdown, 994 const RNA_Alphabet ra, const RNA_Alphabet rb, const double over) const { 995 996 int i, j, k, l; 997 i = alpha2RNA_Alpha(la); 998 j = alpha2RNA_Alpha(ra); 999 k = alpha2RNA_Alpha(lb); 1000 l = alpha2RNA_Alpha(rb); 1001 1002 return s_.bp_indel_score_ + basepairSubstMtrx_[i][j][k][l] + mdown + over; 1003 } 1004 insertPairAndBases(const RNA_Alphabet lb,const double mdown,const RNA_Alphabet rb,const double over)1005 double insertPairAndBases(const RNA_Alphabet lb, const double mdown, 1006 const RNA_Alphabet rb, const double over) const { 1007 1008 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_indel_score_ + over; 1009 double result = s_.bp_indel_score_; 1010 result += s_.b_indel_score_; 1011 result += mdown; 1012 result += s_.b_indel_score_; 1013 result += over; 1014 return result; 1015 } 1016 insertPairAndLeftBase(const RNA_Alphabet lb,const double mdown,const RNA_Alphabet ra,const RNA_Alphabet rb,const double over)1017 double insertPairAndLeftBase(const RNA_Alphabet lb, const double mdown, 1018 const RNA_Alphabet ra, const RNA_Alphabet rb, const double over) const { 1019 1020 int k, l; 1021 k = alpha2RNA_Alpha(lb); 1022 l = alpha2RNA_Alpha(rb); 1023 1024 return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + baseSubstMtrx_[k][l] + over; 1025 } 1026 insertPairAndRightBase(const RNA_Alphabet la,const RNA_Alphabet lb,const double mdown,const RNA_Alphabet rb,const double over)1027 double insertPairAndRightBase(const RNA_Alphabet la, const RNA_Alphabet lb, const double mdown, 1028 const RNA_Alphabet rb, const double over) const { 1029 1030 int i, j; 1031 i = alpha2RNA_Alpha(lb); 1032 j = alpha2RNA_Alpha(rb); 1033 1034 return s_.bp_indel_score_ + baseSubstMtrx_[i][j] + mdown + s_.b_indel_score_ + over; 1035 } 1036 1037 1038 }; 1039 1040 class AffineRIBOSUM8560 : public RNA_AlgebraAffine<double,RNA_Alphabet> { 1041 private: 1042 Score s_; 1043 double baseSubstMtrx_[4][4]; 1044 double basepairSubstMtrx_[4][4][4][4]; 1045 public: empty()1046 double empty() const { 1047 return 0; 1048 }; replacepair(RNA_Alphabet la,RNA_Alphabet lb,double down,RNA_Alphabet ra,RNA_Alphabet rb,double over)1049 double replacepair(RNA_Alphabet la, RNA_Alphabet lb, double down, RNA_Alphabet ra, RNA_Alphabet rb, double over) const { 1050 int i,j,k,l; 1051 i = alpha2RNA_Alpha(la); 1052 j = alpha2RNA_Alpha(ra); 1053 k = alpha2RNA_Alpha(lb); 1054 l = alpha2RNA_Alpha(rb); 1055 1056 return basepairSubstMtrx_[i][j][k][l]+down+over; 1057 }; 1058 replace(RNA_Alphabet a,double down,RNA_Alphabet b,double over)1059 double replace(RNA_Alphabet a,double down, RNA_Alphabet b, double over) const { 1060 assert(!(a==ALPHA_BASEPAIR && b==ALPHA_BASEPAIR)); 1061 1062 if (a==ALPHA_BASEPAIR || b==ALPHA_BASEPAIR) 1063 return INT_MIN; 1064 else { 1065 int i,j; 1066 i = alpha2RNA_Alpha(a); 1067 j = alpha2RNA_Alpha(b); 1068 1069 return baseSubstMtrx_[i][j]+down+over; 1070 } 1071 }; 1072 del(RNA_Alphabet a,double down,double over)1073 double del(RNA_Alphabet a,double down, double over) const { 1074 if (a==ALPHA_BASEPAIR) 1075 return s_.bp_indel_score_+down+over; 1076 else 1077 return s_.b_indel_score_+down+over; 1078 }; 1079 insert(double down,RNA_Alphabet b,double over)1080 double insert(double down,RNA_Alphabet b,double over) const { 1081 if (b==ALPHA_BASEPAIR) 1082 return s_.bp_indel_score_+down+over; 1083 else 1084 return s_.b_indel_score_+down+over; 1085 }; 1086 choice(double a,double b)1087 double choice(double a, double b) const { 1088 return std::max(a,b); 1089 }; 1090 worst_score()1091 double worst_score() const { 1092 return INT_MIN; 1093 }; 1094 AffineRIBOSUM8560(const Score & s)1095 AffineRIBOSUM8560(const Score &s) 1096 : s_(s) { 1097 int i,j,k,l; 1098 1099 // set substitution matrices 1100 1101 // base replacement 1102 baseSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_A]=222; 1103 baseSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_C]=-186; 1104 baseSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_G]=-146; 1105 baseSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U]=-139; 1106 1107 baseSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_C]=116; 1108 baseSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G]=-248; 1109 baseSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_U]=-105; 1110 1111 baseSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_G]=103; 1112 baseSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U]=-174; 1113 1114 baseSubstMtrx_[ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_U]=165; 1115 1116 // copy triangle 1117 for (i=0; i<=ALPHA_PRO_BASE_U; i++) 1118 for (j=0; j<i; j++) 1119 baseSubstMtrx_[i][j]=baseSubstMtrx_[j][i]; 1120 1121 // basepair replacement 1122 1123 // set default score. This score should never be used since the scores for canonical basepairs are defined later 1124 for (i=0; i<=ALPHA_PRO_BASE_U; i++) 1125 for (j=0; j<=ALPHA_PRO_BASE_U; j++) 1126 for (k=i; k<=ALPHA_PRO_BASE_U; k++) 1127 for (l=j; l<=ALPHA_PRO_BASE_U; l++) 1128 basepairSubstMtrx_[i][j][k][l]=-1000; 1129 1130 basepairSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U]=449; 1131 basepairSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G]=167; 1132 basepairSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C]=270; 1133 basepairSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U]=59; 1134 basepairSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A]=161; 1135 basepairSubstMtrx_[ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G]=-51; 1136 1137 basepairSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G]=536; 1138 basepairSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C]=211; 1139 basepairSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U]=-27; 1140 basepairSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A]=275; 1141 basepairSubstMtrx_[ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G]=132; 1142 1143 basepairSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C]=562; 1144 basepairSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U]=121; 1145 basepairSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A]=167; 1146 basepairSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_C][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G]=-8; 1147 1148 basepairSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U]=347; 1149 basepairSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A]=-57; 1150 basepairSubstMtrx_[ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G]=-209; 1151 1152 basepairSubstMtrx_[ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A]=497; 1153 basepairSubstMtrx_[ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_A][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G]=114; 1154 1155 basepairSubstMtrx_[ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G][ALPHA_PRO_BASE_U][ALPHA_PRO_BASE_G]=336; 1156 1157 // copy triangle 1158 for (i=0; i<=ALPHA_PRO_BASE_U; i++) 1159 for (j=0; j<=ALPHA_PRO_BASE_U; j++) 1160 for (k=0; k<=ALPHA_PRO_BASE_U; k++) 1161 for (l=0; l<=ALPHA_PRO_BASE_U; l++) 1162 if (k<i || (k==i && l<j)) 1163 basepairSubstMtrx_[i][j][k][l]=basepairSubstMtrx_[k][l][i][j]; 1164 }; 1165 delO(RNA_Alphabet a,double down,double over)1166 double delO(RNA_Alphabet a,double down, double over) const { 1167 if (a==ALPHA_BASEPAIR) 1168 return s_.bp_indel_open_score_+down+over; 1169 else 1170 return s_.b_indel_open_score_+down+over; 1171 }; 1172 insertO(double down,RNA_Alphabet b,double over)1173 double insertO(double down,RNA_Alphabet b,double over) const { 1174 if (b==ALPHA_BASEPAIR) 1175 return s_.bp_indel_open_score_+down+over; 1176 else 1177 return s_.b_indel_open_score_+down+over; 1178 }; 1179 1180 // Pair indel 1181 1182 // three/four possibilities for scores, choose via choice function: 1183 // delete pairing only, but leave bases 1184 // del p + rep l + mdown + rep r + over 1185 // delete pairing and both bases 1186 // del p + del l + mdown + del r + over 1187 // delete pairing and one base 1188 // del p + del l + mdown + rep r + over 1189 // del p + rep l + mdown + del r + over deletePairOnly(const RNA_Alphabet la,const RNA_Alphabet lb,const double mdown,const RNA_Alphabet ra,const RNA_Alphabet rb,const double over)1190 double deletePairOnly(const RNA_Alphabet la, const RNA_Alphabet lb, const double mdown, 1191 const RNA_Alphabet ra, const RNA_Alphabet rb, const double over) const { 1192 1193 int i, j, k, l; 1194 i = alpha2RNA_Alpha(la); 1195 j = alpha2RNA_Alpha(ra); 1196 k = alpha2RNA_Alpha(lb); 1197 l = alpha2RNA_Alpha(rb); 1198 1199 return s_.bp_indel_score_ + basepairSubstMtrx_[i][j][k][l] + mdown + over; 1200 } deletePairAndBases(const RNA_Alphabet la,const double mdown,const RNA_Alphabet ra,const double over)1201 double deletePairAndBases(const RNA_Alphabet la, const double mdown, 1202 const RNA_Alphabet ra, const double over) const { 1203 1204 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_indel_score_ + over; 1205 double result = s_.bp_indel_score_; 1206 result += s_.b_indel_score_; 1207 result += mdown; 1208 result += s_.b_indel_score_; 1209 result += over; 1210 return result; 1211 } deletePairAndLeftBase(const RNA_Alphabet la,const double mdown,const RNA_Alphabet ra,const RNA_Alphabet rb,const double over)1212 double deletePairAndLeftBase(const RNA_Alphabet la, const double mdown, 1213 const RNA_Alphabet ra, const RNA_Alphabet rb, const double over) const { 1214 1215 int k, l; 1216 k = alpha2RNA_Alpha(ra); 1217 l = alpha2RNA_Alpha(rb); 1218 1219 return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + baseSubstMtrx_[k][l] + over; 1220 } deletePairAndRightBase(const RNA_Alphabet la,const RNA_Alphabet lb,const double mdown,const RNA_Alphabet ra,const double over)1221 double deletePairAndRightBase(const RNA_Alphabet la, const RNA_Alphabet lb, const double mdown, 1222 const RNA_Alphabet ra, const double over) const { 1223 1224 int i, j; 1225 i = alpha2RNA_Alpha(la); 1226 j = alpha2RNA_Alpha(ra); 1227 1228 return s_.bp_indel_score_ + baseSubstMtrx_[i][j] + mdown + s_.b_indel_score_ + over; 1229 } 1230 1231 // three/four possibilities for scores, choose via choice function: 1232 // insert pairing only, but leave bases 1233 // ins p + rep l + mdown + rep r + over 1234 // insert pairing and both bases 1235 // ins p + ins l + mdown + ins r + over 1236 // insert pairing and one base 1237 // ins p + ins l + mdown + rep r + over 1238 // ins p + rep l + mdown + ins r + over insertPairOnly(const RNA_Alphabet la,const RNA_Alphabet lb,const double mdown,const RNA_Alphabet ra,const RNA_Alphabet rb,const double over)1239 double insertPairOnly(const RNA_Alphabet la, const RNA_Alphabet lb, const double mdown, 1240 const RNA_Alphabet ra, const RNA_Alphabet rb, const double over) const { 1241 1242 int i, j, k, l; 1243 i = alpha2RNA_Alpha(la); 1244 j = alpha2RNA_Alpha(ra); 1245 k = alpha2RNA_Alpha(lb); 1246 l = alpha2RNA_Alpha(rb); 1247 1248 return s_.bp_indel_score_ + basepairSubstMtrx_[i][j][k][l] + mdown + over; 1249 } 1250 insertPairAndBases(const RNA_Alphabet lb,const double mdown,const RNA_Alphabet rb,const double over)1251 double insertPairAndBases(const RNA_Alphabet lb, const double mdown, 1252 const RNA_Alphabet rb, const double over) const { 1253 1254 // return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + s_.b_indel_score_ + over; 1255 double result = s_.bp_indel_score_; 1256 result += s_.b_indel_score_; 1257 result += mdown; 1258 result += s_.b_indel_score_; 1259 result += over; 1260 return result; 1261 } 1262 insertPairAndLeftBase(const RNA_Alphabet lb,const double mdown,const RNA_Alphabet ra,const RNA_Alphabet rb,const double over)1263 double insertPairAndLeftBase(const RNA_Alphabet lb, const double mdown, 1264 const RNA_Alphabet ra, const RNA_Alphabet rb, const double over) const { 1265 1266 int k, l; 1267 k = alpha2RNA_Alpha(lb); 1268 l = alpha2RNA_Alpha(rb); 1269 1270 return s_.bp_indel_score_ + s_.b_indel_score_ + mdown + baseSubstMtrx_[k][l] + over; 1271 } 1272 insertPairAndRightBase(const RNA_Alphabet la,const RNA_Alphabet lb,const double mdown,const RNA_Alphabet rb,const double over)1273 double insertPairAndRightBase(const RNA_Alphabet la, const RNA_Alphabet lb, const double mdown, 1274 const RNA_Alphabet rb, const double over) const { 1275 1276 int i, j; 1277 i = alpha2RNA_Alpha(lb); 1278 j = alpha2RNA_Alpha(rb); 1279 1280 return s_.bp_indel_score_ + baseSubstMtrx_[i][j] + mdown + s_.b_indel_score_ + over; 1281 } 1282 1283 }; 1284 1285 1286 1287 1288 // RNA Profile Algebra Classes 1289 1290 // Similarity algebra for RNA profile forests 1291 class DoubleSimiProfileAlgebra : public Algebra<double,RNA_Alphabet_Profile> { 1292 private: 1293 Score s_; 1294 1295 public: empty()1296 double empty() const { 1297 return 0.0; 1298 }; replace(RNA_Alphabet_Profile a,double down,RNA_Alphabet_Profile b,double over)1299 double replace(RNA_Alphabet_Profile a,double down, RNA_Alphabet_Profile b, double over) const { 1300 if (a.p[ALPHA_PRO_BASEPAIR]>0 && b.p[ALPHA_PRO_BASEPAIR]>0) { 1301 // pair replacement 1302 return a.p[ALPHA_PRO_BASEPAIR]*b.p[ALPHA_PRO_BASEPAIR]*s_.bp_rep_score_ + 1303 down+over; 1304 } else { 1305 if (a.p[ALPHA_PRO_BASE]>0 && b.p[ALPHA_PRO_BASE]>0) { 1306 double s=0; 1307 1308 // base replacement 1309 for (int i=ALPHA_PRO_BASE_A; i<=ALPHA_PRO_BASE_U; i++) 1310 for (int j=ALPHA_PRO_BASE_A; j<=ALPHA_PRO_BASE_U; j++) 1311 s+= i==j ? a.p[i]*b.p[j]*s_.b_match_score_ : a.p[i]*b.p[j]*s_.b_rep_score_; 1312 1313 if (s==0) // no sequence information 1314 s=a.p[ALPHA_PRO_BASE]*b.p[ALPHA_PRO_BASE]*s_.b_rep_score_; 1315 1316 return s+down+over; 1317 } else { 1318 // undefined operation (replace base by basepair ??) 1319 return DBL_NEG/4; 1320 } 1321 } 1322 }; 1323 del(RNA_Alphabet_Profile a,double down,double over)1324 double del(RNA_Alphabet_Profile a,double down, double over) const { 1325 if (a.p[ALPHA_PRO_BASEPAIR]>0) 1326 return a.p[ALPHA_PRO_BASEPAIR]*s_.bp_indel_score_+down+over; 1327 else 1328 return a.p[ALPHA_PRO_BASE]*s_.b_indel_score_+down+over; 1329 }; 1330 insert(double down,RNA_Alphabet_Profile b,double over)1331 double insert(double down,RNA_Alphabet_Profile b,double over) const { 1332 if (b.p[ALPHA_PRO_BASEPAIR]>0) 1333 return b.p[ALPHA_PRO_BASEPAIR]*s_.bp_indel_score_+down+over; 1334 else 1335 return b.p[ALPHA_PRO_BASE]*s_.b_indel_score_+down+over; 1336 }; 1337 choice(double a,double b)1338 double choice(double a, double b) const { 1339 return std::max(a,b); 1340 }; 1341 worst_score()1342 double worst_score() const { 1343 return DBL_NEG; 1344 }; 1345 DoubleSimiProfileAlgebra(const Score & s)1346 DoubleSimiProfileAlgebra(const Score &s) 1347 : s_(s) {}; 1348 }; 1349 1350 class AffineDoubleSimiProfileAlgebra : public AlgebraAffine<double,RNA_Alphabet_Profile> { 1351 private: 1352 Score s_; 1353 public: empty()1354 double empty() const { 1355 return 0.0; 1356 }; replace(RNA_Alphabet_Profile a,double down,RNA_Alphabet_Profile b,double over)1357 double replace(RNA_Alphabet_Profile a,double down, RNA_Alphabet_Profile b, double over) const { 1358 if (a.p[ALPHA_PRO_BASEPAIR]>0 && b.p[ALPHA_PRO_BASEPAIR]>0) { 1359 // pair replacement 1360 return a.p[ALPHA_PRO_BASEPAIR]*b.p[ALPHA_PRO_BASEPAIR]*s_.bp_rep_score_ + 1361 down+over; 1362 } else { 1363 if (a.p[ALPHA_PRO_BASE]>0 && b.p[ALPHA_PRO_BASE]>0) { 1364 double s=0; 1365 1366 // base replacement 1367 for (int i=ALPHA_PRO_BASE_A; i<=ALPHA_PRO_BASE_U; i++) 1368 for (int j=ALPHA_PRO_BASE_A; j<=ALPHA_PRO_BASE_U; j++) 1369 s+= i==j ? a.p[i]*b.p[j]*s_.b_match_score_ : a.p[i]*b.p[j]*s_.b_rep_score_; 1370 1371 if (s==0) // no sequence information 1372 s=a.p[ALPHA_PRO_BASE]*b.p[ALPHA_PRO_BASE]*s_.b_rep_score_; 1373 1374 return s+down+over; 1375 } else { 1376 // undefined operation (replace base by basepair ??) 1377 return DBL_NEG/4; 1378 } 1379 } 1380 }; 1381 del(RNA_Alphabet_Profile a,double down,double over)1382 double del(RNA_Alphabet_Profile a,double down, double over) const { 1383 if (a.p[ALPHA_PRO_BASEPAIR]>0) 1384 return a.p[ALPHA_PRO_BASEPAIR]*s_.bp_indel_score_+down+over; 1385 else 1386 return a.p[ALPHA_PRO_BASE]*s_.b_indel_score_+down+over; 1387 }; 1388 insert(double down,RNA_Alphabet_Profile b,double over)1389 double insert(double down,RNA_Alphabet_Profile b,double over) const { 1390 if (b.p[ALPHA_PRO_BASEPAIR]>0) 1391 return b.p[ALPHA_PRO_BASEPAIR]*s_.bp_indel_score_+down+over; 1392 else 1393 return b.p[ALPHA_PRO_BASE]*s_.b_indel_score_+down+over; 1394 }; 1395 choice(double a,double b)1396 double choice(double a, double b) const { 1397 return std::max(a,b); 1398 }; 1399 worst_score()1400 double worst_score() const { 1401 return DBL_NEG; 1402 }; 1403 1404 delO(RNA_Alphabet_Profile a,double down,double over)1405 double delO(RNA_Alphabet_Profile a,double down, double over) const { 1406 if (a.p[ALPHA_PRO_BASEPAIR]>0) 1407 return a.p[ALPHA_PRO_BASEPAIR]*s_.bp_indel_open_score_+down+over; 1408 else 1409 return a.p[ALPHA_PRO_BASE]*s_.b_indel_open_score_+down+over; 1410 }; 1411 insertO(double down,RNA_Alphabet_Profile b,double over)1412 double insertO(double down,RNA_Alphabet_Profile b,double over) const { 1413 if (b.p[ALPHA_PRO_BASEPAIR]>0) 1414 return b.p[ALPHA_PRO_BASEPAIR]*s_.bp_indel_open_score_+down+over; 1415 else 1416 return b.p[ALPHA_PRO_BASE]*s_.b_indel_open_score_+down+over; 1417 }; 1418 AffineDoubleSimiProfileAlgebra(const Score & s)1419 AffineDoubleSimiProfileAlgebra(const Score &s) 1420 : s_(s) {}; 1421 }; 1422 1423 1424 1425 1426 // Distance algebra for RNA profile forests 1427 class DoubleDistProfileAlgebra : public Algebra<double,RNA_Alphabet_Profile> { 1428 private: 1429 Score s_; 1430 1431 public: empty()1432 double empty() const { 1433 return 0.0; 1434 }; replace(RNA_Alphabet_Profile a,double down,RNA_Alphabet_Profile b,double over)1435 double replace(RNA_Alphabet_Profile a,double down, RNA_Alphabet_Profile b, double over) const { 1436 TRACE(DBG_ALGEBRA,"rep","inside!!!"); 1437 1438 if (a.p[ALPHA_PRO_BASEPAIR]>0 && b.p[ALPHA_PRO_BASEPAIR]>0) { 1439 // pair replacement 1440 return a.p[ALPHA_PRO_BASEPAIR]*b.p[ALPHA_PRO_BASEPAIR]*s_.bp_rep_score_ + 1441 down+over; 1442 } else { 1443 if (a.p[ALPHA_PRO_BASE]>0 && b.p[ALPHA_PRO_BASE]>0) { 1444 double s=0; 1445 1446 // base replacement 1447 for (int i=ALPHA_PRO_BASE_A; i<=ALPHA_PRO_BASE_U; i++) 1448 for (int j=ALPHA_PRO_BASE_A; j<=ALPHA_PRO_BASE_U; j++) 1449 s+= i==j ? a.p[i]*b.p[j]*s_.b_match_score_ : a.p[i]*b.p[j]*s_.b_rep_score_; 1450 1451 if (s==0) // no sequence information 1452 s=a.p[ALPHA_PRO_BASE]*b.p[ALPHA_PRO_BASE]*s_.b_rep_score_; 1453 1454 return s+down+over; 1455 } else { 1456 // undefined operation (replace base by basepair ??) 1457 return DBL_POS/4; 1458 } 1459 } 1460 }; 1461 del(RNA_Alphabet_Profile a,double down,double over)1462 double del(RNA_Alphabet_Profile a,double down, double over) const { 1463 if (a.p[ALPHA_PRO_BASEPAIR]>0) 1464 return a.p[ALPHA_PRO_BASEPAIR]*s_.bp_indel_score_+down+over; 1465 else 1466 return a.p[ALPHA_PRO_BASE]*s_.b_indel_score_+down+over; 1467 }; 1468 insert(double down,RNA_Alphabet_Profile b,double over)1469 double insert(double down,RNA_Alphabet_Profile b,double over) const { 1470 if (b.p[ALPHA_PRO_BASEPAIR]>0) 1471 return b.p[ALPHA_PRO_BASEPAIR]*s_.bp_indel_score_+down+over; 1472 else 1473 return b.p[ALPHA_PRO_BASE]*s_.b_indel_score_+down+over; 1474 }; 1475 choice(double a,double b)1476 double choice(double a, double b) const { 1477 return std::min(a,b); 1478 }; 1479 worst_score()1480 double worst_score() const { 1481 return DBL_POS; 1482 }; 1483 DoubleDistProfileAlgebra(const Score & s)1484 DoubleDistProfileAlgebra(const Score &s) 1485 : s_(s) {}; 1486 }; 1487 1488 1489 // Distance algebra for RNA profile forests 1490 class AffineDoubleDistProfileAlgebra : public AlgebraAffine<double,RNA_Alphabet_Profile> { 1491 private: 1492 Score s_; 1493 public: empty()1494 double empty() const { 1495 return 0.0; 1496 }; replace(RNA_Alphabet_Profile a,double down,RNA_Alphabet_Profile b,double over)1497 double replace(RNA_Alphabet_Profile a,double down, RNA_Alphabet_Profile b, double over) const { 1498 TRACE(DBG_ALGEBRA,"rep","inside!!!"); 1499 1500 if (a.p[ALPHA_PRO_BASEPAIR]>0 && b.p[ALPHA_PRO_BASEPAIR]>0) { 1501 // pair replacement 1502 return a.p[ALPHA_PRO_BASEPAIR]*b.p[ALPHA_PRO_BASEPAIR]*s_.bp_rep_score_ + 1503 down+over; 1504 } else { 1505 if (a.p[ALPHA_PRO_BASE]>0 && b.p[ALPHA_PRO_BASE]>0) { 1506 double s=0; 1507 1508 // base replacement 1509 for (int i=ALPHA_PRO_BASE_A; i<=ALPHA_PRO_BASE_U; i++) 1510 for (int j=ALPHA_PRO_BASE_A; j<=ALPHA_PRO_BASE_U; j++) 1511 s+= i==j ? a.p[i]*b.p[j]*s_.b_match_score_ : a.p[i]*b.p[j]*s_.b_rep_score_; 1512 1513 if (s==0) // no sequence information 1514 s=a.p[ALPHA_PRO_BASE]*b.p[ALPHA_PRO_BASE]*s_.b_rep_score_; 1515 1516 return s+down+over; 1517 } else { 1518 // undefined operation (replace base by basepair ??) 1519 return DBL_POS/4; 1520 } 1521 } 1522 }; 1523 del(RNA_Alphabet_Profile a,double down,double over)1524 double del(RNA_Alphabet_Profile a,double down, double over) const { 1525 if (a.p[ALPHA_PRO_BASEPAIR]>0) 1526 return a.p[ALPHA_PRO_BASEPAIR]*s_.bp_indel_score_+down+over; 1527 else 1528 return a.p[ALPHA_PRO_BASE]*s_.b_indel_score_+down+over; 1529 }; 1530 insert(double down,RNA_Alphabet_Profile b,double over)1531 double insert(double down,RNA_Alphabet_Profile b,double over) const { 1532 if (b.p[ALPHA_PRO_BASEPAIR]>0) 1533 return b.p[ALPHA_PRO_BASEPAIR]*s_.bp_indel_score_+down+over; 1534 else 1535 return b.p[ALPHA_PRO_BASE]*s_.b_indel_score_+down+over; 1536 }; 1537 choice(double a,double b)1538 double choice(double a, double b) const { 1539 return std::min(a,b); 1540 }; 1541 worst_score()1542 double worst_score() const { 1543 return DBL_POS; 1544 }; 1545 1546 delO(RNA_Alphabet_Profile a,double down,double over)1547 double delO(RNA_Alphabet_Profile a,double down, double over) const { 1548 if (a.p[ALPHA_PRO_BASEPAIR]>0) 1549 return a.p[ALPHA_PRO_BASEPAIR]*s_.bp_indel_open_score_+down+over; 1550 else 1551 return a.p[ALPHA_PRO_BASE]*s_.b_indel_open_score_+down+over; 1552 }; 1553 insertO(double down,RNA_Alphabet_Profile b,double over)1554 double insertO(double down,RNA_Alphabet_Profile b,double over) const { 1555 if (b.p[ALPHA_PRO_BASEPAIR]>0) 1556 return b.p[ALPHA_PRO_BASEPAIR]*s_.bp_indel_open_score_+down+over; 1557 else 1558 return b.p[ALPHA_PRO_BASE]*s_.b_indel_open_score_+down+over; 1559 }; 1560 AffineDoubleDistProfileAlgebra(const Score & s)1561 AffineDoubleDistProfileAlgebra(const Score &s) 1562 : s_(s) {}; 1563 }; 1564 1565 // SZAlgebra Classes 1566 1567 class IntSimiSZAlgebra : public SZAlgebra<double,RNA_Alphabet> { 1568 private: 1569 Score s_; 1570 1571 public: empty()1572 double empty() const { 1573 return 0; 1574 }; replace(RNA_Alphabet a,double down,RNA_Alphabet b)1575 double replace(RNA_Alphabet a,double down, RNA_Alphabet b) const { 1576 if (a==ALPHA_BASEPAIR && b == ALPHA_BASEPAIR) 1577 return s_.bp_rep_score_+down; 1578 else { 1579 if (a==ALPHA_BASEPAIR || b==ALPHA_BASEPAIR) 1580 return INT_MIN; 1581 else { 1582 if (a==b) 1583 return s_.b_match_score_+down; 1584 else 1585 return s_.b_rep_score_+down; 1586 } 1587 } 1588 }; 1589 del(RNA_Alphabet a,double down)1590 double del(RNA_Alphabet a,double down) const { 1591 if (a==ALPHA_BASEPAIR) 1592 return s_.bp_indel_score_+down; 1593 else 1594 return s_.b_indel_score_+down; 1595 }; 1596 insert(double down,RNA_Alphabet b)1597 double insert(double down,RNA_Alphabet b) const { 1598 if (b==ALPHA_BASEPAIR) 1599 return s_.bp_indel_score_+down; 1600 else 1601 return s_.b_indel_score_+down; 1602 }; 1603 choice(double a,double b)1604 double choice(double a, double b) const { 1605 return std::max(a,b); 1606 }; 1607 worst_score()1608 double worst_score() const { 1609 return INT_MIN; 1610 }; 1611 IntSimiSZAlgebra(const Score & s)1612 IntSimiSZAlgebra(const Score &s) 1613 : s_(s) {}; 1614 }; 1615 1616 1617 class IntDistSZAlgebra : public SZAlgebra<double,RNA_Alphabet> { 1618 private: 1619 Score s_; 1620 1621 public: empty()1622 double empty() const { 1623 return 0; 1624 }; replace(RNA_Alphabet a,double down,RNA_Alphabet b)1625 double replace(RNA_Alphabet a,double down, RNA_Alphabet b) const { 1626 if (a==ALPHA_BASEPAIR && b == ALPHA_BASEPAIR) 1627 return s_.bp_rep_score_+down; 1628 else { 1629 if (a==ALPHA_BASEPAIR || b==ALPHA_BASEPAIR) 1630 return INT_MAX; 1631 else { 1632 if (a==b) 1633 return s_.b_match_score_+down; 1634 else 1635 return s_.b_rep_score_+down; 1636 } 1637 } 1638 }; 1639 del(RNA_Alphabet a,double down)1640 double del(RNA_Alphabet a,double down) const { 1641 if (a==ALPHA_BASEPAIR) 1642 return s_.bp_indel_score_+down; 1643 else 1644 return s_.b_indel_score_+down; 1645 }; 1646 insert(double down,RNA_Alphabet b)1647 double insert(double down,RNA_Alphabet b) const { 1648 if (b==ALPHA_BASEPAIR) 1649 return s_.bp_indel_score_+down; 1650 else 1651 return s_.b_indel_score_+down; 1652 }; 1653 choice(double a,double b)1654 double choice(double a, double b) const { 1655 return std::min(a,b); 1656 }; 1657 worst_score()1658 double worst_score() const { 1659 return INT_MAX; 1660 }; 1661 IntDistSZAlgebra(const Score & s)1662 IntDistSZAlgebra(const Score &s) 1663 : s_(s) {}; 1664 }; 1665 1666 // General Algebra Classes 1667 1668 // Distance algebra for forests 1669 class DoubleDist_Algebra : public Algebra<double,RNA_Alphabet> { 1670 private: 1671 Score s_; 1672 1673 public: empty()1674 double empty() const { 1675 return 0; 1676 }; 1677 replace(RNA_Alphabet a,double down,RNA_Alphabet b,double over)1678 double replace(RNA_Alphabet a,double down, RNA_Alphabet b, double over) const { 1679 if (a==ALPHA_BASEPAIR && b == ALPHA_BASEPAIR) 1680 return s_.bp_rep_score_+down+over; 1681 else { 1682 if (a==ALPHA_BASEPAIR || b==ALPHA_BASEPAIR) 1683 return INT_MAX; 1684 else { 1685 if (a==b) 1686 return s_.b_match_score_+down+over; 1687 else 1688 return s_.b_rep_score_+down+over; 1689 } 1690 } 1691 }; 1692 del(RNA_Alphabet a,double down,double over)1693 double del(RNA_Alphabet a,double down,double over) const { 1694 if (a==ALPHA_BASEPAIR) 1695 return s_.bp_indel_score_+down+over; 1696 else 1697 return s_.b_indel_score_+down+over; 1698 }; 1699 insert(double down,RNA_Alphabet b,double over)1700 double insert(double down,RNA_Alphabet b,double over) const { 1701 if (b==ALPHA_BASEPAIR) 1702 return s_.bp_indel_score_+down+over; 1703 else 1704 return s_.b_indel_score_+down+over; 1705 }; 1706 choice(double a,double b)1707 double choice(double a, double b) const { 1708 return std::min(a,b); 1709 }; 1710 worst_score()1711 double worst_score() const { 1712 return INT_MAX; 1713 }; 1714 DoubleDist_Algebra(const Score & s)1715 DoubleDist_Algebra(const Score &s) 1716 : s_(s) {}; 1717 }; 1718 1719 1720 // Distance algebra for forests 1721 class AffineDoubleDist_Algebra : public AlgebraAffine<double,RNA_Alphabet> { 1722 private: 1723 Score s_; 1724 1725 public: empty()1726 double empty() const { 1727 return 0; 1728 }; 1729 replace(RNA_Alphabet a,double down,RNA_Alphabet b,double over)1730 double replace(RNA_Alphabet a,double down, RNA_Alphabet b, double over) const { 1731 if (a==ALPHA_BASEPAIR && b == ALPHA_BASEPAIR) 1732 return s_.bp_rep_score_+down+over; 1733 else { 1734 if (a==ALPHA_BASEPAIR || b==ALPHA_BASEPAIR) 1735 return INT_MAX; 1736 else { 1737 if (a==b) 1738 return s_.b_match_score_+down+over; 1739 else 1740 return s_.b_rep_score_+down+over; 1741 } 1742 } 1743 }; 1744 del(RNA_Alphabet a,double down,double over)1745 double del(RNA_Alphabet a,double down,double over) const { 1746 if (a==ALPHA_BASEPAIR) 1747 return s_.bp_indel_score_+down+over; 1748 else 1749 return s_.b_indel_score_+down+over; 1750 }; 1751 insert(double down,RNA_Alphabet b,double over)1752 double insert(double down,RNA_Alphabet b,double over) const { 1753 if (b==ALPHA_BASEPAIR) 1754 return s_.bp_indel_score_+down+over; 1755 else 1756 return s_.b_indel_score_+down+over; 1757 }; 1758 choice(double a,double b)1759 double choice(double a, double b) const { 1760 return std::min(a,b); 1761 }; 1762 worst_score()1763 double worst_score() const { 1764 return INT_MAX; 1765 }; 1766 delO(RNA_Alphabet a,double down,double over)1767 double delO(RNA_Alphabet a,double down,double over) const { 1768 if (a==ALPHA_BASEPAIR) 1769 return s_.bp_indel_open_score_+down+over; 1770 else 1771 return s_.b_indel_open_score_+down+over; 1772 }; 1773 insertO(double down,RNA_Alphabet b,double over)1774 double insertO(double down,RNA_Alphabet b,double over) const { 1775 if (b==ALPHA_BASEPAIR) 1776 return s_.bp_indel_open_score_+down+over; 1777 else 1778 return s_.b_indel_open_score_+down+over; 1779 }; 1780 AffineDoubleDist_Algebra(const Score & s)1781 AffineDoubleDist_Algebra(const Score &s) 1782 : s_(s) {}; 1783 }; 1784 1785 1786 #endif 1787