4 * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
10 #include "ObjectiveScore.h"
11 #include "Alignment.h"
12 #include "../general/userparams.h"
13 #include "../general/debuglogObject.h"
18 ObjectiveScore::ObjectiveScore()
24 gapPos1(userParameters->getGapPos1()),
25 gapPos2(userParameters->getGapPos2())
30 long ObjectiveScore::getScore(const Alignment* alnToScore)
33 if(logObject && DEBUGLOG)
35 logObject->logMsg("In getScore function");
39 alignToScore = alnToScore;
40 // If it doesnt point to any object return 0;
46 int maxRes = subMatrix->getAlnScoreMatrix(matrix);
49 utilityObject->error("Matrix for alignment scoring not found\n");
53 const vector<int>* seqWeight = alignToScore->getSeqWeights();
54 vector<float> normalisedSeqWeights;
56 calcNormalisedSeqWeights(seqWeight, &normalisedSeqWeights);
63 float pwScoreLetters = 0, pwScoreGaps = 0, pwScore = 0;
64 float scoreTotal = 0.0;
65 int numSeqs = alignToScore->getNumSeqs();
66 int sizeNormalSeqWeight = normalisedSeqWeights.size();
68 for (seq1 = 1; seq1 <= numSeqs && seq1 <= sizeNormalSeqWeight; seq1++)
70 w1 = normalisedSeqWeights[seq1 - 1];
72 for (seq2 = seq1 + 1; seq2 <= numSeqs && seq2 <= sizeNormalSeqWeight; seq2++)
74 w2 = normalisedSeqWeights[seq2 - 1];
77 pwScoreLetters = scoreLetters(seq1, seq2);
78 pwScoreGaps = scoreGaps(seq1, seq2);
79 pwScore = pwScoreLetters + pwScoreGaps;
81 scoreTotal += weight * pwScore;
83 if(logObject && DEBUGLOG)
86 outs << " weight = " << weight
87 << " scoreLetters = " << pwScoreLetters << " scoreGaps = "
88 << pwScoreGaps << " scorePair = "
90 << "scoreTotal = " << scoreTotal << "\n";
91 logObject->logMsg(outs.str());
97 score = static_cast<long>(scoreTotal);
99 if(logObject && DEBUGLOG)
102 outs << " score = " << score;
103 logObject->logMsg(outs.str());
107 utilityObject->info("Alignment Score %d\n", score);
111 float ObjectiveScore::scoreLetters(int seq1, int seq2)
117 const unsigned numColsSeq1 = alignToScore->getSeqLength(seq1);
118 const unsigned numColsSeq2 = alignToScore->getSeqLength(seq2);
120 if(numColsSeq1 != numColsSeq2)
122 return 0; // The sequences should be the same length after alignment.
125 float scoreLetters = 0;
126 unsigned colStart = 1;
129 for(unsigned col = 1; col < numColsSeq1; col++)
131 gap1 = alignToScore->isGap(seq1, col);
132 gap2 = alignToScore->isGap(seq2, col);
141 unsigned colEnd = numColsSeq1;
143 for(unsigned col = numColsSeq1; col >= 1; --col)
145 gap1 = alignToScore->isGap(seq1, col);
146 gap2 = alignToScore->isGap(seq2, col);
155 const SeqArray* seqArray = alignToScore->getSeqArray();
159 for(unsigned col = colStart; col <= colEnd; col++)
161 int res1 = (*seqArray)[seq1][col];
162 int res2 = (*seqArray)[seq2][col];
163 scoreMatch = matrix[res1][res2];
164 scoreLetters += scoreMatch;
170 float ObjectiveScore::scoreGaps(int seq1, int seq2)
176 const unsigned numColsSeq1 = alignToScore->getSeqLength(seq1);
177 const unsigned numColsSeq2 = alignToScore->getSeqLength(seq2);
179 if(numColsSeq1 != numColsSeq2)
181 return 0; // The sequences should be the same length after alignment.
184 unsigned colStart = 1;
187 for(unsigned col = 1; col < numColsSeq1; col++)
189 gap1 = alignToScore->isGap(seq1, col);
190 gap2 = alignToScore->isGap(seq2, col);
199 unsigned colEnd = numColsSeq1;
201 for(unsigned col = numColsSeq1; col >= 1; --col)
203 gap1 = alignToScore->isGap(seq1, col);
204 gap2 = alignToScore->isGap(seq2, col);
215 float gapOpen = userParameters->getGapOpen();
216 float gapExtend = userParameters->getGapExtend();
219 for(unsigned col = colStart; col <= colEnd; col++)
221 gap1 = alignToScore->isGap(seq1, col);
222 gap2 = alignToScore->isGap(seq2, col);
232 // NOTE I left out the option of having different end gaps stuff.
233 scoreGaps += gapOpen;
234 inGap1 = true; // Opening a gap in seq1
239 scoreGaps += gapExtend;
247 // NOTE I left out the option of having different end gaps stuff.
248 scoreGaps += gapOpen;
249 inGap2 = true; // Opening a gap in seq2
254 scoreGaps += gapExtend;
258 inGap1 = inGap2 = false;
263 void ObjectiveScore::calcNormalisedSeqWeights(const vector<int>* seqWeight,
264 vector<float>* normSeqWeight)
266 if(!seqWeight || !normSeqWeight)
272 for(int i = 0; i < (int)seqWeight->size() - 1; i++)
274 sumWeights += (*seqWeight)[i];
277 normSeqWeight->resize(seqWeight->size());
278 for(int i = 0; i < (int)seqWeight->size() - 1; i++)
280 (*normSeqWeight)[i] = static_cast<float>((*seqWeight)[i]) /
281 static_cast<float>(sumWeights);
283 if(logObject && DEBUGLOG)
286 outs << " normSeqWeight[i] = " << (*normSeqWeight)[i];
287 logObject->logMsg(outs.str());