| 1 | |
|---|
| 2 | #ifndef MCCAKILL_H |
|---|
| 3 | #define MCCAKILL_H |
|---|
| 4 | |
|---|
| 5 | //#define NDEBUG |
|---|
| 6 | |
|---|
| 7 | #include <string> |
|---|
| 8 | #include <iostream> |
|---|
| 9 | #include "nrutil.h" |
|---|
| 10 | #include <cassert> |
|---|
| 11 | #include "Util.hpp" |
|---|
| 12 | #include "Beta.hpp" |
|---|
| 13 | |
|---|
| 14 | #include "energy_param.hpp" |
|---|
| 15 | //#include "energy_param3.hpp" |
|---|
| 16 | |
|---|
| 17 | #include "ScoreType.hpp" |
|---|
| 18 | |
|---|
| 19 | |
|---|
| 20 | using namespace std; |
|---|
| 21 | using namespace ProbCons; |
|---|
| 22 | |
|---|
| 23 | namespace MXSCARNA { |
|---|
| 24 | class McCaskill { |
|---|
| 25 | char *seq; |
|---|
| 26 | int *numSeq; |
|---|
| 27 | int n_seq; |
|---|
| 28 | static energy_param energyParam; |
|---|
| 29 | |
|---|
| 30 | Trimat<float> a, q1, ab, am, am1, p; |
|---|
| 31 | TriVertMat<float> q1v, abv, amv, am1v, pv; |
|---|
| 32 | Trimat<int> typeMat; |
|---|
| 33 | |
|---|
| 34 | void printExpEnergy(); |
|---|
| 35 | void initParameter(); |
|---|
| 36 | void Inside(); |
|---|
| 37 | void Outside(); |
|---|
| 38 | void convertProbability(); |
|---|
| 39 | |
|---|
| 40 | inline float expHairpinEnergy(const int type, const int l, const int i, const int j); |
|---|
| 41 | inline float expLoopEnergy(int u1, int u2, int type, int type2, |
|---|
| 42 | int si1, int sj1, int sp1, int sq1); |
|---|
| 43 | |
|---|
| 44 | inline float compQb(int i, int j, int tmpType); |
|---|
| 45 | inline float compQm1(int i, int j, int tmpType); |
|---|
| 46 | inline float compQm(int i, int j); |
|---|
| 47 | inline float compQ1(int i, int j, int tmpType); |
|---|
| 48 | inline float compQ(int i, int j); |
|---|
| 49 | inline float compP(int h, int l, int tmpType); |
|---|
| 50 | inline float compPm(int i, int l); |
|---|
| 51 | inline float compPm1(int i, int l); |
|---|
| 52 | inline float beginStemScore(const int i, const int j) const; |
|---|
| 53 | inline float endStemScore(const int i, const int j) const; |
|---|
| 54 | |
|---|
| 55 | static const float GASCONST; |
|---|
| 56 | static const float T; |
|---|
| 57 | static const int MAXLOOP; |
|---|
| 58 | static const int TETRA_ENTH37; |
|---|
| 59 | static const int NBPAIRS; |
|---|
| 60 | static const int SCALE; |
|---|
| 61 | static const int TURN; |
|---|
| 62 | static float *exphairpin; |
|---|
| 63 | static float expninio[5][32]; |
|---|
| 64 | static float expdangle5[6][4]; |
|---|
| 65 | static float expdangle3[6][4]; |
|---|
| 66 | static float expinternalLoop[31]; |
|---|
| 67 | static float expbulge[31]; |
|---|
| 68 | static char exptriLoop[2][6]; |
|---|
| 69 | static float exptriLoopEnergy[2]; |
|---|
| 70 | static char exptetraLoop[30][7]; |
|---|
| 71 | static float exptetraLoopEnergy[30]; |
|---|
| 72 | static float expStack[6][6]; |
|---|
| 73 | static float expTstackH[6][16]; |
|---|
| 74 | static float expTstackI[6][16]; |
|---|
| 75 | static float expint11[6][6][16]; |
|---|
| 76 | static float expint21[6][6][16][4]; |
|---|
| 77 | static float expint22[6][6][16][16]; |
|---|
| 78 | static float expMLclosing; |
|---|
| 79 | static float expMLintern[8]; |
|---|
| 80 | static float expTermAU; |
|---|
| 81 | static float expMLbase[31]; |
|---|
| 82 | |
|---|
| 83 | public: |
|---|
| 84 | |
|---|
| 85 | McCaskill(int n, const char *mySeq) { |
|---|
| 86 | seq = new char[n + 1]; |
|---|
| 87 | numSeq = new int[n + 1]; |
|---|
| 88 | n_seq = 0; |
|---|
| 89 | |
|---|
| 90 | |
|---|
| 91 | for(int i = 0; i < n; i++) { |
|---|
| 92 | if (mySeq[i] == 'a' || mySeq[i] == 'A') { seq[n_seq] = 'A'; numSeq[n_seq] = Beta::A_CODE; n_seq++; } |
|---|
| 93 | else if(mySeq[i] == 't' || mySeq[i] == 'T' || |
|---|
| 94 | mySeq[i] == 'u' || mySeq[i] == 'U') { seq[n_seq] = 'U'; numSeq[n_seq] = Beta::U_CODE; n_seq++; } |
|---|
| 95 | else if(mySeq[i] == 'g' || mySeq[i] == 'G') { seq[n_seq] = 'G'; numSeq[n_seq] = Beta::G_CODE; n_seq++; } |
|---|
| 96 | else if(mySeq[i] == 'c' || mySeq[i] == 'C') { seq[n_seq] = 'C'; numSeq[n_seq] = Beta::C_CODE; n_seq++; } |
|---|
| 97 | else if(mySeq[i] == 'n' || mySeq[i] == 'N') { seq[n_seq] = 'N'; numSeq[n_seq] = Beta::N_CODE; n_seq++; } |
|---|
| 98 | else if(mySeq[i] == '.' || mySeq[i] == '-') { seq[n_seq] = '-'; numSeq[n_seq] = Beta::GAP_CODE; n_seq++; } |
|---|
| 99 | else { seq[n_seq] = mySeq[i]; numSeq[n_seq] = Beta::INVALID_CODE; n_seq++; } |
|---|
| 100 | } |
|---|
| 101 | seq[n_seq] = '\0'; |
|---|
| 102 | a.Allocator(n_seq); |
|---|
| 103 | q1.Allocator(n_seq); |
|---|
| 104 | ab.Allocator(n_seq); |
|---|
| 105 | am.Allocator(n_seq); |
|---|
| 106 | am1.Allocator(n_seq); |
|---|
| 107 | p.Allocator(n_seq); |
|---|
| 108 | q1v.Allocator(n_seq); |
|---|
| 109 | abv.Allocator(n_seq); |
|---|
| 110 | amv.Allocator(n_seq); |
|---|
| 111 | am1v.Allocator(n_seq); |
|---|
| 112 | pv.Allocator(n_seq); |
|---|
| 113 | typeMat.Allocator(n_seq); |
|---|
| 114 | |
|---|
| 115 | if(n_seq > 31) { |
|---|
| 116 | exphairpin = new float[n_seq + 1]; |
|---|
| 117 | } |
|---|
| 118 | else { |
|---|
| 119 | exphairpin = new float[31]; |
|---|
| 120 | } |
|---|
| 121 | |
|---|
| 122 | for(int i = 0; i < n_seq; i++) { |
|---|
| 123 | for(int j = i; j < n_seq; j++) { |
|---|
| 124 | a.ref(i,j) = q1.ref(i,j) = IMPOSSIBLE; |
|---|
| 125 | q1v.ref(i,j) = IMPOSSIBLE; |
|---|
| 126 | } |
|---|
| 127 | } |
|---|
| 128 | |
|---|
| 129 | for(int i = 0; i < n_seq; i++) { |
|---|
| 130 | a.ref(i,i) = 0.0; |
|---|
| 131 | q1.ref(i,i) = IMPOSSIBLE; |
|---|
| 132 | q1v.ref(i,i) = IMPOSSIBLE; |
|---|
| 133 | } |
|---|
| 134 | |
|---|
| 135 | for(int i = 0; i < n_seq-1; i++) { |
|---|
| 136 | a.ref(i,i+1) = 0.0; |
|---|
| 137 | q1.ref(i,i+1) = IMPOSSIBLE; |
|---|
| 138 | q1v.ref(i,i+1) = IMPOSSIBLE; |
|---|
| 139 | } |
|---|
| 140 | |
|---|
| 141 | for(int i = 0; i < n_seq-2; i++) { |
|---|
| 142 | a.ref(i,i+2) = 0.0; |
|---|
| 143 | q1.ref(i,i+2) = IMPOSSIBLE; |
|---|
| 144 | q1v.ref(i,i+2) = IMPOSSIBLE; |
|---|
| 145 | } |
|---|
| 146 | |
|---|
| 147 | for(int i = 0; i < n_seq-3; i++) { |
|---|
| 148 | a.ref(i,i+3) = 0.0; |
|---|
| 149 | q1.ref(i,i+3) = IMPOSSIBLE; |
|---|
| 150 | q1v.ref(i,i+3) = IMPOSSIBLE; |
|---|
| 151 | |
|---|
| 152 | } |
|---|
| 153 | |
|---|
| 154 | for(int i = 0; i < n_seq; i++) { |
|---|
| 155 | for(int j = i; j < n_seq; j++) { |
|---|
| 156 | ab.ref(i,j) = am.ref(i,j) = am1.ref(i,j) = p.ref(i,j) = IMPOSSIBLE; |
|---|
| 157 | abv.ref(i,j) = amv.ref(i,j) = am1v.ref(i,j) = pv.ref(i,j) = IMPOSSIBLE; |
|---|
| 158 | } |
|---|
| 159 | } |
|---|
| 160 | |
|---|
| 161 | /* the type of base pair */ |
|---|
| 162 | /* C <-> G : type 1 */ |
|---|
| 163 | /* G <-> C : type 2 */ |
|---|
| 164 | /* G <-> U : type 3 */ |
|---|
| 165 | /* U <-> G : type 5 */ |
|---|
| 166 | /* A <-> U : type 0 */ |
|---|
| 167 | /* U <-> A : type 4 */ |
|---|
| 168 | /* ? <-> ? : type 6 */ |
|---|
| 169 | for(int i = 0; i < n_seq; i++) { |
|---|
| 170 | for(int j = i; j < n_seq; j++) { |
|---|
| 171 | typeMat.ref(i,j) = Beta::getReducedPairCode(numSeq[i], numSeq[j]); |
|---|
| 172 | } |
|---|
| 173 | } |
|---|
| 174 | |
|---|
| 175 | } |
|---|
| 176 | |
|---|
| 177 | /*------------------------------------------------------------------------*/ |
|---|
| 178 | /* dangling ends should never be destabilizing, i.e. expdangle>=1 */ |
|---|
| 179 | /* specific heat needs smooth function (2nd derivative) */ |
|---|
| 180 | /* we use a*(sin(x+b)+1)^2, with a=2/(3*sqrt(3)), b=Pi/6-sqrt(3)/2, */ |
|---|
| 181 | /* in the interval b<x<sqrt(3)/2 */ |
|---|
| 182 | float SMOOTH(float X) { |
|---|
| 183 | return ((X)/SCALE<-1.2283697)?0:(((X)/SCALE>0.8660254)?(X): |
|---|
| 184 | SCALE*0.38490018*(sin((X)/SCALE-0.34242663)+1)*(sin((X)/SCALE-0.34242663)+1)); |
|---|
| 185 | } |
|---|
| 186 | |
|---|
| 187 | ~McCaskill() { |
|---|
| 188 | delete[] seq; |
|---|
| 189 | delete[] numSeq; |
|---|
| 190 | delete[] exphairpin; |
|---|
| 191 | } |
|---|
| 192 | |
|---|
| 193 | void calcPartitionFunction(); |
|---|
| 194 | void printProbMat(); |
|---|
| 195 | |
|---|
| 196 | inline float getProb(const int i, const int j) const { |
|---|
| 197 | // 0 origin : 0..(n-1) |
|---|
| 198 | return p.ref(i, j); |
|---|
| 199 | } |
|---|
| 200 | }; |
|---|
| 201 | } |
|---|
| 202 | #endif // MCCASKILL_H |
|---|
