2 #include "McCaskill.hpp"
3 //#include "energy_param3.hpp"
8 energy_param McCaskill::energyParam;
10 float *McCaskill::exphairpin;
11 float McCaskill::expninio[5][32];
12 float McCaskill::expdangle5[6][4];
13 float McCaskill::expdangle3[6][4];
14 float McCaskill::expinternalLoop[31];
15 float McCaskill::expbulge[31];
16 char McCaskill::exptriLoop[2][6];
17 float McCaskill::exptriLoopEnergy[2];
18 char McCaskill::exptetraLoop[30][7];
19 float McCaskill::exptetraLoopEnergy[30];
20 float McCaskill::expStack[6][6];
21 float McCaskill::expTstackH[6][16];
22 float McCaskill::expTstackI[6][16];
23 float McCaskill::expint11[6][6][16];
24 float McCaskill::expint21[6][6][16][4];
25 float McCaskill::expint22[6][6][16][16];
26 float McCaskill::expMLclosing;
27 float McCaskill::expMLintern[8];
28 float McCaskill::expTermAU;
29 float McCaskill::expMLbase[31];
31 const int McCaskill::TURN = 3;
32 const float McCaskill::GASCONST = 1.98717;
33 const float McCaskill::T = 37 + 273.15;
34 const int McCaskill::MAXLOOP = 30;
35 const int McCaskill::TETRA_ENTH37 = -400;
36 const int McCaskill::NBPAIRS = 7;
37 const int McCaskill::SCALE = 10;
42 calcPartitionFunction()
50 for(int i = 0; i < n_seq; i++) {
51 for(int j = i; j < n_seq; j++) {
52 cout << getProb(i, j) << " ";
63 float *pPointer = p.getPointer(0, 0);
64 float *abPointer = ab.getPointer(0,0);
65 for(int i = 0; i < n_seq; i++) {
66 for(int j = i; j < n_seq; j++) {
67 *pPointer += *abPointer++;
68 *pPointer++ = EXP(*pPointer);
78 for (int i = n_seq - TURN - 2; i >= 0; i--) {
79 float *abPointer = ab.getPointer(i, i + TURN + 1);
80 float *am1Pointer = am1.getPointer(i, i + TURN + 1);
81 float *amPointer = am.getPointer(i, i + TURN + 1);
82 float *q1Pointer = q1.getPointer(i, i + TURN + 1);
83 float *aPointer = a.getPointer(i, i + TURN + 1);
84 int *typePointer = typeMat.getPointer(i, i+TURN+1);
85 for (int j = i + TURN + 1; j < n_seq; j++) {
86 int tmpType = *typePointer++;
87 *abPointer++ = compQb(i, j, tmpType);
88 am1v.ref(i,j) = *am1Pointer++ = compQm1(i,j, tmpType);
89 amv.ref(i,j) = *amPointer++ = compQm(i,j);
90 q1v.ref(i,j) = *q1Pointer++ = compQ1(i,j, tmpType);
91 *aPointer++ = compQ(i,j);
98 compQb(int i, int j, int tmpType)
105 if(Beta::isCanonicalReducedPairCode(type)) {
108 tmpAb = expHairpinEnergy(type, u, i + 1, j - 1);
110 // base pair, bulge, interior loop energy
111 for(int h = i + 1; h <= MIN(i + MAXLOOP + 1, j - TURN - 2); h++) {
113 int max = MAX(h + TURN + 1, j-1-MAXLOOP+u1);
114 float *abPointer = ab.getPointer(h, max - 1);
115 const int *typeMatPointer = typeMat.getPointer(h, max);
117 for(int l = max; l < j; l++) {
118 int type2 = *typeMatPointer++;
120 if(!Beta::isCanonicalReducedPairCode(type2)) continue;
122 assert(h >= 0 && h < n_seq && l >= 0 && l < n_seq);
123 type2 = Beta::flipReducedPairCode(type2);
124 assert(h-i-1 >= 0); assert(j-l-1 >= 0);
125 float loopE = *abPointer;
126 loopE += expLoopEnergy(u1, j-l-1, tmpType, type2, i+1, j-1, h-1, l+1);
127 tmpAb = LOG_ADD(tmpAb, loopE);
132 float tmpQm = IMPOSSIBLE;
133 float *amPointer = am.getPointer(i + 1, j - TURN - 3);
134 float *am1Pointer = am1v.getPointer(j-TURN-2, j-1);
135 for(int h = j - TURN - 2; h >= i + TURN + 3; h--) {
136 assert(h >= 0 && h < n_seq);
137 float tmpScore = *amPointer--;
138 tmpScore += *am1Pointer--;
139 tmpQm = LOG_ADD(tmpQm, tmpScore);
141 tmpQm += expMLclosing + expMLintern[type];
142 tmpQm += endStemScore(i, j);
143 tmpAb = LOG_ADD(tmpAb, tmpQm);
151 //F = a:ML_closing + b:ML_intern*k + c:ML_BASE*u
155 compQm1(int i, int j, int tmpType)
157 float tmpQm1 = IMPOSSIBLE;
160 if (i + TURN + 1 <= l) {
161 int type = typeMat.ref(i,l);
162 if(Beta::isCanonicalReducedPairCode(type)) {
163 float tmpScore = ab.ref(i,l);
164 tmpScore += beginStemScore(i, l);
165 //tmpScore += expMLbase[1]*(j-l) + expMLintern[type];
166 tmpScore += expMLintern[tmpType];
167 tmpQm1 = LOG_ADD(tmpQm1, tmpScore);
171 tmpQm1 = LOG_ADD(tmpQm1, am1.ref(i,j-1));
181 float tmpQm = IMPOSSIBLE;
182 float *amPointer = am.getPointer(i,j-TURN-2);
183 float *am1Pointer = am1v.getPointer(j-TURN-1, j);
184 for(int h = j - TURN - 1; h >= i ; h--) {
186 float tmpAm1 = *am1Pointer--;
189 tmpQm = LOG_ADD(tmpQm, tmpScore);
190 tmpScore = *amPointer--;
192 tmpQm = LOG_ADD(tmpQm, tmpScore);
200 compQ1(int i, int j, int tmpType)
202 float tmpQ1 = IMPOSSIBLE;
204 if(Beta::isCanonicalReducedPairCode(tmpType)) {
205 float tmpScore = ab.ref(i, j);
206 tmpScore += beginStemScore(i, j);
207 tmpQ1 = LOG_ADD(tmpQ1, tmpScore);
209 tmpQ1 = LOG_ADD(tmpQ1, q1.ref(i, j - 1));
219 tmpQ = LOG_ADD(tmpQ, q1.ref(i,j));
221 float *aPointer = a.getPointer(i,j-TURN-2);
222 float *q1Pointer = q1v.getPointer(j-TURN-1, j);
223 for(int h = j - TURN - 1; h >= i + 1; h--) {
224 float tmpScore = *aPointer--;
225 tmpScore += *q1Pointer--;
226 tmpQ = LOG_ADD(tmpQ, tmpScore);
234 beginStemScore(const int i, const int j) const
237 int type = typeMat.ref(i,j);
239 if(0 < i) { temp += expdangle5[type][numSeq[i-1]]; }
240 if(j < n_seq-1) { temp += expdangle3[type][numSeq[j+1]]; }
241 if(type != Beta::REDUCED_CG_CODE && type != Beta::REDUCED_GC_CODE) { temp += expTermAU; }
247 endStemScore(const int i, const int j) const
250 int type = typeMat.ref(i,j);
252 type = Beta::flipReducedPairCode(type);
254 if(i < n_seq-1) { temp += expdangle3[type][numSeq[i+1]]; }
255 if(j > 0) { temp += expdangle5[type][numSeq[j-1]]; }
256 if(type != Beta::REDUCED_CG_CODE && type != Beta::REDUCED_GC_CODE) { temp += expTermAU; }
262 compP(int h, int l, int tmpType)
264 float prob = IMPOSSIBLE;
267 if(Beta::isCanonicalReducedPairCode(type)) {
271 tmp_p -= a.ref(0,n_seq-1);
273 tmp_p += a.ref(0,h-1);
276 tmp_p += a.ref(l+1, n_seq-1);
278 tmp_p += beginStemScore(h, l);
279 prob = LOG_ADD(prob, tmp_p);
281 assert(IMPOSSIBLE <= prob && prob <= 0);
285 int tt = Beta::flipReducedPairCode(tmpType);
286 int max = MAX(0,h-MAXLOOP-1);
287 for(int i = max; i <= h - 1; i++) {
288 float min = MIN(l+MAXLOOP-h+i+2, n_seq-1);
289 int *typeMatPointer = typeMat.getPointer(i,l+1);
290 float *pPointer = p.getPointer(i,l);
291 for(int j = l + 1; j <= min; j++) {
292 int type2 = *typeMatPointer++;
294 if(!Beta::isCanonicalReducedPairCode(type2)) continue;
295 assert(i >= 0 && i < n_seq && j >= 0 && j < n_seq);
297 float tmpScore = *pPointer;
298 tmpScore += expLoopEnergy(h-i-1, j-l-1, type2, tt, i+1, j-1, h-1, l+1);
299 tmp_p = LOG_ADD(tmp_p, tmpScore);
302 prob = LOG_ADD(prob, tmp_p);
303 assert(IMPOSSIBLE <= prob && prob <= 0);
307 float tmp_begin = beginStemScore(h, l);
308 float *q1Pointer = q1v.getPointer(0, l);
309 float *am1Pointer = am1v.getPointer(0, l);
310 float *amPointer = amv.getPointer(1,h-1);
311 for(int i = 0; i <= h-TURN-1; i++) {
312 float tmpq1 = *q1Pointer++;
313 float tmpam = *amPointer++;
314 float tmpScore = *am1Pointer++;
317 tmpScore += tmp_begin;
318 tmpScore += expMLclosing + expMLintern[tt];
319 tmp_p = LOG_ADD(tmp_p, tmpScore);
323 tmpScore += tmp_begin;
324 tmpScore += expMLclosing + expMLintern[tt];
325 tmp_p = LOG_ADD(tmp_p, tmpScore);
328 tmpScore += tmp_begin;
329 tmpScore += expMLclosing + expMLintern[tt];
330 tmp_p = LOG_ADD(tmp_p, tmpScore);
333 assert(IMPOSSIBLE <= tmp_p && tmp_p <= 0);
334 prob = LOG_ADD(prob, tmp_p);
337 for(int i = h-TURN; i <= h-1; i++) {
339 float tmpScore = q1.ref(i,l);
340 tmpScore += tmp_begin;
341 tmpScore += expMLclosing + expMLintern[tt];
342 tmp_p = LOG_ADD(tmp_p, tmpScore);
345 assert(IMPOSSIBLE <= tmp_p && tmp_p <= 0);
346 prob = LOG_ADD(prob, tmp_p);
359 float tmpPm = IMPOSSIBLE;
361 int *typeMatPointer = typeMat.getPointer(i,n_seq-1);
362 float *pPointer = p.getPointer(i,n_seq);
363 float *amPointer = am.getPointer(l+1,n_seq-1);
364 float *abPointer = ab.getPointer(i, n_seq);
365 for(int j = n_seq - 1; j >= l + TURN + 1; j--) {
366 int type = *typeMatPointer--;
370 if(Beta::isCanonicalReducedPairCode(type)) {
371 float tmp = *pPointer;
373 tmp += endStemScore(i, j);
374 tmpPm = LOG_ADD(tmpPm, tmp);
377 tmpPm += expMLintern[1];
384 compPm1(int i, int l)
386 float tmpPm1 = IMPOSSIBLE;
390 int type = typeMat.ref(i,j);
391 if(Beta::isCanonicalReducedPairCode(type)) {
392 float tmp = p.ref(i,j);
393 tmp += endStemScore(i, j);
396 tmpPm1 += expMLintern[1];
398 if(l+1 <= n_seq - 1) {
399 tmpPm1 = LOG_ADD(tmpPm1, am1.ref(i, l+1));
409 for(int h = 0; h <= n_seq - TURN - 2; h++) {
410 float *pPointer = p.getPointer(h, n_seq-1);
411 float *q1Pointer = q1.getPointer(h, n_seq-1);
412 float *am1Pointer = am1.getPointer(h, n_seq-1);
413 int *typePointer = typeMat.getPointer(h, n_seq-1);
414 for(int l = n_seq-1; l >= h + TURN + 1; l--) {
415 int tmpType = *typePointer--;
416 pv.ref(h,l) = *pPointer-- = compP(h,l,tmpType);
417 q1v.ref(h,l) = *q1Pointer-- = compPm(h,l);
418 am1v.ref(h,l) = *am1Pointer-- = compPm1(h,l);
420 assert(p.ref(h,l) <= 0);
430 for(int i = 0; i < n_seq; i++) cout << " " << seq[i];
432 for(int i = 0; i < n_seq; i++) {
436 for(int j = 0; j <= i-1; j++) {
437 if(j != i-1) cout << " ";
440 if(i != 0 && i != n_seq-1) {
444 for(int j = i; j < n_seq; j++) {
445 if(p.ref(i,j) > 0.01) {
447 int type = Beta::getReducedPairCode(numSeq[i], numSeq[j]);
449 if(!Beta::isCanonicalReducedPairCode(type)) {
450 cout << "\n" << seq[i] << " " << seq[j] << " " << exp(p.ref(i,j)) << endl;
474 cout << seq[m++] << endl;
476 if(i == n_seq - 1) cout << endl;
478 for(int i = 0; i < n_seq; i++) cout << " " << seq[i];
487 float RT_KCAL_MOL = McCaskill::T*McCaskill::GASCONST;
490 for(int i = 0; i < len; i++) {
491 GT = energyParam.getHairpin(i);
492 exphairpin[i] = -GT*10/RT_KCAL_MOL;
495 for (int i = 0; i < len; i++) {
496 GT = energyParam.getBulge(i);
497 expbulge[i] = -GT*10/RT_KCAL_MOL;
498 GT = energyParam.getInternalLoop(i);
499 expinternalLoop[i] = -GT*10/RT_KCAL_MOL;
501 expinternalLoop[2] = -80*10/RT_KCAL_MOL; /* special case of size 2 interior loops (single mismatch) */
503 // float lxc = energy_param3::lxc37;
504 for (int i = 31; i < n_seq; i++) {
505 GT = energyParam.getHairpin(30) + (107.856*LOG((float)i/30));
506 exphairpin[i] = -GT*10/RT_KCAL_MOL;
509 for(int i = 0; i < 5; i++) {
510 GT = energyParam.getNinio(i);
511 for(int j = 0; j <= MAXLOOP; j++) {
512 expninio[i][j] = -MIN(energyParam.getMaxNinio(), j*GT)*10/RT_KCAL_MOL;
516 for(int i = 0; i < 30; i++) {
517 GT = energyParam.getTetraLoopEnergy(i);
518 exptetraLoopEnergy[i] = -GT*10/RT_KCAL_MOL;
521 /*no parameters for Triloop*/
522 for(int i = 0; i < 2; i++) {
524 exptriLoopEnergy[i] = -GT*10/RT_KCAL_MOL;
527 GT = energyParam.getMLclosing();
528 expMLclosing = -GT*10/RT_KCAL_MOL;
530 for(int i = 0; i <= NBPAIRS; i++) {
531 GT = energyParam.getMLintern();
532 expMLintern[i] = -GT*10/RT_KCAL_MOL;
535 expTermAU = -energyParam.getTerminalAU()*10/RT_KCAL_MOL;
537 GT = energyParam.getMLBase();
538 for(int i = 0; i < len; i++) {
539 expMLbase[i] = -GT*10*(float)i/RT_KCAL_MOL;
543 if danlges = 0 just set their energy to 0,
544 don't let dangle energyies become > 0 (at large temps)
545 but make sure go smoothly to 0
547 for(int i = 0; i < 6; i++) {
548 for(int j =0; j < 4; j++) {
549 GT = energyParam.getDangle5(i,j);
550 expdangle5[i][j] = -GT*10/RT_KCAL_MOL;
551 GT = energyParam.getDangle3(i,j);
552 expdangle3[i][j] = -GT*10/RT_KCAL_MOL;
556 /* stacking energies */
557 for(int i = 0; i < 6; i++) {
558 for(int j = 0; j < 6; j++) {
559 GT = energyParam.getStack(i,j);
560 expStack[i][j] = -GT*10/RT_KCAL_MOL;
564 /* mismatch energies */
565 for (int i = 0; i < 6; i++) {
566 for (int j = 0; j < 16; j++) {
567 GT = energyParam.getTstackI(i, j);
568 // cout << i << " " << " " << j << " " << GT << endl;
569 expTstackI[i][j] = -GT*10/RT_KCAL_MOL;
570 GT = energyParam.getTstackH(i, j);
571 expTstackH[i][j] = -GT*10/RT_KCAL_MOL;
575 /* interior loops of length 2*/
576 for(int i = 0; i < 6; i++) {
577 for(int j = 0; j < 6; j++) {
578 for(int k = 0; k < 16; k++) {
579 GT = energyParam.getInt11(i, j, k);
580 expint11[i][j][k] = -GT*10/RT_KCAL_MOL;
585 /* interior 2*1 loops */
586 for(int i = 0; i < 6; i++) {
587 for(int j =0; j < 6; j++) {
588 for(int k = 0; k < 16; k++) {
589 for(int l = 0; l < 4; l++) {
590 GT = energyParam.getInt21(i,j,k,l);
591 expint21[i][j][k][l] = -GT*10/RT_KCAL_MOL;
597 /* interior 2*2 loops */
598 for (int i = 0; i < 6; i++) {
599 for(int j = 0; j < 6; j++) {
600 for(int k = 0; k < 16; k++) {
601 for(int l = 0; l < 16; l++) {
602 GT = energyParam.getInt22(i,j,k,l);
603 expint22[i][j][k][l] = -GT*10/RT_KCAL_MOL;
613 expHairpinEnergy(const int type, const int l, const int i, const int j)
624 for(int iter = i - 1; iter < i + 5; iter++) {
625 temp_seq[iter - (i-1)] = seq[iter];
629 for(k = 0; k < 30; k++) {
630 if(strcmp(temp_seq, energyParam.getTetraLoop(k)) == 0) break;
633 q += exptetraLoopEnergy[k];
641 for(int iter = i - 1; iter < i + 4; iter++) {
642 temp_seq[iter - (i-1)] = seq[iter];
645 for(k = 0; k < 2; k++) {
646 if(strcmp(temp_seq, energyParam.getTriLoop(k)) == 0) break;
649 q *= exptriLoopEnergy[k];
653 if(type != Beta::REDUCED_CG_CODE && type != Beta::REDUCED_GC_CODE) q += expTermAU;
656 int type2 = Beta::getPairCode(numSeq[i], numSeq[j]);
657 q += expTstackH[type][type2];
666 expLoopEnergy(int u1, int u2, int type, int type2,
667 int si1, int sj1, int sp1, int sq1)
671 if((u1 == 0) && (u2 == 0)) { z = expStack[type][type2]; }
673 if((u1 == 0) || (u2 == 0)) {
675 if(u1 == 0) { u = u2; }
678 if(u1 + u2 == 1) z += expStack[type][type2];
680 if (type != Beta::REDUCED_CG_CODE && type != Beta::REDUCED_GC_CODE) z += expTermAU;
681 if (type2 != Beta::REDUCED_CG_CODE && type2 != Beta::REDUCED_GC_CODE) z += expTermAU;
686 z = expint11[type][type2][Beta::getPairCode(numSeq[si1], numSeq[sj1])];
688 else if((u1 == 1) && (u2 == 2)) {
689 z = expint21[type][type2][Beta::getPairCode(numSeq[si1], numSeq[sj1])][numSeq[sq1]];
691 else if((u1 == 2) && (u2 == 1)) {
692 z = expint21[type2][type][Beta::getPairCode(numSeq[sq1], numSeq[sp1])][numSeq[si1]];
694 else if((u1 == 2) && (u2 == 2)) {
695 z = expint22[type][type2][Beta::getPairCode(numSeq[si1], numSeq[sj1])][Beta::getPairCode(numSeq[sp1], numSeq[sq1])];
698 z = expinternalLoop[u1 + u2] +
699 expTstackI[type][Beta::getPairCode(numSeq[si1], numSeq[sj1])]
700 + expTstackI[type2][Beta::getPairCode(numSeq[sq1], numSeq[sp1])];
701 z += expninio[2][abs(u1-u2)];
713 cout << "exphairpin:" << endl;
714 for(int i = 0; i < 31; i++) {
715 cout << exphairpin[i] << endl;
718 cout << "expninio[5][32]:" << endl;
719 for(int i = 0; i < 5; i++) {
720 for(int j = 0; j < 32; j++) {
721 cout << expninio[i][j] << " ";
726 cout << "expdangle5[6][4]:" << endl;
727 for(int i = 0; i < 6; i++) {
728 for(int j = 0; j < 4; j++) {
729 cout << expdangle5[i][j] << " ";
734 cout << "expdangle3[6][4]:" << endl;
735 for(int i = 0; i < 6; i++) {
736 for(int j = 0; j < 4; j++) {
737 cout << expdangle3[i][j] << " ";
742 cout << "expinternalLoop[31]:" << endl;
743 for(int i = 0; i < 31; i++) {
744 cout << i << ":" << expinternalLoop[i] << endl;
746 cout << "expbulge[31]:" << endl;
747 for(int i = 0; i < 31; i++) {
748 cout << i << ":" << expbulge[i] << endl;
751 cout << "exptriLoopEnergy[2]:" << endl;
752 for(int i = 0; i < 2; i++) {
753 cout << i << ":" << exptriLoopEnergy[i] << endl;
756 cout << "exptetraLoopEnergy[15]" << endl;
757 for(int i = 0; i < 15; i++) {
758 cout << i << ":" << exptetraLoopEnergy[i] << endl;
761 cout << "expStack[6][6]:" << endl;
762 for(int i = 0; i < 6; i++) {
763 for(int j = 0; j < 6; j++) {
764 cout << expStack[i][j] << " ";
769 cout << "expTstackH[6][16]:" << endl;
770 for(int i = 0; i < 6; i++) {
771 for(int j = 0; j < 16; j++) {
772 cout << expTstackH[i][j] << " ";
777 cout << "expTstackI[6][16]:" << endl;
778 for(int i = 0; i < 6; i++) {
779 for(int j = 0; j < 16; j++) {
780 cout << expTstackI[i][j] << " ";
785 cout << "expMLclosing=" << expMLclosing << endl;
786 cout << "expMLintern:" << endl;
787 for(int i = 0; i < 8; i++) {
788 cout << expMLintern[i] << " ";
792 cout << "expMLbase[31]:";
793 for(int i = 0; i < 31; i++) {
794 cout << i << ":" << expMLbase[i] << endl;
797 cout << "expint11[6][6][16]:";
798 for(int i = 0; i < 6; i++) {
799 for(int j = 0; j < 6; j++) {
800 for(int k = 0; k < 16; k++) {
801 cout << expint11[i][j][k] << " ";
808 cout << "expint21[6][6][16][4]:" << endl;
809 for(int i = 0; i < 6; i++) {
810 for(int j = 0; j < 6; j++) {
811 for(int k = 0; k < 16; k++) {
812 for(int l = 0; l < 4; l++) {
813 cout << expint21[i][j][k][l] << " ";
823 cout << "expint22[6][6][16][16]:" << endl;
824 for(int i = 0; i < 6; i++) {
825 for(int j = 0; j < 6; j++) {
826 for(int k = 0; k < 16; k++) {
827 for(int l = 0; l < 16; l++) {
828 cout << expint22[i][j][k][l] << " ";