source: trunk/GDE/MUSCLE/src/anchors.cpp

#include "muscle.h"
#include "msa.h"
#include "objscore.h"

#define TRACE   0

static void WindowSmooth(const SCORE Score[], unsigned uCount, unsigned uWindowLength,
  SCORE SmoothScore[], double dCeil)
        {
#define Ceil(x) ((SCORE) ((x) > dCeil ? dCeil : (x)))

        if (1 != uWindowLength%2)
                Quit("WindowSmooth=%u must be odd", uWindowLength);

        if (uCount <= uWindowLength)
                {
                for (unsigned i = 0; i < uCount; ++i)
                        SmoothScore[i] = 0;
                return;
                }

        const unsigned w2 = uWindowLength/2;
        for (unsigned i = 0; i < w2; ++i)
                {
                SmoothScore[i] = 0;
                SmoothScore[uCount - i - 1] = 0;
                }

        SCORE scoreWindowTotal = 0;
        for (unsigned i = 0; i < uWindowLength; ++i)
                {
                scoreWindowTotal += Ceil(Score[i]);
                }

        for (unsigned i = w2; ; ++i)
                {
                SmoothScore[i] = scoreWindowTotal/uWindowLength;
                if (i == uCount - w2 - 1)
                        break;

                scoreWindowTotal -= Ceil(Score[i - w2]);
                scoreWindowTotal += Ceil(Score[i + w2 + 1]);
                }
#undef Ceil
        }

// Find columns that score above the given threshold.
// A range of scores is defined between the average
// and the maximum. The threshold is a fraction 0.0 .. 1.0
// within that range, where 0.0 is the average score
// and 1.0 is the maximum score.
// "Grade" is by analogy with grading on a curve.
static void FindBestColsGrade(const SCORE Score[], unsigned uCount,
  double dThreshold, unsigned BestCols[], unsigned *ptruBestColCount)
        {
        SCORE scoreTotal = 0;
        for (unsigned uIndex = 0; uIndex < uCount; ++uIndex)
                scoreTotal += Score[uIndex];
        const SCORE scoreAvg = scoreTotal / uCount;

        SCORE scoreMax = MINUS_INFINITY;
        for (unsigned uIndex = 0; uIndex < uCount; ++uIndex)
                if (Score[uIndex] > scoreMax)
                        scoreMax = Score[uIndex];

        unsigned uBestColCount = 0;
        for (unsigned uIndex = 0; uIndex < uCount; ++uIndex)
                {
                const SCORE s = Score[uIndex];
                const double dHeight = (s - scoreAvg)/(scoreMax - scoreAvg);
                if (dHeight >= dThreshold)
                        {
                        BestCols[uBestColCount] = uIndex;
                        ++uBestColCount;
                        }
                }
        *ptruBestColCount = uBestColCount;
        }

// Best col only if all following criteria satisfied:
// (1) Score >= min
// (2) Smoothed score >= min
// (3) No gaps.
static void FindBestColsCombo(const MSA &msa, const SCORE Score[],
  const SCORE SmoothScore[], double dMinScore, double dMinSmoothScore,
  unsigned BestCols[], unsigned *ptruBestColCount)
        {
        const unsigned uColCount = msa.GetColCount();

        unsigned uBestColCount = 0;
        for (unsigned uIndex = 0; uIndex < uColCount; ++uIndex)
                {
                if (Score[uIndex] < dMinScore)
                        continue;
                if (SmoothScore[uIndex] < dMinSmoothScore)
                        continue;
                if (msa.ColumnHasGap(uIndex))
                        continue;
                BestCols[uBestColCount] = uIndex;
                ++uBestColCount;
                }
        *ptruBestColCount = uBestColCount;
        }

