source: tags/arb-6.0/GDE/MAFFT/mafft-7.055-with-extensions/extensions/mxscarna_src/Globaldp.cpp

Last change on this file was 10371, checked in by aboeckma, 11 years ago

updated mafft version. Added extensions (no svn ignore, yet)
File size: 15.2 KB
Line 
1#include "Globaldp.hpp"
2
3namespace MXSCARNA {
4
5double Globaldp::ribosum_matrix[16][16]
6= { { -2.49,  -7.04,  -8.24,  -4.32,  -8.84,  -14.37,  -4.68, -12.64, -6.86, -5.03, -8.39,  -5.84, -4.01, -11.32, -6.16, -9.05 },
7    { -7.04,  -2.11,  -8.89,  -2.04,  -9.37,  -9.08,   -5.86, -10.45, -9.73, -3.81, -11.05, -4.72, -5.33, -8.67,  -6.93, -7.83 },
8    { -8.24,  -8.89,  -0.80,  -5.13,  -10.41, -14.53,  -4.57, -10.14, -8.61, -5.77, -5.38,  -6.60, -5.43, -8.87,  -5.94, -11.07 },
9    { -4.32,  -2.04,  -5.13,   4.49,  -5.56,  -6.71,    1.67, -5.17,  -5.33,  2.70, -5.61,   0.59,  1.61, -4.81,  -0.51, -2.98 },
10    { -8.84,  -9.37,  -10.41, -5.56,  -5.13,  -10.45,  -3.57, -8.49,  -7.98, -5.95, -11.36, -7.93, -2.42, -7.08,  -5.63, -8.39 },
11    { -14.37, -9.08,  -14.53, -6.71,  -10.45, -3.59,   -5.71, -5.77, -12.43, -3.70, -12.58, -7.88, -6.88, -7.40,  -8.41, -5.41 },
12    { -4.68,  -5.86,  -4.57,   1.67,  -3.57,  -5.71,    5.36, -4.96,  -6.00,  2.11, -4.66,  -0.27,  2.75, -4.91,   1.32, -3.67 },
13    { -12.64, -10.45, -10.14, -5.17,  -8.49,  -5.77,   -4.96, -2.28,  -7.71, -5.84, -13.69, -5.61, -4.72, -3.83,  -7.36, -5.21 },
14    { -6.86,  -9.73,  -8.61,  -5.33,  -7.98,  -12.43,  -6.00, -7.71,  -1.05, -4.88, -8.67,  -6.10, -5.85, -6.63,  -7.55, -11.54 },
15    { -5.03,  -3.81,  -5.77,   2.70,  -5.95,  -3.70,    2.11, -5.84,  -4.88,  5.62, -4.13,   1.21,  1.60, -4.49,  -0.08, -3.90 },
16    { -8.39,  -11.05, -5.38,  -5.61,  -11.36, -12.58,  -4.66, -13.69, -8.67, -4.13, -1.98,  -5.77, -5.75, -12.01, -4.27, -10.79 },
17    { -5.84,  -4.72,  -6.60,   0.59,  -7.93,  -7.88,   -0.27, -5.61,  -6.10,  1.21, -5.77,   3.47, -0.57, -5.30,  -2.09, -4.45 },
18    { -4.01,  -5.33,  -5.43,   1.61,  -2.42,  -6.88,    2.75, -4.72,  -5.85,  1.60, -5.75,  -0.57,  4.97, -2.98,   1.14, -3.39 },
19    { -11.32, -8.67,  -8.87,  -4.81,  -7.08,  -7.40,   -4.91, -3.83,  -6.63, -4.49, -12.01, -5.30, -2.98, -3.21,  -4.76, -5.97 },
20    { -6.16,  -6.93,  -5.94,  -0.51,  -5.63,  -8.41,    1.32, -7.36,  -7.55, -0.08, -4.27,  -2.09,  1.14, -4.76,   3.36, -4.28 },
21    { -9.05,  -7.83,  -11.07, -2.98,  -8.39,  -5.41,   -3.67, -5.21, -11.54, -3.90, -10.79, -4.45, -3.39, -5.97,  -4.28, -0.02 }
22};
23
24
25Trimat<float>* 
26Globaldp::
27makeProfileBPPMatrix(const MultiSequence *Sequences)
28{
29    int length = Sequences->GetSequence(0)->GetLength();
30    float thr  = THR;
31    Trimat<float> *consBppMat = new Trimat<float>(length + 1);
32    fill(consBppMat->begin(), consBppMat->end(), 0);
33
34    int number = Sequences->GetNumSequences();
35    for(int seqNum = 0; seqNum < number; seqNum++) {
36        SafeVector<int> *tmpMap = Sequences->GetSequence(seqNum)->GetMappingNumber();
37        int label = Sequences->GetSequence(seqNum)->GetLabel();
38        BPPMatrix *tmpBppMatrix = BPPMatrices[label];
39           
40        for(int i = 1; i <= length ; i++) {
41            int originI = tmpMap->at(i);
42            for(int j = i; j <= length; j++) {
43                int originJ = tmpMap->at(j);
44                if(originI != 0 && originJ != 0) {
45                    float tmpProb = tmpBppMatrix->GetProb(originI, originJ);
46
47                    if(tmpProb >= thr) {
48                        consBppMat->ref(i, j) += tmpProb;
49                    }
50                }
51            }
52        }
53    }
54
55        /* compute the mean of base pairing probability  */
56    for(int i = 1; i <= length; i++) {
57        for(int j = i; j <= length; j++) {
58            consBppMat->ref(i,j) = consBppMat->ref(i,j)/(float)number;
59        }
60    }
61
62    return consBppMat;   
63}
64
65float 
66Globaldp::
67incrementalScorePSCPair(const StemCandidate &sc1, const StemCandidate &sc2)
68{
69    int wordLength = WORDLENGTH;
70
71    int pos1, rvpos1;
72    if (sc1.GetDistance() > 0) {
73        pos1   = sc1.GetPosition();
74        rvpos1 = sc1.GetRvposition();
75    }
76    else {
77        pos1   = sc1.GetRvposition();
78        rvpos1 = sc1.GetPosition();
79    }
80
81    int pos2, rvpos2;
82    if (sc2.GetDistance() > 0) {
83        pos2   = sc2.GetPosition();
84        rvpos2 = sc2.GetRvposition();
85    }
86    else {
87        pos2   = sc2.GetRvposition();
88        rvpos2 = sc2.GetPosition();
89    }
90/*
91    cout << "1:" << i << " " << j << " " << sc1.GetDistance() << " " << sc2.GetDistance() << endl;
92    cout << sc1.GetPosition() << " " << sc1.GetRvposition() << " " << sc2.GetPosition() << " " << sc2.GetRvposition() << endl;
93*/ 
94    float score = 
95          posterior[sc1.GetPosition() + wordLength - 1][sc2.GetPosition() + wordLength - 1]
96//      * sc1.GetBaseScore(wordLength - 1)
97        * profileBPPMat1->ref(pos1 + wordLength - 1, rvpos1)
98        * posterior[sc1.GetRvposition()][sc2.GetRvposition()]
99//      * sc2.GetBaseScore(wordLength - 1);
100        * profileBPPMat2->ref(pos2 + wordLength - 1, rvpos2);
101       
102       
103/*
104          incrementalWordSimilarity(sc1, sc2, i, j)
105        + incrementalProbabilityScore(sc1, sc2) * multiDeltaScore
106        + incrementalStackingScore(sc1, sc2) * multiDeltaStacking;
107*/
108
109    return score*sc1.GetNumSeq()*sc2.GetNumSeq();
110}
111
112float
113Globaldp::
114incrementalProbabilityScore(const StemCandidate &sc1, const StemCandidate &sc2)
115{
116    int wordLength = WORDLENGTH;
117    return sc1.GetBaseScore(wordLength-1) + sc2.GetBaseScore(wordLength-1);
118}
119
120float
121Globaldp::
122incrementalStackingScore(const StemCandidate &sc1, const StemCandidate &sc2)
123{
124    return - (sc1.GetStacking() + sc2.GetStacking());
125}
126
127float
128Globaldp::
129incrementalWordSimilarity(const StemCandidate &sc1, const StemCandidate &sc2, int i, int j)
130{
131    int numSeq1 = sc1.GetNumSeq();
132    int numSeq2 = sc2.GetNumSeq();
133
134    float score = 0;
135
136    for(int k = 0; k < 16; k++) {
137        for(int l = 0; l < 16; l++) {
138            score += countContConp1[k][i]*countContConp2[l][j]*ribosum_matrix[k][l];
139        }
140    }
141
142    return score/(numSeq1*numSeq2);
143}
144
145float 
146Globaldp::
147scorePSCPair(const StemCandidate &sc1, const StemCandidate &sc2)
148{
149
150
151    int wordLength = WORDLENGTH;
152    float score = 0;
153   
154    int pos1, rvpos1;
155    if (sc1.GetDistance() > 0) {
156        pos1   = sc1.GetPosition();
157        rvpos1 = sc1.GetRvposition();
158    }
159    else {
160        pos1   = sc1.GetRvposition();
161        rvpos1 = sc1.GetPosition();
162    }
163
164    int pos2, rvpos2;
165    if (sc2.GetDistance() > 0) {
166        pos2   = sc2.GetPosition();
167        rvpos2 = sc2.GetRvposition();
168    }
169    else {
170        pos2   = sc2.GetRvposition();
171        rvpos2 = sc2.GetPosition();
172    }
173   
174    for (int k = 0; k < wordLength; k++) {
175        score +=
176              posterior[sc1.GetPosition() + k][sc2.GetPosition() + k]
177            * profileBPPMat1->ref(pos1 + k, rvpos1 + wordLength - 1 - k)
178//          * sc1.GetBaseScore(k)
179            * posterior[sc1.GetRvposition() + wordLength - 1 - k][sc2.GetRvposition() + wordLength - 1 - k]
180//          * sc2.GetBaseScore(k);
181            * profileBPPMat2->ref(pos2 + k, rvpos2 + wordLength - 1 - k);
182    }
183    // validation code
184    /*
185    if (profileBPPMat1->ref(pos1 , rvpos1 + wordLength - 1) != sc1.GetBaseScore(0)) {
186        cout << "1 " << profileBPPMat1->ref(pos1, rvpos1 + wordLength - 1) << " " << sc1.GetBaseScore(0) << endl;
187    }
188    if (profileBPPMat2->ref(pos2, rvpos2 + wordLength - 1) != sc2.GetBaseScore(0)) {
189        cout << profileBPPMat2->ref(pos2, rvpos2 + wordLength - 1) << " " << sc2.GetBaseScore(0) << endl;
190    }
191    if (profileBPPMat1->ref(pos1 + 1, rvpos1 + wordLength - 1 - 1) != sc1.GetBaseScore(1)) {
192        cout << "1 " << profileBPPMat1->ref(pos1 + 1, rvpos1 + wordLength - 1 - 1) << " " << sc1.GetBaseScore(1) << endl;
193    }
194    if (profileBPPMat2->ref(pos2 + 1, rvpos2 + wordLength - 1 - 1) != sc2.GetBaseScore(1)) {
195        cout << profileBPPMat2->ref(pos2 + 1, rvpos2 + wordLength - 1 - 1) << " " << sc2.GetBaseScore(1) << endl;
196        }*/
197
198//    cout << sc1.GetPosition() << " " << sc1.GetRvposition() << " " << sc1.GetDistance() << endl;
199//    cout << sc2.GetPosition() << " " << sc2.GetRvposition() << " " << sc2.GetDistance() << endl;
200/*
201          wordSimilarity(sc1, sc2, i, j)
202        + probabilityScore(sc1, sc2) * multiScore
203        + stackingScore(sc1, sc2) * multiStacking
204
205*/ 
206/*
207    if (sc1.GetDistance() < 0 && sc2.GetDistance() < 0) {
208        cout << "2:" << i << " " << j << " " << sc1.GetDistance() << " " << sc2.GetDistance() << endl;
209        cout << sc1.GetPosition() << " " << sc1.GetRvposition() << " " << sc2.GetPosition() << " " << sc2.GetRvposition() << endl;
210        cout << "2:score" << score << endl;
211
212    }
213*/
214    return score*sc1.GetNumSeq()*sc2.GetNumSeq();
215}
216
217float
218Globaldp::
219wordSimilarity(const StemCandidate &sc1, const StemCandidate &sc2, int i, int j)
220{
221    int wordLength = WORDLENGTH;
222
223    int numSeq1 = sc1.GetNumSeq();
224    int numSeq2 = sc2.GetNumSeq();
225
226    float score = 0;
227
228    for(int k = 0; k < 16; k++) {
229        for(int l = 0; l < 16; l++) {
230            for(int iter = 0; iter < wordLength; iter++) {
231                score += countConp1[iter][k][i]*countConp2[iter][l][j]*ribosum_matrix[k][l];
232            }
233        }
234    }
235
236    return score/(numSeq1*numSeq2);
237}
238
239int
240Globaldp::
241returnBasePairType(char s, char r)
242{
243    if  (s == 'A') {
244        if      (r == 'A') return 0;
245        else if (r == 'C') return 1;
246        else if (r == 'G') return 2;
247        else if (r == 'U') return 3;
248    }
249    else if (s == 'C') {
250        if      (r == 'A') return 4;
251        else if (r == 'C') return 5;
252        else if (r == 'G') return 6;
253        else if (r == 'U') return 7;
254    }
255    else if (s == 'G') {
256        if      (r == 'A') return 8;
257        else if (r == 'C') return 9;
258        else if (r == 'G') return 10;
259        else if (r == 'U') return 11;
260    }
261    else if (s == 'U') {
262        if      (r == 'A') return 12;
263        else if (r == 'C') return 13;
264        else if (r == 'G') return 14;
265        else if (r == 'U') return 15;
266    }
267
268    return 16;
269}
270
271float
272Globaldp::
273probabilityScore(const StemCandidate &sc1, const StemCandidate &sc2)
274{
275    return sc1.GetScore() + sc2.GetScore();
276}
277
278float
279Globaldp::
280stackingScore(const StemCandidate &sc1, const StemCandidate &sc2)
281{
282    return - (sc1.GetStemStacking() + sc2.GetStemStacking());
283}
284
285float
286Globaldp::
287distancePenalty(const StemCandidate &sc1, const StemCandidate &sc2)
288{
289    return std::sqrt((float)std::abs(sc1.GetDistance() - sc2.GetDistance()));
290}
291
292float
293Globaldp::
294Run()
295{
296    Initialize();
297    DP();
298    float score = traceBack();
299   
300    // printMatch();
301    //cout << "score=" << score << endl;
302    return score;
303}
304
305void
306Globaldp::
307Initialize()
308{
309    int nPscs1 = pscs1->size();
310    int nPscs2 = pscs2->size();
311
312
313    for(int i = 0; i < nPscs1; i++) {
314        for(int j = 0; j < nPscs2; j++) {
315            VM[i][j] = 0;
316            VG[i][j] = 0;
317        }
318    }
319
320    VM[0][0] = 0;
321    VG[0][0] = IMPOSSIBLE;
322
323    for (int i = 1; i < nPscs1; i++) {
324        VM[i][0] = IMPOSSIBLE; VG[i][0] = 0;
325    }
326    for (int j = 1; j < nPscs2; j++) {
327        VM[0][j] = IMPOSSIBLE; VG[0][j] = 0;
328    }
329
330    for (int i = 0; i < nPscs1; i++) {
331        for (int j = 0; j < nPscs2; j++) {
332            traceMI[i][j] = 0; traceMJ[i][j] = 0;
333            traceGI[i][j] = 0; traceGJ[i][j] = 0;
334        }
335    }
336
337    int wordLength = WORDLENGTH;
338    int size1   = pscs1->size();
339    int size2   = pscs2->size();
340
341    for(int i = 0; i < wordLength; i++) {
342        for(int j = 0; j < 17; j++) {
343            for(int k = 0; k < (int)pscs1->size(); k++) {
344                countConp1[i][j][k] = 0;
345            }
346        }
347    }
348    for(int i = 0; i < wordLength; i++) {
349        for(int j = 0; j < 17; j++) {
350            for(int k = 0; k < (int)pscs2->size(); k++) {
351                countConp2[i][j][k] = 0;
352            }
353        }
354    }
355    for(int i = 0; i < 17; i++) {
356        for(int j = 0; j < (int)pscs1->size()+1; j++) {
357            countContConp1[i][j] = 0;
358        }
359    }
360    for(int i = 0; i < 17; i++) {
361        for(int j = 0; j < (int)pscs2->size()+1; j++) {
362            countContConp2[i][j] = 0;
363        }
364    }
365
366    for(int iter = 1; iter < size1; iter++) {
367
368        const StemCandidate &sc1 = pscs1->at(iter);
369        int numSeq1 = sc1.GetNumSeq();
370        for (int i = 0; i < wordLength; i++) {
371       
372            for(int k = 0; k < numSeq1; k++) {
373//              const Sequence *seq = seqs1->GetSequence(k);
374                string substr = sc1.GetSubstr(k);
375                string rvstr  = sc1.GetRvstr(k);
376           
377                char sChar = substr[i];
378                char rChar = rvstr[wordLength - 1 -i];
379           
380                int type = returnBasePairType(sChar, rChar);
381                ++countConp1[i][type][iter];
382            }
383        }
384        for(int k = 0; k < numSeq1; k++) {
385//          const Sequence *seq = seqs1->GetSequence(k);
386            string substr = sc1.GetSubstr(k);
387            string rvstr  = sc1.GetRvstr(k);
388
389            char sChar = substr[wordLength-1];
390            char rChar = rvstr[0];
391           
392            int type = returnBasePairType(sChar, rChar);
393            ++countContConp1[type][iter];
394
395        }
396    }
397    for(int iter = 1; iter < size2; iter++) {
398        const StemCandidate &sc2 = pscs2->at(iter);
399        int numSeq2 = sc2.GetNumSeq();
400        for (int i = 0; i < wordLength; i++) {
401       
402            for(int k = 0; k < numSeq2; k++) {
403//              const Sequence *seq = seqs2->GetSequence(k);
404                string substr = sc2.GetSubstr(k);
405                string rvstr  = sc2.GetRvstr(k);
406           
407                char sChar = substr[i];
408                char rChar = rvstr[wordLength - 1 -i];
409           
410                int type = returnBasePairType(sChar, rChar);
411                ++countConp2[i][type][iter];
412            }
413        }
414        for(int k = 0; k < numSeq2; k++) {
415//          const Sequence *seq = seqs2->GetSequence(k);
416            string substr = sc2.GetSubstr(k);
417            string rvstr  = sc2.GetRvstr(k);
418
419            char sChar = substr[wordLength-1];
420            char rChar = rvstr[0];
421           
422            int type = returnBasePairType(sChar, rChar);
423            ++countContConp2[type][iter];
424
425        }
426    }
427    profileBPPMat1 = makeProfileBPPMatrix(seqs1);
428    profileBPPMat2 = makeProfileBPPMatrix(seqs2);
429}
430
431void
432Globaldp::
433DP()
434{
435    int nPscs1 = pscs1->size();
436    int nPscs2 = pscs2->size();
437   
438    for(int i = 1; i < nPscs1; i++) {
439        const StemCandidate &sc1 = pscs1->at(i);
440
441        for(int j = 1; j < nPscs2; j++) {
442            const StemCandidate &sc2 = pscs2->at(j);
443           
444            float max = IMPOSSIBLE;
445            int traceI = 0; 
446            int traceJ = 0;
447            int alpha = sc1.GetContPos(), beta = sc2.GetContPos();
448            if( (alpha > 0) && (beta > 0) ) {
449                max = VM[alpha][beta] + incrementalScorePSCPair(sc1, sc2);
450                traceI = alpha;
451                traceJ = beta;
452            }
453           
454            float similarity = scorePSCPair(sc1, sc2);
455            int p = sc1.GetBeforePos(), q = sc2.GetBeforePos();
456            float tmpScore = VM[p][q] + similarity;
457//          cout << i << " " << j << " " << p << " " << q << " " << tmpScore << " " << VM[p][q] << " " << similarity << " " << endl;
458               
459            if(max <= tmpScore) {
460                max = tmpScore;
461                traceI = p;
462                traceJ = q;
463            }
464
465            tmpScore = VG[p][q] + similarity;
466            if(max <= tmpScore) {
467                max = tmpScore;
468                traceI = traceGI[p][q];
469                traceJ = traceGJ[p][q];
470            }
471           
472            VM[i][j]      = max;
473            traceMI[i][j] = traceI;
474            traceMJ[i][j] = traceJ;
475           
476            max = VM[i][j-1];
477            traceI   = i;
478            traceJ   = j-1;
479
480            tmpScore = VM[i-1][j];
481            if(max <= tmpScore) {
482                max    = tmpScore;
483                traceI = i-1;
484                traceJ = j;
485            }
486
487            tmpScore = VG[i-1][j];
488            if(max <= tmpScore) {
489                max    = tmpScore;
490                traceI = traceGI[i-1][j];
491                traceJ = traceGJ[i-1][j];
492            }
493           
494            tmpScore = VG[i][j-1];
495            if(max <= tmpScore) {
496                max = tmpScore;
497                traceI = traceGI[i][j-1];
498                traceJ = traceGJ[i][j-1];
499            }
500
501            VG[i][j]      = max;
502            traceGI[i][j] = traceI;
503            traceGJ[i][j] = traceJ;
504        }
505    }
506}
507
508float
509Globaldp::
510traceBack()
511{
512    int nPscs1 = pscs1->size();
513    int nPscs2 = pscs2->size();
514   
515    float max = IMPOSSIBLE;
516    int traceI, traceJ;
517    for (int i = 1; i < nPscs1; i++) {
518        for (int j = 1; j < nPscs2; j++) {
519            if(max < VM[i][j]) {
520                max = VM[i][j];
521                traceI = i;
522                traceJ = j;
523            }
524        }
525    }
526
527    while ( (traceI != 0) &&  (traceJ != 0) ) {
528        matchPSCS1->push_back(traceI); matchPSCS2->push_back(traceJ);
529        int tmpI = traceMI[traceI][traceJ];
530        int tmpJ = traceMJ[traceI][traceJ];
531        traceI = tmpI; traceJ = tmpJ;
532    }
533   
534    if(traceI != 0 && traceJ != 0) {
535        std::cout << "error in trace back of pscs:" << traceI << " " << traceJ << endl;
536    }
537
538    return max;
539}
540
541void
542Globaldp::
543printMatch()
544{
545    int size = matchPSCS1->size();
546    for(int iter = 0; iter < size; iter++) {
547        int i = matchPSCS1->at(iter);
548        int j = matchPSCS2->at(iter);
549       
550        const StemCandidate &sc1 = pscs1->at(i);
551        const StemCandidate &sc2 = pscs2->at(j);
552
553        std::cout << iter << "\t" << sc1.GetPosition()  << "\t" << sc1.GetRvposition() << "\t" << sc2.GetPosition() << "\t" << sc2.GetRvposition() << endl;
554    }
555
556}
557}
Note: See TracBrowser for help on using the repository browser.