source: tags/arb-6.0/GDE/MAFFT/mafft-7.055-with-extensions/extensions/mxscarna_src/seq2scs.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: 19.2 KB
Line 
1///////////////////////////////////////////////////////////////
2// seq2scs.cpp
3//
4// make SCS(Stem Candidate Sequence) from the profile
5//////////////////////////////////////////////////////////////
6
7#include "scarna.hpp"
8#include "SafeVector.h"
9#include "StemCandidate.hpp"
10#include "Sequence.h"
11#include "MultiSequence.h"
12#include "BPPMatrix.hpp"
13#include "nrutil.h"
14#include <vector>
15#include <algorithm>
16#include <stdio.h>
17#include <cstring>
18#include <stdio.h>
19#include <math.h>
20
21using namespace std;
22using namespace::MXSCARNA;
23
24// for alipfold
25/*
26#include "utils.h"
27#include "fold_vars.h"
28#include "fold.h"
29#include "part_func.h"
30#include "inverse.h"
31#include "RNAstruct.h"
32#include "treedist.h"
33#include "stringdist.h"
34#include "profiledist.h"
35#include "alifold.h"
36#include "aln_util.h"
37#include "dist_vars.h"
38*/
39double Stacking_Energy[36] ={
40 -0.9,-2.1,-1.7,-0.5,-0.9,-1.0,
41 -1.8,-2.9,-2.0,-1.2,-1.7,-1.9,
42 -2.3,-3.4,-2.9,-1.4,-2.1,-2.1,
43 -1.1,-2.1,-1.9,-0.4,-1.0,1.5,
44 -1.1,-2.3,-1.8,-0.8,-0.9,-1.1,
45 -0.8,-1.4,-1.2,-0.2,-0.5,-0.4 };
46
47static Trimat<float>* makeProfileBPPMatrix(MultiSequence *Sequences, SafeVector<BPPMatrix*> &BPPMatrices);
48static int countSCS(MultiSequence *Sequences, Trimat<float>* consBppMat, int BandWidth);
49static std::vector<StemCandidate>* makeProfileScs(std::vector<StemCandidate> *pscs, MultiSequence *Sequences, Trimat<float>* consBppMat, int BandWidth);
50static void printScs(std::vector<StemCandidate> *pscs);
51static std::vector<StemCandidate>* doubleScs(std::vector<StemCandidate> *pscs);
52static std::vector<StemCandidate>* findRelations(std::vector<StemCandidate> *pscs);
53static std::vector<StemCandidate>* findCorresponding(std::vector<StemCandidate>* pscs);
54static std::vector<StemCandidate>* calculateStackingEnergy(std::vector<StemCandidate>* pscs);
55
56//float alipf_fold(char **sequences, char *structure, pair_info **pi);
57
58struct SortCmp {
59    bool operator()(const StemCandidate &sc1, const StemCandidate &sc2) const {
60        if      (sc1.GetPosition() > sc2.GetPosition()) return false;
61        else if (sc1.GetPosition() < sc2.GetPosition()) return true;
62        else if (sc1.GetDistance() > sc2.GetDistance()) return false;
63        else return true;
64    }
65};
66
67
68vector<StemCandidate>*
69seq2scs(MultiSequence *Sequences, SafeVector<BPPMatrix*> &BPPMatrices, int BandWidth)
70{
71
72    Trimat<float> *consBppMat = makeProfileBPPMatrix(Sequences, BPPMatrices);
73
74    int numberScs = countSCS(Sequences, consBppMat, BandWidth);
75
76    std::vector<StemCandidate> *pscs = new std::vector<StemCandidate>(); // Profile Stem Candidate Sequence
77//    cout << "numberScs=" << numberScs << endl;
78    pscs->resize(numberScs+1);
79
80    pscs = makeProfileScs(pscs, Sequences, consBppMat, BandWidth);
81
82    pscs = doubleScs(pscs);
83
84    std::vector<StemCandidate>::iterator startIter = pscs->begin();
85    std::vector<StemCandidate>::iterator endIter   = pscs->end();
86    ++startIter;
87    std::sort(startIter, endIter, SortCmp());
88   
89//    printScs(pscs);
90    pscs = findRelations(pscs);
91
92    pscs = findCorresponding(pscs);
93
94    pscs = calculateStackingEnergy(pscs);
95
96//    findStemRelation()
97   
98//    exit(1);
99    delete consBppMat;
100
101    return pscs;
102}
103
104static Trimat<float>*
105makeProfileBPPMatrix(MultiSequence *Sequences, SafeVector<BPPMatrix*> &BPPMatrices)
106{
107    int length = Sequences->GetSequence(0)->GetLength();
108//    float thr  = BaseProbThreshold;
109    Trimat<float> *consBppMat = new Trimat<float>(length + 1);
110    fill(consBppMat->begin(), consBppMat->end(), 0);
111
112// gabage
113//    for(int i = 0; i <= length; i++)
114//      for(int j = i; j <= length; j++)
115//          cout << "i=" << i << " j=" << j << " " << consBppMat->ref(i,j) << endl;
116//          consBppMat->ref(i,j) = 0;
117
118
119    int number = Sequences->GetNumSequences();
120//    if( number == 1) {
121    for(int seqNum = 0; seqNum < number; seqNum++) {
122        SafeVector<int> *tmpMap = Sequences->GetSequence(seqNum)->GetMappingNumber();
123        int label = Sequences->GetSequence(seqNum)->GetLabel();
124        BPPMatrix *tmpBppMatrix = BPPMatrices[label];
125           
126        for(int i = 1; i <= length ; i++) {
127            int originI = tmpMap->at(i);
128            for(int j = i + 3; j <= length; j++) {
129                int originJ = tmpMap->at(j);
130                if(originI != 0 && originJ != 0) {
131                    float tmpProb = tmpBppMatrix->GetProb(originI, originJ);
132
133//                  if(tmpProb >= thr) {
134                        consBppMat->ref(i, j) += tmpProb;
135//                      cout << i << " " << j << " " << consBppMat->ref(i,j) << endl;
136//                  }
137                }
138            }
139        }
140    }
141
142        /* compute the mean of base pairing probability  */
143    for(int i = 1; i <= length; i++) {
144        for(int j = i + 3; j <= length; j++) {
145            consBppMat->ref(i,j) = consBppMat->ref(i,j)/(float)number;
146            //consBppMat->ref(i,j) = std::sqrt[number](consBppMat->ref(i,j));
147            //      cout << i << " " << j <<  " " << consBppMat->ref(i,j) << endl;
148        }
149    }
150
151
152/*
153    else {
154        char **Seqs;
155        Seqs = (char **) malloc(sizeof(char*) * number);
156
157        for(int i = 0; i < number; i++) {
158            Seqs[i] = (char *) malloc(sizeof(char) * (length + 1));
159            for(int j = 1; j <= length; j++) {
160                Seqs[i][j-1] = Sequences->GetSequence(i)->GetPosition(j);
161            }
162            Seqs[i][length] = '\0';
163        }
164
165       
166        char *structure = NULL;
167        pair_info *pi;
168           
169        alipf_fold(Seqs, structure, &pi, number);
170
171        for(int iter = 0; iter < length; iter++) {
172            if(pi[iter].i == 0) break;
173            consBppMat->ref(pi[iter].i, pi[iter].j) = pi[iter].p;
174        }
175
176        for(int i = 0; i < number; i++) {
177            free (Seqs[i]);
178        }
179        free (Seqs);
180       
181//      free_alifold_arrays(void);
182    }
183*/ 
184
185    return consBppMat;
186}
187 
188static int 
189countSCS(MultiSequence *Sequences, Trimat<float> *consBppMat, int BandWidth)
190{
191
192    int length = Sequences->GetSequence(0)->GetLength();
193    int Word_Length = scsLength;
194
195    int i, j, k, s;
196    int count;
197    int sum;
198
199    sum = 0;
200    for(k = 1; k <= length; k++) {
201        count = 0;
202
203        for(i = k, j = k; i >= 1 && j <= (k + BandWidth - 1) && j <= length; i--, j++) {
204            if(consBppMat->ref(i, j) >= BaseProbThreshold) {
205                count++;
206            }
207            else if(count >= Word_Length) {
208                for(s = 0; s <= count - Word_Length; s++) 
209                  sum++;
210
211                count = 0;
212            }
213            else {
214                count = 0;
215            }
216            if(consBppMat->ref(i, j) >= BaseProbThreshold && count >= Word_Length && (i == 1 || j == length || j == (k + BandWidth - 1))) {
217                for(s = 0; s <= count - Word_Length; s++) 
218                    sum++;
219
220                count = 0;
221            }
222        }
223   
224        count = 0;
225        for(i = k, j = k + 1; i >= 1 && j <= (k + BandWidth - 1) && j <= length; i--, j++) {
226            if(consBppMat->ref(i, j) >= BaseProbThreshold) {
227                count++;
228            }
229            else if(count >= Word_Length) {
230                for(s = 0; s <= count - Word_Length; s++) sum++;
231
232                count = 0;
233            }
234            else {
235                count = 0;
236            }
237
238            if(consBppMat->ref(i, j) >= BaseProbThreshold && count >= Word_Length && (i == 1 || j == length || j == (k + BandWidth - 1))) {
239                for(s = 0; s <= count - Word_Length; s++) sum++;
240
241                count = 0;
242            }
243        }
244    }
245
246    return 2 * sum;
247}
248
249static std::vector<StemCandidate>* 
250makeProfileScs(std::vector<StemCandidate> *pscs, MultiSequence *Sequences, Trimat<float>* consBppMat, int BandWidth)
251{
252
253    int length = Sequences->GetSequence(0)->GetLength();
254    int Word_Length = scsLength;
255    float Thr = BaseProbThreshold;
256    int i, j, k, s, t, m, n, l;
257    int count;
258    int sum;
259   
260    sum = 0;
261    for(k = 1; k <= length; k++) {
262        count = 0;
263        for(i = k, j = k; i >= 1 && j <= (k + BandWidth - 1) && j <= length; i--, j++) {
264            if(consBppMat->ref(i,j) >= Thr) {
265                count++;
266            }
267            else if(count >= Word_Length) {
268                for(s = 0; s <= count- Word_Length; s++) {
269                    sum++;
270                    pscs->at(sum).SetLength(Word_Length);
271                    pscs->at(sum).SetPosition(i+1+s);
272                    pscs->at(sum).SetRvposition(j-count+(count-Word_Length-s));
273                    pscs->at(sum).SetDistance((j-count+count-Word_Length-s) - (i+1+s+Word_Length));
274                    pscs->at(sum).SetNumSeq(Sequences->GetNumSequences());
275                    pscs->at(sum).SetNumSubstr(Sequences->GetNumSequences());
276                    pscs->at(sum).SetNumRvstr(Sequences->GetNumSequences());
277                    for(m = i + 1 + s, n = j - count + (count-Word_Length-s), l = j - 1 - s, t = 0; n < j-s; m++, n++, l--, t++) {
278                        for(int num = 0; num < Sequences->GetNumSequences(); num++) {
279                            Sequence *seq = Sequences->GetSequence(num);
280//                          cout << num << "; " << m << ":" << seq->GetPosition(m) << " " << n << ":" << seq->GetPosition(n) << endl;
281                            pscs->at(sum).AddSubstr(num, seq->GetPosition(m));
282                            pscs->at(sum).AddRvstr(num, seq->GetPosition(n));
283                        }
284                        //          assert(pr[iindx[m]-l] > Thr);
285//                      cout << "prob=" << consBppMat->ref(m,l) << endl;
286                        pscs->at(sum).AddScore(consBppMat->ref(m,l));
287                        pscs->at(sum).AddBaseScore(consBppMat->ref(m, l));
288                    }
289                    for(int num = 0; num < Sequences->GetNumSequences(); num++) {
290                        pscs->at(sum).AddSubstr(num, '\0'); 
291                        pscs->at(sum).AddRvstr(num, '\0');
292                    }
293                }
294                count = 0;
295            }
296            else {
297                count = 0;
298            }
299            if(consBppMat->ref(i,j) >= Thr && count >= Word_Length && (i == 1 || j == length || j == (k + BandWidth - 1))) {
300                for(s = 0; s <= count- Word_Length; s++) {
301                    sum++;
302                    pscs->at(sum).SetLength(Word_Length);
303                    pscs->at(sum).SetPosition(i+s);
304                    pscs->at(sum).SetRvposition(j-count+1+(count-Word_Length-s));
305                    pscs->at(sum).SetDistance((j-count+1+count-Word_Length-s) - (i+s+Word_Length));
306                    pscs->at(sum).SetNumSeq(Sequences->GetNumSequences());
307                    pscs->at(sum).SetNumSubstr(Sequences->GetNumSequences());
308                    pscs->at(sum).SetNumRvstr(Sequences->GetNumSequences());
309                    for(m = i + s, n = j - count + 1 + (count-Word_Length-s), l = j - s, t = 0; n <= j-s; m++, n++, l--, t++) {
310                        for(int num = 0; num < Sequences->GetNumSequences(); num++) {
311                            Sequence *seq = Sequences->GetSequence(num);
312                            pscs->at(sum).AddSubstr(num, seq->GetPosition(m));
313                            pscs->at(sum).AddRvstr(num, seq->GetPosition(n));
314                        }
315
316                        pscs->at(sum).AddScore(consBppMat->ref(m,l));
317                        pscs->at(sum).AddBaseScore(consBppMat->ref(m, l));
318                    }
319                    for(int num = 0; num < Sequences->GetNumSequences(); num++) {
320                        pscs->at(sum).AddSubstr(num, '\0'); 
321                        pscs->at(sum).AddRvstr(num, '\0');
322                    }
323                }
324                count = 0;
325            }
326        }
327        count = 0;
328        for(i = k, j = k + 1; i >= 1 && j <= (k + BandWidth - 1) && j <= length; i--, j++) {
329            if(consBppMat->ref(i,j) >= Thr) {
330                count++;
331            }
332            else if(count >= Word_Length) {
333                for(s = 0; s <= count- Word_Length; s++) {
334                    sum++;
335                    pscs->at(sum).SetLength(Word_Length);
336                    pscs->at(sum).SetPosition(i+1+s);
337                    pscs->at(sum).SetRvposition( j-count+(count-Word_Length-s));
338                    pscs->at(sum).SetDistance((j-count+count-Word_Length-s) - (i+1+s+Word_Length));
339                    pscs->at(sum).SetNumSeq(Sequences->GetNumSequences());
340                    pscs->at(sum).SetNumSubstr(Sequences->GetNumSequences());
341                    pscs->at(sum).SetNumRvstr(Sequences->GetNumSequences());
342                    for(m = i + 1 + s, n = j - count + (count-Word_Length-s), l = j - 1 - s, t = 0; n < j-s; m++, n++, l--, t++) {
343                        for(int num = 0; num < Sequences->GetNumSequences(); num++) {
344                            Sequence *seq = Sequences->GetSequence(num);
345                            pscs->at(sum).AddSubstr(num, seq->GetPosition(m));
346                            pscs->at(sum).AddRvstr(num, seq->GetPosition(n));
347                        }
348
349                        pscs->at(sum).AddScore(consBppMat->ref(m,l));
350                        pscs->at(sum).AddBaseScore(consBppMat->ref(m, l));
351                    }
352                    for(int num = 0; num < Sequences->GetNumSequences(); num++) {
353                        pscs->at(sum).AddSubstr(num, '\0'); 
354                        pscs->at(sum).AddRvstr(num, '\0');
355                    }
356                }
357                count = 0;
358            }
359            else {
360                count = 0;
361            }
362            if(consBppMat->ref(i,j) >= Thr && count >= Word_Length && (i == 1 || j == length || j == (k + BandWidth - 1))) {
363                for(s = 0; s <= count - Word_Length; s++) {
364                    sum++;
365                    pscs->at(sum).SetLength(Word_Length);
366                    pscs->at(sum).SetPosition(i+s);
367                    pscs->at(sum).SetRvposition(j-count+1+(count-Word_Length-s));
368                    pscs->at(sum).SetDistance((j-count+1+count-Word_Length-s) - (i+s+Word_Length));
369//                  pscs->at(sum).SetDistance((j-count+count-Word_Length-s) - (i+1+s+Word_Length));
370                    pscs->at(sum).SetNumSeq(Sequences->GetNumSequences());
371                    pscs->at(sum).SetNumSubstr(Sequences->GetNumSequences());
372                    pscs->at(sum).SetNumRvstr(Sequences->GetNumSequences());
373                    for(m = i + s, n = j - count + 1 + (count-Word_Length-s), l = j - s, t=0; n <= j-s; m++, n++, l--, t++) {
374                        for(int num = 0; num < Sequences->GetNumSequences(); num++) {
375                            Sequence *seq = Sequences->GetSequence(num);
376                            pscs->at(sum).AddSubstr(num, seq->GetPosition(m));
377                            pscs->at(sum).AddRvstr(num, seq->GetPosition(n));
378                        }
379
380                        pscs->at(sum).AddScore(consBppMat->ref(m,l));
381                        pscs->at(sum).AddBaseScore(consBppMat->ref(m, l));
382                    }
383                    for(int num = 0; num < Sequences->GetNumSequences(); num++) {
384                        pscs->at(sum).AddSubstr(num, '\0'); 
385                        pscs->at(sum).AddRvstr(num, '\0');
386                    }
387                }
388                count = 0;
389            }
390        }
391    }
392
393    return pscs;
394}
395
396static std::vector<StemCandidate>* 
397doubleScs(std::vector<StemCandidate> *pscs)
398{
399    int num = pscs->size()/2;
400   
401    for(int i = 1; i <= num; i++) {
402        int latter = num + i;
403        //cout << i << " " << latter << endl;
404        StemCandidate &tmpScs = pscs->at(i);
405        pscs->at(latter).SetLength(tmpScs.GetLength());
406        pscs->at(latter).SetPosition(tmpScs.GetRvposition());
407        pscs->at(latter).SetRvposition(tmpScs.GetPosition());
408        pscs->at(latter).SetDistance(-tmpScs.GetDistance());
409        pscs->at(latter).SetNumSeq(tmpScs.GetNumSeq());
410        pscs->at(latter).SetNumSubstr(tmpScs.GetNumSeq());
411        pscs->at(latter).SetNumRvstr(tmpScs.GetNumSeq());
412
413        pscs->at(latter).SetScore(tmpScs.GetScore());
414        for(int num = 0; num < tmpScs.GetNumSeq(); num++) {
415            string tmpSubstr = tmpScs.GetSubstr(num);
416            string tmpRvstr  = tmpScs.GetRvstr(num);
417
418            for(int k = 0; k < tmpScs.GetLength(); k++) {
419                pscs->at(latter).AddSubstr(num, tmpSubstr[k]);
420                pscs->at(latter).AddRvstr(num, tmpRvstr[k]);
421            }
422        }
423        for(int k = 0; k < tmpScs.GetLength(); k++) {
424            pscs->at(latter).AddBaseScore(tmpScs.GetBaseScore(k));
425        }
426    }
427    return pscs;
428}
429
430
431static void 
432printScs(std::vector<StemCandidate> *pscs)
433{
434    int num = pscs->size();
435//    std::cout << "size = " << num << endl;
436    for(int i = 1; i < num; i++) {
437        StemCandidate &sc = pscs->at(i);
438       
439        std::cout << i << "\t" << sc.GetLength() << "\t" << sc.GetPosition() << "\t" << 
440            sc.GetRvposition() << "\t" << sc.GetDistance() << "\t" << sc.GetNumSeq() << 
441            "\t" << sc.GetScore() << "\t" << sc.GetContPos() << "\t" << sc.GetBeforePos() << 
442            "\t" << sc.GetRvscnumber() << "\t" << sc.GetStacking() << "\t" << sc.GetStemStacking() << 
443            "\t" << sc.GetBaseScore(0) << "\t" << sc.GetBaseScore(1) << endl;
444        cout << "substr:" << endl;
445        for(int k = 0; k < sc.GetNumSeq(); k++) {
446            cout << k << "\t" << sc.GetSubstr(k) << "\t" << sc.GetRvstr(k) << "\t" << endl;
447        }
448       
449    }
450
451}
452
453static std::vector<StemCandidate>*
454findRelations(std::vector<StemCandidate> *pscs)
455{
456    int num = pscs->size();
457   
458    for(int i = 1; i < num; i++) {
459        int pt = i-1; 
460        StemCandidate &sc = pscs->at(i);
461        sc.SetContPos(-1);
462        while(sc.GetPosition() == pscs->at(pt).GetPosition()) { pt--; }
463       
464        while((sc.GetPosition() == pscs->at(pt).GetPosition() + 1) && (pt > 0)) {
465            if(sc.GetRvposition() == pscs->at(pt).GetRvposition() - 1) {
466                sc.SetContPos(pt);
467                break;
468            }
469            --pt;
470        }
471        while((sc.GetPosition() < pscs->at(pt).GetPosition() + pscs->at(pt).GetLength())&&(pt > 0)) { pt--; }
472        sc.SetBeforePos(pt);
473    }
474   
475    return pscs;
476}
477
478static std::vector<StemCandidate>*
479findCorresponding(std::vector<StemCandidate>* pscs)
480{
481    int num = pscs->size();
482   
483    for(int i = 1; i < num; i++) { pscs->at(i).SetRvscnumber(0); }
484   
485    for(int i = 1; i < num; i++) {
486        if(pscs->at(i).GetDistance() > 0) {
487            for(int j = i + 1; j < num; j++) {
488                if ( (pscs->at(j).GetPosition() == pscs->at(i).GetRvposition())
489                     && (pscs->at(i).GetPosition() == pscs->at(j).GetRvposition()) ) {
490                    pscs->at(i).SetRvscnumber(j);
491                    pscs->at(j).SetRvscnumber(i);
492                    break;
493                }
494            }
495        }
496        if(pscs->at(i).GetRvscnumber() == 0) {
497            std::cerr << "error in findCorresponding" << " i=" << i << endl;
498//          exit(1);
499        }
500    }
501
502    return pscs;
503}
504
505static std::vector<StemCandidate>*
506calculateStackingEnergy(std::vector<StemCandidate>* pscs)
507{
508    int num = pscs->size();
509    int wordLength = scsLength;
510
511    for(int i = 1; i < num; i++) {
512        StemCandidate &scI = pscs->at(i);
513        int j = pscs->at(i).GetContPos();
514       
515        if(j > 0) {
516
517            StemCandidate &scJ = pscs->at(j);
518            float stacking = 0;
519            int profNum =  scJ.GetNumSeq();
520            for(int k = 0; k < profNum; k++) {
521                string substr = scJ.GetSubstr(k);
522                string rvstr  = scJ.GetRvstr(k);
523                int index = 0;
524                switch(substr[wordLength - 1]) {
525                        case 'A': if(     rvstr[0]=='U' ) {index += 0;}
526                                  else{ index = -1000; }
527                                  break;
528                        case 'C': if(     rvstr[0]=='G' ) {index += 6;}
529                                  else{ index = -1000; }
530                                  break;
531                        case 'G': if(     rvstr[0]=='C'){index += 12;}
532                                  else if(rvstr[0]=='U'){index += 18;}
533                                  else{ index = -1000; }
534                                  break;
535                        case 'U': if(     rvstr[0]=='A'){index += 24;}
536                                  else if(rvstr[0]=='G'){index += 30;}
537                                  else{ index = - 1000; }
538                                  break;
539                }
540                substr = scI.GetSubstr(k);
541                rvstr  = scI.GetRvstr(k);
542                switch(substr[wordLength - 1]){
543                        case 'A': if(     rvstr[0]=='U'){index += 0;}
544                                  else{ index = -1000; }
545                                  break;
546                        case 'C': if(     rvstr[0]=='G'){index += 1;}
547                                  else{ index = -1000; }
548                                  break;
549                        case 'G': if(     rvstr[0]=='C'){index += 2;}
550                                  else if(rvstr[0]=='U'){index += 3;}
551                                  else{ index = -1000; }
552                                  break;
553                        case 'U': if(     rvstr[0]=='A'){index += 4;}
554                                  else if(rvstr[0]=='G'){index += 5;}
555                                  else{ index = -1000; }
556                                  break;
557                }
558                if(index > 0) {
559                    stacking += Stacking_Energy[index];
560                }
561            }
562            scI.SetStacking(stacking/(float)profNum);
563        }
564        else {
565            scI.SetStacking(1000);
566        }
567    }
568
569    for(int i = 1; i < num; i++) {
570        StemCandidate &sc = pscs->at(i);
571        float stemStacking = 0;
572        int profNum = sc.GetNumSeq();
573        for(int k = 0; k < profNum; k++) {
574            string substr = sc.GetSubstr(k);
575            string rvstr  = sc.GetRvstr(k);
576            for(int j = 0; j < wordLength-1; j++) {
577                int index = 0;
578
579                switch(substr[j]) {
580                    case 'A': if(   rvstr[wordLength-1-j]=='U'){index += 0;}
581                              else{ index = -1000; }
582                              break;
583                    case 'C': if(   rvstr[wordLength-1-j]=='G'){index += 6;}
584                              else{ index = -1000; }
585                              break;
586                    case 'G': if(   rvstr[wordLength-1-j]=='C'){index += 12;}
587                              else if(rvstr[wordLength-1-j]=='U'){index += 18;}
588                              else{ index = -1000; }
589                              break;
590                    case 'U': if(   rvstr[wordLength-1-j]=='A'){index += 24;}
591                              else if(rvstr[wordLength-1-j]=='G'){index += 30;}
592                              else{ index = -1000; }
593                              break;
594                }
595                switch(substr[j+1]){
596                    case 'A': if(   rvstr[wordLength-1-(j+1)]=='U'){index += 0;}
597                              else{ index = -1000; }
598                              break;
599                    case 'C': if(   rvstr[wordLength-1-(j+1)]=='G'){index += 1;}
600                              else{ index = -1000; }
601                              break;
602                    case 'G': if(     rvstr[wordLength-1-(j+1)]=='C'){index += 2;}
603                              else if(rvstr[wordLength-1-(j+1)]=='U'){index += 3;}
604                              else{ index = -1000; }
605                              break;
606                    case 'U': if(     rvstr[wordLength-1-(j+1)]=='A'){index += 4;}
607                              else if(rvstr[wordLength-1-(j+1)]=='G'){index += 5;}
608                              else{ index = -1000; }
609                              break;
610                }
611                if(index > 0) {
612                    stemStacking += Stacking_Energy[index];
613                }
614            }
615            sc.SetStemStacking(stemStacking/(float)profNum);
616        }
617    }
618    return pscs;
619}
620
Note: See TracBrowser for help on using the repository browser.