+++ /dev/null
-/**
- * Author: Mark Larkin
- *
- * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
- */
-#ifdef HAVE_CONFIG_H
- #include "config.h"
-#endif
-#include <sstream>
-#include "ObjectiveScore.h"
-#include "Alignment.h"
-#include "../general/userparams.h"
-#include "../general/debuglogObject.h"
-
-namespace clustalw
-{
-
-ObjectiveScore::ObjectiveScore()
- : score(0),
- alignToScore(0),
- scale(100000),
- sagaGapEx(12),
- sagaGapOp(8),
- gapPos1(userParameters->getGapPos1()),
- gapPos2(userParameters->getGapPos2())
-{
-
-}
-
-long ObjectiveScore::getScore(const Alignment* alnToScore)
-{
- #if DEBUGFULL
- if(logObject && DEBUGLOG)
- {
- logObject->logMsg("In getScore function");
- }
- #endif
-
- alignToScore = alnToScore;
- // If it doesnt point to any object return 0;
- if(!alignToScore)
- {
- return 0;
- }
-
- int maxRes = subMatrix->getAlnScoreMatrix(matrix);
- if (maxRes == 0)
- {
- utilityObject->error("Matrix for alignment scoring not found\n");
- return 0;
- }
-
- const vector<int>* seqWeight = alignToScore->getSeqWeights();
- vector<float> normalisedSeqWeights;
-
- calcNormalisedSeqWeights(seqWeight, &normalisedSeqWeights);
-
- int seq1, seq2;
-
- float w1, w2 = 1.0;
- float weight = 1.0;
- score = 0;
- float pwScoreLetters = 0, pwScoreGaps = 0, pwScore = 0;
- float scoreTotal = 0.0;
- int numSeqs = alignToScore->getNumSeqs();
- int sizeNormalSeqWeight = normalisedSeqWeights.size();
-
- for (seq1 = 1; seq1 <= numSeqs && seq1 <= sizeNormalSeqWeight; seq1++)
- {
- w1 = normalisedSeqWeights[seq1 - 1];
-
- for (seq2 = seq1 + 1; seq2 <= numSeqs && seq2 <= sizeNormalSeqWeight; seq2++)
- {
- w2 = normalisedSeqWeights[seq2 - 1];
- weight = w1 * w2;
-
- pwScoreLetters = scoreLetters(seq1, seq2);
- pwScoreGaps = scoreGaps(seq1, seq2);
- pwScore = pwScoreLetters + pwScoreGaps;
-
- scoreTotal += weight * pwScore;
- #if DEBUGFULL
- if(logObject && DEBUGLOG)
- {
- ostringstream outs;
- outs << " weight = " << weight
- << " scoreLetters = " << pwScoreLetters << " scoreGaps = "
- << pwScoreGaps << " scorePair = "
- << pwScore << "\n"
- << "scoreTotal = " << scoreTotal << "\n";
- logObject->logMsg(outs.str());
- }
- #endif
- }
- }
-
- score = static_cast<long>(scoreTotal);
- #if DEBUGFULL
- if(logObject && DEBUGLOG)
- {
- ostringstream outs;
- outs << " score = " << score;
- logObject->logMsg(outs.str());
- }
- #endif
-
- utilityObject->info("Alignment Score %d\n", score);
- return score;
-}
-
-float ObjectiveScore::scoreLetters(int seq1, int seq2)
-{
- if(!alignToScore)
- {
- return 0;
- }
- const unsigned numColsSeq1 = alignToScore->getSeqLength(seq1);
- const unsigned numColsSeq2 = alignToScore->getSeqLength(seq2);
-
- if(numColsSeq1 != numColsSeq2)
- {
- return 0; // The sequences should be the same length after alignment.
- }
-
- float scoreLetters = 0;
- unsigned colStart = 1;
- bool gap1, gap2;
-
- for(unsigned col = 1; col < numColsSeq1; col++)
- {
- gap1 = alignToScore->isGap(seq1, col);
- gap2 = alignToScore->isGap(seq2, col);
-
- if(!gap1 || !gap2)
- {
- colStart = col;
- break;
- }
- }
-
- unsigned colEnd = numColsSeq1;
-
- for(unsigned col = numColsSeq1; col >= 1; --col)
- {
- gap1 = alignToScore->isGap(seq1, col);
- gap2 = alignToScore->isGap(seq2, col);
-
- if(!gap1 || !gap2)
- {
- colEnd = col;
- break;
- }
- }
-
- const SeqArray* seqArray = alignToScore->getSeqArray();
- int scoreMatch = 0;
-
-
- for(unsigned col = colStart; col <= colEnd; col++)
- {
- int res1 = (*seqArray)[seq1][col];
- int res2 = (*seqArray)[seq2][col];
- scoreMatch = matrix[res1][res2];
- scoreLetters += scoreMatch;
- }
-
- return scoreLetters;
-}
-
-float ObjectiveScore::scoreGaps(int seq1, int seq2)
-{
- if(!alignToScore)
- {
- return 0;
- }
- const unsigned numColsSeq1 = alignToScore->getSeqLength(seq1);
- const unsigned numColsSeq2 = alignToScore->getSeqLength(seq2);
-
- if(numColsSeq1 != numColsSeq2)
- {
- return 0; // The sequences should be the same length after alignment.
- }
-
- unsigned colStart = 1;
- bool gap1, gap2;
-
- for(unsigned col = 1; col < numColsSeq1; col++)
- {
- gap1 = alignToScore->isGap(seq1, col);
- gap2 = alignToScore->isGap(seq2, col);
-
- if(!gap1 || !gap2)
- {
- colStart = col;
- break;
- }
- }
-
- unsigned colEnd = numColsSeq1;
-
- for(unsigned col = numColsSeq1; col >= 1; --col)
- {
- gap1 = alignToScore->isGap(seq1, col);
- gap2 = alignToScore->isGap(seq2, col);
-
- if(!gap1 || !gap2)
- {
- colEnd = col;
- break;
- }
- }
-
- bool inGap1 = false;
- bool inGap2 = false;
- float gapOpen = userParameters->getGapOpen();
- float gapExtend = userParameters->getGapExtend();
-
- float scoreGaps = 0;
- for(unsigned col = colStart; col <= colEnd; col++)
- {
- gap1 = alignToScore->isGap(seq1, col);
- gap2 = alignToScore->isGap(seq2, col);
-
- if(gap1 && gap2)
- {
- continue;
- }
- if(gap1)
- {
- if(!inGap1)
- {
- // NOTE I left out the option of having different end gaps stuff.
- scoreGaps += gapOpen;
- inGap1 = true; // Opening a gap in seq1
- }
- else
- {
- // Already in a gap
- scoreGaps += gapExtend;
- }
- continue;
- }
- else if(gap2)
- {
- if(!inGap2)
- {
- // NOTE I left out the option of having different end gaps stuff.
- scoreGaps += gapOpen;
- inGap2 = true; // Opening a gap in seq2
- }
- else
- {
- // Already in a gap
- scoreGaps += gapExtend;
- }
- continue;
- }
- inGap1 = inGap2 = false;
- }
- return scoreGaps;
-}
-
-void ObjectiveScore::calcNormalisedSeqWeights(const vector<int>* seqWeight,
- vector<float>* normSeqWeight)
-{
- if(!seqWeight || !normSeqWeight)
- {
- return;
- }
-
- int sumWeights = 0;
- for(int i = 0; i < (int)seqWeight->size() - 1; i++)
- {
- sumWeights += (*seqWeight)[i];
- }
-
- normSeqWeight->resize(seqWeight->size());
- for(int i = 0; i < (int)seqWeight->size() - 1; i++)
- {
- (*normSeqWeight)[i] = static_cast<float>((*seqWeight)[i]) /
- static_cast<float>(sumWeights);
- #if DEBUGFULL
- if(logObject && DEBUGLOG)
- {
- ostringstream outs;
- outs << " normSeqWeight[i] = " << (*normSeqWeight)[i];
- logObject->logMsg(outs.str());
- }
- #endif
- }
-
-}
-
-}
-
git://source.jalview.org / jabaws.git / blobdiff