static void ListBestCols(const MSA &msa, const SCORE Score[], const SCORE SmoothScore[],
  unsigned BestCols[], unsigned uBestColCount)
        {
        const unsigned uColCount = msa.GetColCount();
        const unsigned uSeqCount = msa.GetSeqCount();

        Log("Col  ");
        for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)
                Log("%u", uSeqIndex%10);
        Log("  ");

        for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)
                {
                Log("%3u  ", uColIndex);
                for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)
                        Log("%c", msa.GetChar(uSeqIndex, uColIndex));

                Log("  %10.3f", Score[uColIndex]);
                Log("  %10.3f", SmoothScore[uColIndex]);

                for (unsigned i = 0; i < uBestColCount; ++i)
                        if (BestCols[i] == uColIndex)
                                Log(" <-- Best");
                Log("\n");
                }
        }

// If two best columns are found within a window, choose
// the highest-scoring. If more than two, choose the one
// closest to the center of the window.
static void MergeBestCols(const SCORE Scores[], const unsigned BestCols[],
  unsigned uBestColCount, unsigned uWindowLength, unsigned AnchorCols[],
  unsigned *ptruAnchorColCount)
        {
        unsigned uAnchorColCount = 0;
        for (unsigned n = 0; n < uBestColCount; /* update inside loop */)
                {
                unsigned uBestColIndex = BestCols[n];
                unsigned uCountWithinWindow = 0;
                for (unsigned i = n + 1; i < uBestColCount; ++i)
                        {
                        unsigned uBestColIndex2 = BestCols[i];
                        if (uBestColIndex2 - uBestColIndex >= uWindowLength)
                                break;
                        ++uCountWithinWindow;
                        }
                unsigned uAnchorCol = uBestColIndex;
                if (1 == uCountWithinWindow)
                        {
                        unsigned uBestColIndex2 = BestCols[n+1];
                        if (Scores[uBestColIndex] > Scores[uBestColIndex2])
                                uAnchorCol = uBestColIndex;
                        else
                                uAnchorCol = uBestColIndex2;
                        }
                else if (uCountWithinWindow > 1)
                        {
                        unsigned uWindowCenter = uBestColIndex + uWindowLength/2;
                        int iClosestDist = uWindowLength;
                        unsigned uClosestCol = uBestColIndex;
                        for (unsigned i = n + 1; i < n + uCountWithinWindow; ++i)
                                {
                                unsigned uColIndex = BestCols[i];
                                int iDist = uColIndex - uBestColIndex;
                                if (iDist < 0)
                                        iDist = -iDist;
                                if (iDist < iClosestDist)
                                        {
                                        uClosestCol = uColIndex;
                                        iClosestDist = iDist;
                                        }
                                }
                        uAnchorCol = uClosestCol;
                        }
                AnchorCols[uAnchorColCount] = uAnchorCol;
                ++uAnchorColCount;
                n += uCountWithinWindow + 1;
                }
        *ptruAnchorColCount = uAnchorColCount;
        }

void FindAnchorCols(const MSA &msa, unsigned AnchorCols[],
  unsigned *ptruAnchorColCount)
        {
        const unsigned uColCount = msa.GetColCount();
        if (uColCount < 16)
                {
                *ptruAnchorColCount = 0;
                return;
                }

        SCORE *MatchScore = new SCORE[uColCount];
        SCORE *SmoothScore = new SCORE[uColCount];
        unsigned *BestCols = new unsigned[uColCount];

        GetLetterScores(msa, MatchScore);
        WindowSmooth(MatchScore, uColCount, g_uSmoothWindowLength, SmoothScore,
          g_dSmoothScoreCeil);

        unsigned uBestColCount;
        FindBestColsCombo(msa, MatchScore, SmoothScore, g_dMinBestColScore, g_dMinSmoothScore,
          BestCols, &uBestColCount);

#if     TRACE
        ListBestCols(msa, MatchScore, SmoothScore, BestCols, uBestColCount);
#endif

        MergeBestCols(MatchScore, BestCols, uBestColCount, g_uAnchorSpacing, AnchorCols,
          ptruAnchorColCount);

        delete[] MatchScore;
        delete[] SmoothScore;
        delete[] BestCols;
        }