binaries/src/clustalw/src/multipleAlign/ProfileWithSub.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 "ProfileWithSub.h"
  10
  11 namespace clustalw
  12 {
  13
  14 /**
  15  *
  16  * @param prfLen
  17  * @param firstS
  18  * @param lastS
  19  * @return
  20  */
  21 ProfileWithSub::ProfileWithSub(int prfLen, int firstS, int lastS)
  22  : ProfileBase(prfLen, firstS, lastS)
  23 {
  24 }
  25
  26 /**
  27  *
  28  */
  29 void ProfileWithSub::resetPrf1()
  30 {
  31     profile.clear();
  32 }
  33
  34
  35 /**
  36  *
  37  * @param seqArray
  38  * @param gaps
  39  * @param matrix[][]
  40  * @param seqWeight
  41  */
  42 void ProfileWithSub::calcProfileWithSub(SeqArray* seqArray, vector<int>* gaps,
  43                               int matrix[NUMRES][NUMRES], vector<int>* seqWeight)
  44 {
  45     vector<vector<int> > weighting;
  46     int sum2, aa, seq, res;
  47     int _numSeq;
  48     int col, pos;
  49     int f;
  50     float scale;
  51     int _maxAA = userParameters->getMaxAA();
  52     int _gapPos1 = userParameters->getGapPos1();
  53     int _gapPos2 = userParameters->getGapPos2();
  54
  55     weighting.resize(NUMRES + 2, vector<int>(prfLength + 2));
  56
  57     _numSeq = lastSeq - firstSeq;
  58
  59     sum2 = 0;
  60     for (seq = firstSeq; seq < lastSeq; seq++)
  61     {
  62         sum2 += (*seqWeight)[seq];
  63     }
  64
  65     for (col = 0; col < prfLength; col++)
  66     {
  67         for (aa = 0; aa <= _maxAA; aa++)
  68         {
  69             weighting[aa][col] = 0;
  70
  71             for (seq = firstSeq; seq < lastSeq; seq++)
  72                 if (aa == (*seqArray)[seq][col])
  73                 {
  74                     weighting[aa][col] += (*seqWeight)[seq];
  75                 }
  76         }
  77         weighting[_gapPos1][col] = 0;
  78
  79         for (seq = firstSeq; seq < lastSeq; seq++)
  80         {
  81             if (_gapPos1 == (*seqArray)[seq][col])
  82             {
  83                 weighting[_gapPos1][col] += (*seqWeight)[seq];
  84             }
  85         }
  86
  87         weighting[_gapPos2][col] = 0;
  88
  89         for (seq = firstSeq; seq < lastSeq; seq++)
  90         {
  91             if (_gapPos2 == (*seqArray)[seq][col])
  92             {
  93                 weighting[_gapPos2][col] += (*seqWeight)[seq];
  94             }
  95         }
  96     }
  97
  98     for (pos = 0; pos < prfLength; pos++)
  99     {
 100         if ((*gaps)[pos] == _numSeq) // If all gaps
 101         {
 102             for (res = 0; res <= _maxAA; res++)
 103             {
 104                 profile[pos + 1][res] = matrix[res][_gapPos1];
 105             }
 106             profile[pos + 1][_gapPos1] = matrix[_gapPos1][_gapPos1];
 107             profile[pos + 1][_gapPos2] = matrix[_gapPos2][_gapPos1];
 108         }
 109         else
 110         {
 111             scale = (float)(_numSeq - (*gaps)[pos]) / (float)_numSeq;
 112             for (res = 0; res <= _maxAA; res++)
 113             {
 114                 f = 0;
 115
 116                 for (aa = 0; aa <= _maxAA; aa++)
 117                 {
 118                     f += (weighting[aa][pos] * matrix[aa][res]);
 119                 }
 120
 121                 f += (weighting[_gapPos1][pos] * matrix[_gapPos1][res]);
 122                 f += (weighting[_gapPos2][pos] * matrix[_gapPos2][res]);
 123                 profile[pos + 1][res] = (int)(((float)f / (float)sum2) * scale);
 124             }
 125             f = 0;
 126
 127             for (aa = 0; aa <= _maxAA; aa++)
 128             {
 129                 f += (weighting[aa][pos] * matrix[aa][_gapPos1]);
 130             }
 131
 132             f += (weighting[_gapPos1][pos] * matrix[_gapPos1][_gapPos1]);
 133             f += (weighting[_gapPos2][pos] * matrix[_gapPos2][_gapPos1]);
 134             profile[pos + 1][_gapPos1] = (int)(((float)f / (float)sum2) * scale);
 135             f = 0;
 136
 137             for (aa = 0; aa <= _maxAA; aa++)
 138             {
 139                 f += (weighting[aa][pos] * matrix[aa][_gapPos2]);
 140             }
 141
 142             f += (weighting[_gapPos1][pos] * matrix[_gapPos1][_gapPos2]);
 143             f += (weighting[_gapPos2][pos] * matrix[_gapPos2][_gapPos2]);
 144             profile[pos + 1][_gapPos2] = (int)(((float)f / (float)sum2) * scale);
 145         }
 146     }
 147 }
 148
 149
 150 }