binaries/src/clustalw/src/alignment/ObjectiveScore.cpp

   1 /**
   2  * Author: Mark Larkin
   3  *
   4  * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
   5  */
   6 #ifdef HAVE_CONFIG_H
   7     #include "config.h"
   8 #endif
   9 #include <sstream>
  10 #include "ObjectiveScore.h"
  11 #include "Alignment.h"
  12 #include "../general/userparams.h"
  13 #include "../general/debuglogObject.h"
  14
  15 namespace clustalw
  16 {
  17
  18 ObjectiveScore::ObjectiveScore()
  19  : score(0),
  20    alignToScore(0),
  21    scale(100000),
  22    sagaGapEx(12),
  23    sagaGapOp(8),
  24    gapPos1(userParameters->getGapPos1()),
  25    gapPos2(userParameters->getGapPos2())
  26 {
  27
  28 }
  29
  30 long ObjectiveScore::getScore(const Alignment* alnToScore)
  31 {
  32     #if DEBUGFULL
  33         if(logObject && DEBUGLOG)
  34         {
  35             logObject->logMsg("In getScore function");
  36         }
  37     #endif
  38
  39     alignToScore = alnToScore;
  40     // If it doesnt point to any object return 0;
  41     if(!alignToScore)
  42     {
  43         return 0;
  44     }
  45
  46     int maxRes = subMatrix->getAlnScoreMatrix(matrix);
  47     if (maxRes == 0)
  48     {
  49         utilityObject->error("Matrix for alignment scoring not found\n");
  50         return 0;
  51     }
  52
  53     const vector<int>* seqWeight = alignToScore->getSeqWeights();
  54     vector<float> normalisedSeqWeights;
  55
  56     calcNormalisedSeqWeights(seqWeight, &normalisedSeqWeights);
  57
  58     int seq1, seq2;
  59
  60     float w1, w2 = 1.0;
  61     float weight = 1.0;
  62     score = 0;
  63     float pwScoreLetters = 0, pwScoreGaps = 0, pwScore = 0;
  64     float scoreTotal = 0.0;
  65     int numSeqs = alignToScore->getNumSeqs();
  66     int sizeNormalSeqWeight = normalisedSeqWeights.size();
  67
  68     for (seq1 = 1; seq1 <= numSeqs && seq1 <= sizeNormalSeqWeight; seq1++)
  69     {
  70         w1 = normalisedSeqWeights[seq1 - 1];
  71
  72         for (seq2 = seq1 + 1; seq2 <= numSeqs && seq2 <= sizeNormalSeqWeight; seq2++)
  73         {
  74             w2 = normalisedSeqWeights[seq2 - 1];
  75             weight = w1 * w2;
  76
  77             pwScoreLetters = scoreLetters(seq1, seq2);
  78             pwScoreGaps = scoreGaps(seq1, seq2);
  79             pwScore = pwScoreLetters + pwScoreGaps;
  80
  81             scoreTotal += weight * pwScore;
  82             #if DEBUGFULL
  83                 if(logObject && DEBUGLOG)
  84                 {
  85                     ostringstream outs;
  86                     outs << " weight = " << weight
  87                          << " scoreLetters = " << pwScoreLetters << " scoreGaps = "
  88                          << pwScoreGaps << " scorePair = "
  89                          << pwScore << "\n"
  90                          << "scoreTotal = " << scoreTotal << "\n";
  91                     logObject->logMsg(outs.str());
  92                 }
  93             #endif
  94         }
  95     }
  96
  97     score = static_cast<long>(scoreTotal);
  98     #if DEBUGFULL
  99         if(logObject && DEBUGLOG)
 100         {
 101             ostringstream outs;
 102             outs << " score = " << score;
 103             logObject->logMsg(outs.str());
 104         }
 105     #endif
 106
 107     utilityObject->info("Alignment Score %d\n", score);
 108     return score;
 109 }
 110
 111 float ObjectiveScore::scoreLetters(int seq1, int seq2)
 112 {
 113     if(!alignToScore)
 114     {
 115         return 0;
 116     }
 117     const unsigned numColsSeq1 = alignToScore->getSeqLength(seq1);
 118     const unsigned numColsSeq2 = alignToScore->getSeqLength(seq2);
 119
 120     if(numColsSeq1 != numColsSeq2)
 121     {
 122         return 0; // The sequences should be the same length after alignment.
 123     }
 124
 125     float scoreLetters = 0;
 126     unsigned colStart = 1;
 127     bool gap1, gap2;
 128
 129     for(unsigned col = 1; col < numColsSeq1; col++)
 130     {
 131         gap1 = alignToScore->isGap(seq1, col);
 132         gap2 = alignToScore->isGap(seq2, col);
 133
 134         if(!gap1 || !gap2)
 135         {
 136             colStart = col;
 137             break;
 138         }
 139     }
 140
 141     unsigned colEnd = numColsSeq1;
 142
 143     for(unsigned col = numColsSeq1; col >= 1; --col)
 144     {
 145         gap1 = alignToScore->isGap(seq1, col);
 146         gap2 = alignToScore->isGap(seq2, col);
 147
 148         if(!gap1 || !gap2)
 149         {
 150             colEnd = col;
 151             break;
 152         }
 153     }
 154
 155     const SeqArray* seqArray = alignToScore->getSeqArray();
 156     int scoreMatch = 0;
 157
 158
 159     for(unsigned col = colStart; col <= colEnd; col++)
 160     {
 161         int res1 = (*seqArray)[seq1][col];
 162         int res2 = (*seqArray)[seq2][col];
 163         scoreMatch = matrix[res1][res2];
 164         scoreLetters += scoreMatch;
 165     }
 166
 167     return scoreLetters;
 168 }
 169
 170 float ObjectiveScore::scoreGaps(int seq1, int seq2)
 171 {
 172     if(!alignToScore)
 173     {
 174         return 0;
 175     }
 176     const unsigned numColsSeq1 = alignToScore->getSeqLength(seq1);
 177     const unsigned numColsSeq2 = alignToScore->getSeqLength(seq2);
 178
 179     if(numColsSeq1 != numColsSeq2)
 180     {
 181         return 0; // The sequences should be the same length after alignment.
 182     }
 183
 184     unsigned colStart = 1;
 185     bool gap1, gap2;
 186
 187     for(unsigned col = 1; col < numColsSeq1; col++)
 188     {
 189         gap1 = alignToScore->isGap(seq1, col);
 190         gap2 = alignToScore->isGap(seq2, col);
 191
 192         if(!gap1 || !gap2)
 193         {
 194             colStart = col;
 195             break;
 196         }
 197     }
 198
 199     unsigned colEnd = numColsSeq1;
 200
 201     for(unsigned col = numColsSeq1; col >= 1; --col)
 202     {
 203         gap1 = alignToScore->isGap(seq1, col);
 204         gap2 = alignToScore->isGap(seq2, col);
 205
 206         if(!gap1 || !gap2)
 207         {
 208             colEnd = col;
 209             break;
 210         }
 211     }
 212
 213     bool inGap1 = false;
 214     bool inGap2 = false;
 215     float gapOpen = userParameters->getGapOpen();
 216     float gapExtend = userParameters->getGapExtend();
 217
 218     float scoreGaps = 0;
 219     for(unsigned col = colStart; col <= colEnd; col++)
 220     {
 221         gap1 = alignToScore->isGap(seq1, col);
 222         gap2 = alignToScore->isGap(seq2, col);
 223
 224         if(gap1 && gap2)
 225         {
 226             continue;
 227         }
 228         if(gap1)
 229         {
 230             if(!inGap1)
 231             {
 232                 // NOTE I left out the option of having different end gaps stuff.
 233                 scoreGaps += gapOpen;
 234                 inGap1 = true; // Opening a gap in seq1
 235             }
 236             else
 237             {
 238                 // Already in a gap
 239                 scoreGaps += gapExtend;
 240             }
 241             continue;
 242         }
 243         else if(gap2)
 244         {
 245             if(!inGap2)
 246             {
 247                 // NOTE I left out the option of having different end gaps stuff.
 248                 scoreGaps += gapOpen;
 249                 inGap2 = true; // Opening a gap in seq2
 250             }
 251             else
 252             {
 253                 // Already in a gap
 254                 scoreGaps += gapExtend;
 255             }
 256             continue;
 257         }
 258         inGap1 = inGap2 = false;
 259     }
 260     return scoreGaps;
 261 }
 262
 263 void ObjectiveScore::calcNormalisedSeqWeights(const vector<int>* seqWeight,
 264                                               vector<float>* normSeqWeight)
 265 {
 266     if(!seqWeight || !normSeqWeight)
 267     {
 268         return;
 269     }
 270
 271     int sumWeights = 0;
 272     for(int i = 0; i < (int)seqWeight->size() - 1; i++)
 273     {
 274         sumWeights += (*seqWeight)[i];
 275     }
 276
 277     normSeqWeight->resize(seqWeight->size());
 278     for(int i = 0; i < (int)seqWeight->size() - 1; i++)
 279     {
 280         (*normSeqWeight)[i] = static_cast<float>((*seqWeight)[i]) /
 281                               static_cast<float>(sumWeights);
 282     #if DEBUGFULL
 283         if(logObject && DEBUGLOG)
 284         {
 285             ostringstream outs;
 286             outs << " normSeqWeight[i] = " << (*normSeqWeight)[i];
 287             logObject->logMsg(outs.str());
 288         }
 289     #endif
 290     }
 291
 292 }
 293
 294 }
 295