--- /dev/null
+#include "Globaldp.hpp"
+
+namespace MXSCARNA {
+
+double Globaldp::ribosum_matrix[16][16]
+= { { -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.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.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 },
+ { -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 },
+ { -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 },
+ { -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 },
+ { -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 },
+ { -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 },
+ { -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 },
+ { -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 },
+ { -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 },
+ { -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 },
+ { -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 },
+ { -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 },
+ { -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 },
+ { -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 }
+};
+
+
+Trimat<float>*
+Globaldp::
+makeProfileBPPMatrix(const MultiSequence *Sequences)
+{
+ int length = Sequences->GetSequence(0)->GetLength();
+ float thr = THR;
+ Trimat<float> *consBppMat = new Trimat<float>(length + 1);
+ fill(consBppMat->begin(), consBppMat->end(), 0);
+
+ int number = Sequences->GetNumSequences();
+ for(int seqNum = 0; seqNum < number; seqNum++) {
+ SafeVector<int> *tmpMap = Sequences->GetSequence(seqNum)->GetMappingNumber();
+ int label = Sequences->GetSequence(seqNum)->GetLabel();
+ BPPMatrix *tmpBppMatrix = BPPMatrices[label];
+
+ for(int i = 1; i <= length ; i++) {
+ int originI = tmpMap->at(i);
+ for(int j = i; j <= length; j++) {
+ int originJ = tmpMap->at(j);
+ if(originI != 0 && originJ != 0) {
+ float tmpProb = tmpBppMatrix->GetProb(originI, originJ);
+
+ if(tmpProb >= thr) {
+ consBppMat->ref(i, j) += tmpProb;
+ }
+ }
+ }
+ }
+ }
+
+ /* compute the mean of base pairing probability */
+ for(int i = 1; i <= length; i++) {
+ for(int j = i; j <= length; j++) {
+ consBppMat->ref(i,j) = consBppMat->ref(i,j)/(float)number;
+ }
+ }
+
+ return consBppMat;
+}
+
+float
+Globaldp::
+incrementalScorePSCPair(const StemCandidate &sc1, const StemCandidate &sc2)
+{
+ int wordLength = WORDLENGTH;
+
+ int pos1, rvpos1;
+ if (sc1.GetDistance() > 0) {
+ pos1 = sc1.GetPosition();
+ rvpos1 = sc1.GetRvposition();
+ }
+ else {
+ pos1 = sc1.GetRvposition();
+ rvpos1 = sc1.GetPosition();
+ }
+
+ int pos2, rvpos2;
+ if (sc2.GetDistance() > 0) {
+ pos2 = sc2.GetPosition();
+ rvpos2 = sc2.GetRvposition();
+ }
+ else {
+ pos2 = sc2.GetRvposition();
+ rvpos2 = sc2.GetPosition();
+ }
+/*
+ cout << "1:" << i << " " << j << " " << sc1.GetDistance() << " " << sc2.GetDistance() << endl;
+ cout << sc1.GetPosition() << " " << sc1.GetRvposition() << " " << sc2.GetPosition() << " " << sc2.GetRvposition() << endl;
+*/
+ float score =
+ posterior[sc1.GetPosition() + wordLength - 1][sc2.GetPosition() + wordLength - 1]
+// * sc1.GetBaseScore(wordLength - 1)
+ * profileBPPMat1->ref(pos1 + wordLength - 1, rvpos1)
+ * posterior[sc1.GetRvposition()][sc2.GetRvposition()]
+// * sc2.GetBaseScore(wordLength - 1);
+ * profileBPPMat2->ref(pos2 + wordLength - 1, rvpos2);
+
+
+/*
+ incrementalWordSimilarity(sc1, sc2, i, j)
+ + incrementalProbabilityScore(sc1, sc2) * multiDeltaScore
+ + incrementalStackingScore(sc1, sc2) * multiDeltaStacking;
+*/
+
+ return score*sc1.GetNumSeq()*sc2.GetNumSeq();
+}
+
+float
+Globaldp::
+incrementalProbabilityScore(const StemCandidate &sc1, const StemCandidate &sc2)
+{
+ int wordLength = WORDLENGTH;
+ return sc1.GetBaseScore(wordLength-1) + sc2.GetBaseScore(wordLength-1);
+}
+
+float
+Globaldp::
+incrementalStackingScore(const StemCandidate &sc1, const StemCandidate &sc2)
+{
+ return - (sc1.GetStacking() + sc2.GetStacking());
+}
+
+float
+Globaldp::
+incrementalWordSimilarity(const StemCandidate &sc1, const StemCandidate &sc2, int i, int j)
+{
+ int numSeq1 = sc1.GetNumSeq();
+ int numSeq2 = sc2.GetNumSeq();
+
+ float score = 0;
+
+ for(int k = 0; k < 16; k++) {
+ for(int l = 0; l < 16; l++) {
+ score += countContConp1[k][i]*countContConp2[l][j]*ribosum_matrix[k][l];
+ }
+ }
+
+ return score/(numSeq1*numSeq2);
+}
+
+float
+Globaldp::
+scorePSCPair(const StemCandidate &sc1, const StemCandidate &sc2)
+{
+
+
+ int wordLength = WORDLENGTH;
+ float score = 0;
+
+ int pos1, rvpos1;
+ if (sc1.GetDistance() > 0) {
+ pos1 = sc1.GetPosition();
+ rvpos1 = sc1.GetRvposition();
+ }
+ else {
+ pos1 = sc1.GetRvposition();
+ rvpos1 = sc1.GetPosition();
+ }
+
+ int pos2, rvpos2;
+ if (sc2.GetDistance() > 0) {
+ pos2 = sc2.GetPosition();
+ rvpos2 = sc2.GetRvposition();
+ }
+ else {
+ pos2 = sc2.GetRvposition();
+ rvpos2 = sc2.GetPosition();
+ }
+
+ for (int k = 0; k < wordLength; k++) {
+ score +=
+ posterior[sc1.GetPosition() + k][sc2.GetPosition() + k]
+ * profileBPPMat1->ref(pos1 + k, rvpos1 + wordLength - 1 - k)
+// * sc1.GetBaseScore(k)
+ * posterior[sc1.GetRvposition() + wordLength - 1 - k][sc2.GetRvposition() + wordLength - 1 - k]
+// * sc2.GetBaseScore(k);
+ * profileBPPMat2->ref(pos2 + k, rvpos2 + wordLength - 1 - k);
+ }
+ // validation code
+ /*
+ if (profileBPPMat1->ref(pos1 , rvpos1 + wordLength - 1) != sc1.GetBaseScore(0)) {
+ cout << "1 " << profileBPPMat1->ref(pos1, rvpos1 + wordLength - 1) << " " << sc1.GetBaseScore(0) << endl;
+ }
+ if (profileBPPMat2->ref(pos2, rvpos2 + wordLength - 1) != sc2.GetBaseScore(0)) {
+ cout << profileBPPMat2->ref(pos2, rvpos2 + wordLength - 1) << " " << sc2.GetBaseScore(0) << endl;
+ }
+ if (profileBPPMat1->ref(pos1 + 1, rvpos1 + wordLength - 1 - 1) != sc1.GetBaseScore(1)) {
+ cout << "1 " << profileBPPMat1->ref(pos1 + 1, rvpos1 + wordLength - 1 - 1) << " " << sc1.GetBaseScore(1) << endl;
+ }
+ if (profileBPPMat2->ref(pos2 + 1, rvpos2 + wordLength - 1 - 1) != sc2.GetBaseScore(1)) {
+ cout << profileBPPMat2->ref(pos2 + 1, rvpos2 + wordLength - 1 - 1) << " " << sc2.GetBaseScore(1) << endl;
+ }*/
+
+// cout << sc1.GetPosition() << " " << sc1.GetRvposition() << " " << sc1.GetDistance() << endl;
+// cout << sc2.GetPosition() << " " << sc2.GetRvposition() << " " << sc2.GetDistance() << endl;
+/*
+ wordSimilarity(sc1, sc2, i, j)
+ + probabilityScore(sc1, sc2) * multiScore
+ + stackingScore(sc1, sc2) * multiStacking
+
+*/
+/*
+ if (sc1.GetDistance() < 0 && sc2.GetDistance() < 0) {
+ cout << "2:" << i << " " << j << " " << sc1.GetDistance() << " " << sc2.GetDistance() << endl;
+ cout << sc1.GetPosition() << " " << sc1.GetRvposition() << " " << sc2.GetPosition() << " " << sc2.GetRvposition() << endl;
+ cout << "2:score" << score << endl;
+
+ }
+*/
+ return score*sc1.GetNumSeq()*sc2.GetNumSeq();
+}
+
+float
+Globaldp::
+wordSimilarity(const StemCandidate &sc1, const StemCandidate &sc2, int i, int j)
+{
+ int wordLength = WORDLENGTH;
+
+ int numSeq1 = sc1.GetNumSeq();
+ int numSeq2 = sc2.GetNumSeq();
+
+ float score = 0;
+
+ for(int k = 0; k < 16; k++) {
+ for(int l = 0; l < 16; l++) {
+ for(int iter = 0; iter < wordLength; iter++) {
+ score += countConp1[iter][k][i]*countConp2[iter][l][j]*ribosum_matrix[k][l];
+ }
+ }
+ }
+
+ return score/(numSeq1*numSeq2);
+}
+
+int
+Globaldp::
+returnBasePairType(char s, char r)
+{
+ if (s == 'A') {
+ if (r == 'A') return 0;
+ else if (r == 'C') return 1;
+ else if (r == 'G') return 2;
+ else if (r == 'U') return 3;
+ }
+ else if (s == 'C') {
+ if (r == 'A') return 4;
+ else if (r == 'C') return 5;
+ else if (r == 'G') return 6;
+ else if (r == 'U') return 7;
+ }
+ else if (s == 'G') {
+ if (r == 'A') return 8;
+ else if (r == 'C') return 9;
+ else if (r == 'G') return 10;
+ else if (r == 'U') return 11;
+ }
+ else if (s == 'U') {
+ if (r == 'A') return 12;
+ else if (r == 'C') return 13;
+ else if (r == 'G') return 14;
+ else if (r == 'U') return 15;
+ }
+
+ return 16;
+}
+
+float
+Globaldp::
+probabilityScore(const StemCandidate &sc1, const StemCandidate &sc2)
+{
+ return sc1.GetScore() + sc2.GetScore();
+}
+
+float
+Globaldp::
+stackingScore(const StemCandidate &sc1, const StemCandidate &sc2)
+{
+ return - (sc1.GetStemStacking() + sc2.GetStemStacking());
+}
+
+float
+Globaldp::
+distancePenalty(const StemCandidate &sc1, const StemCandidate &sc2)
+{
+ return std::sqrt((float)std::abs(sc1.GetDistance() - sc2.GetDistance()));
+}
+
+float
+Globaldp::
+Run()
+{
+ Initialize();
+ DP();
+ float score = traceBack();
+
+ // printMatch();
+ //cout << "score=" << score << endl;
+ return score;
+}
+
+void
+Globaldp::
+Initialize()
+{
+ int nPscs1 = pscs1->size();
+ int nPscs2 = pscs2->size();
+
+
+ for(int i = 0; i < nPscs1; i++) {
+ for(int j = 0; j < nPscs2; j++) {
+ VM[i][j] = 0;
+ VG[i][j] = 0;
+ }
+ }
+
+ VM[0][0] = 0;
+ VG[0][0] = IMPOSSIBLE;
+
+ for (int i = 1; i < nPscs1; i++) {
+ VM[i][0] = IMPOSSIBLE; VG[i][0] = 0;
+ }
+ for (int j = 1; j < nPscs2; j++) {
+ VM[0][j] = IMPOSSIBLE; VG[0][j] = 0;
+ }
+
+ for (int i = 0; i < nPscs1; i++) {
+ for (int j = 0; j < nPscs2; j++) {
+ traceMI[i][j] = 0; traceMJ[i][j] = 0;
+ traceGI[i][j] = 0; traceGJ[i][j] = 0;
+ }
+ }
+
+ int wordLength = WORDLENGTH;
+ int size1 = pscs1->size();
+ int size2 = pscs2->size();
+
+ for(int i = 0; i < wordLength; i++) {
+ for(int j = 0; j < 17; j++) {
+ for(int k = 0; k < (int)pscs1->size(); k++) {
+ countConp1[i][j][k] = 0;
+ }
+ }
+ }
+ for(int i = 0; i < wordLength; i++) {
+ for(int j = 0; j < 17; j++) {
+ for(int k = 0; k < (int)pscs2->size(); k++) {
+ countConp2[i][j][k] = 0;
+ }
+ }
+ }
+ for(int i = 0; i < 17; i++) {
+ for(int j = 0; j < (int)pscs1->size()+1; j++) {
+ countContConp1[i][j] = 0;
+ }
+ }
+ for(int i = 0; i < 17; i++) {
+ for(int j = 0; j < (int)pscs2->size()+1; j++) {
+ countContConp2[i][j] = 0;
+ }
+ }
+
+ for(int iter = 1; iter < size1; iter++) {
+
+ const StemCandidate &sc1 = pscs1->at(iter);
+ int numSeq1 = sc1.GetNumSeq();
+ for (int i = 0; i < wordLength; i++) {
+
+ for(int k = 0; k < numSeq1; k++) {
+// const Sequence *seq = seqs1->GetSequence(k);
+ string substr = sc1.GetSubstr(k);
+ string rvstr = sc1.GetRvstr(k);
+
+ char sChar = substr[i];
+ char rChar = rvstr[wordLength - 1 -i];
+
+ int type = returnBasePairType(sChar, rChar);
+ ++countConp1[i][type][iter];
+ }
+ }
+ for(int k = 0; k < numSeq1; k++) {
+// const Sequence *seq = seqs1->GetSequence(k);
+ string substr = sc1.GetSubstr(k);
+ string rvstr = sc1.GetRvstr(k);
+
+ char sChar = substr[wordLength-1];
+ char rChar = rvstr[0];
+
+ int type = returnBasePairType(sChar, rChar);
+ ++countContConp1[type][iter];
+
+ }
+ }
+ for(int iter = 1; iter < size2; iter++) {
+ const StemCandidate &sc2 = pscs2->at(iter);
+ int numSeq2 = sc2.GetNumSeq();
+ for (int i = 0; i < wordLength; i++) {
+
+ for(int k = 0; k < numSeq2; k++) {
+// const Sequence *seq = seqs2->GetSequence(k);
+ string substr = sc2.GetSubstr(k);
+ string rvstr = sc2.GetRvstr(k);
+
+ char sChar = substr[i];
+ char rChar = rvstr[wordLength - 1 -i];
+
+ int type = returnBasePairType(sChar, rChar);
+ ++countConp2[i][type][iter];
+ }
+ }
+ for(int k = 0; k < numSeq2; k++) {
+// const Sequence *seq = seqs2->GetSequence(k);
+ string substr = sc2.GetSubstr(k);
+ string rvstr = sc2.GetRvstr(k);
+
+ char sChar = substr[wordLength-1];
+ char rChar = rvstr[0];
+
+ int type = returnBasePairType(sChar, rChar);
+ ++countContConp2[type][iter];
+
+ }
+ }
+ profileBPPMat1 = makeProfileBPPMatrix(seqs1);
+ profileBPPMat2 = makeProfileBPPMatrix(seqs2);
+}
+
+void
+Globaldp::
+DP()
+{
+ int nPscs1 = pscs1->size();
+ int nPscs2 = pscs2->size();
+
+ for(int i = 1; i < nPscs1; i++) {
+ const StemCandidate &sc1 = pscs1->at(i);
+
+ for(int j = 1; j < nPscs2; j++) {
+ const StemCandidate &sc2 = pscs2->at(j);
+
+ float max = IMPOSSIBLE;
+ int traceI = 0;
+ int traceJ = 0;
+ int alpha = sc1.GetContPos(), beta = sc2.GetContPos();
+ if( (alpha > 0) && (beta > 0) ) {
+ max = VM[alpha][beta] + incrementalScorePSCPair(sc1, sc2);
+ traceI = alpha;
+ traceJ = beta;
+ }
+
+ float similarity = scorePSCPair(sc1, sc2);
+ int p = sc1.GetBeforePos(), q = sc2.GetBeforePos();
+ float tmpScore = VM[p][q] + similarity;
+// cout << i << " " << j << " " << p << " " << q << " " << tmpScore << " " << VM[p][q] << " " << similarity << " " << endl;
+
+ if(max <= tmpScore) {
+ max = tmpScore;
+ traceI = p;
+ traceJ = q;
+ }
+
+ tmpScore = VG[p][q] + similarity;
+ if(max <= tmpScore) {
+ max = tmpScore;
+ traceI = traceGI[p][q];
+ traceJ = traceGJ[p][q];
+ }
+
+ VM[i][j] = max;
+ traceMI[i][j] = traceI;
+ traceMJ[i][j] = traceJ;
+
+ max = VM[i][j-1];
+ traceI = i;
+ traceJ = j-1;
+
+ tmpScore = VM[i-1][j];
+ if(max <= tmpScore) {
+ max = tmpScore;
+ traceI = i-1;
+ traceJ = j;
+ }
+
+ tmpScore = VG[i-1][j];
+ if(max <= tmpScore) {
+ max = tmpScore;
+ traceI = traceGI[i-1][j];
+ traceJ = traceGJ[i-1][j];
+ }
+
+ tmpScore = VG[i][j-1];
+ if(max <= tmpScore) {
+ max = tmpScore;
+ traceI = traceGI[i][j-1];
+ traceJ = traceGJ[i][j-1];
+ }
+
+ VG[i][j] = max;
+ traceGI[i][j] = traceI;
+ traceGJ[i][j] = traceJ;
+ }
+ }
+}
+
+float
+Globaldp::
+traceBack()
+{
+ int nPscs1 = pscs1->size();
+ int nPscs2 = pscs2->size();
+
+ float max = IMPOSSIBLE;
+ int traceI, traceJ;
+ for (int i = 1; i < nPscs1; i++) {
+ for (int j = 1; j < nPscs2; j++) {
+ if(max < VM[i][j]) {
+ max = VM[i][j];
+ traceI = i;
+ traceJ = j;
+ }
+ }
+ }
+
+ while ( (traceI != 0) && (traceJ != 0) ) {
+ matchPSCS1->push_back(traceI); matchPSCS2->push_back(traceJ);
+ int tmpI = traceMI[traceI][traceJ];
+ int tmpJ = traceMJ[traceI][traceJ];
+ traceI = tmpI; traceJ = tmpJ;
+ }
+
+ if(traceI != 0 && traceJ != 0) {
+ std::cout << "error in trace back of pscs:" << traceI << " " << traceJ << endl;
+ }
+
+ return max;
+}
+
+void
+Globaldp::
+printMatch()
+{
+ int size = matchPSCS1->size();
+ for(int iter = 0; iter < size; iter++) {
+ int i = matchPSCS1->at(iter);
+ int j = matchPSCS2->at(iter);
+
+ const StemCandidate &sc1 = pscs1->at(i);
+ const StemCandidate &sc2 = pscs2->at(j);
+
+ std::cout << iter << "\t" << sc1.GetPosition() << "\t" << sc1.GetRvposition() << "\t" << sc2.GetPosition() << "\t" << sc2.GetRvposition() << endl;
+ }
+
+}
+}
git://source.jalview.org / jabaws.git / blobdiff