4 * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
9 #include "ProfileWithSub.h"
21 ProfileWithSub::ProfileWithSub(int prfLen, int firstS, int lastS)
22 : ProfileBase(prfLen, firstS, lastS)
29 void ProfileWithSub::resetPrf1()
42 void ProfileWithSub::calcProfileWithSub(SeqArray* seqArray, vector<int>* gaps,
43 int matrix[NUMRES][NUMRES], vector<int>* seqWeight)
45 vector<vector<int> > weighting;
46 int sum2, aa, seq, res;
51 int _maxAA = userParameters->getMaxAA();
52 int _gapPos1 = userParameters->getGapPos1();
53 int _gapPos2 = userParameters->getGapPos2();
55 weighting.resize(NUMRES + 2, vector<int>(prfLength + 2));
57 _numSeq = lastSeq - firstSeq;
60 for (seq = firstSeq; seq < lastSeq; seq++)
62 sum2 += (*seqWeight)[seq];
65 for (col = 0; col < prfLength; col++)
67 for (aa = 0; aa <= _maxAA; aa++)
69 weighting[aa][col] = 0;
71 for (seq = firstSeq; seq < lastSeq; seq++)
72 if (aa == (*seqArray)[seq][col])
74 weighting[aa][col] += (*seqWeight)[seq];
77 weighting[_gapPos1][col] = 0;
79 for (seq = firstSeq; seq < lastSeq; seq++)
81 if (_gapPos1 == (*seqArray)[seq][col])
83 weighting[_gapPos1][col] += (*seqWeight)[seq];
87 weighting[_gapPos2][col] = 0;
89 for (seq = firstSeq; seq < lastSeq; seq++)
91 if (_gapPos2 == (*seqArray)[seq][col])
93 weighting[_gapPos2][col] += (*seqWeight)[seq];
98 for (pos = 0; pos < prfLength; pos++)
100 if ((*gaps)[pos] == _numSeq) // If all gaps
102 for (res = 0; res <= _maxAA; res++)
104 profile[pos + 1][res] = matrix[res][_gapPos1];
106 profile[pos + 1][_gapPos1] = matrix[_gapPos1][_gapPos1];
107 profile[pos + 1][_gapPos2] = matrix[_gapPos2][_gapPos1];
111 scale = (float)(_numSeq - (*gaps)[pos]) / (float)_numSeq;
112 for (res = 0; res <= _maxAA; res++)
116 for (aa = 0; aa <= _maxAA; aa++)
118 f += (weighting[aa][pos] * matrix[aa][res]);
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);
127 for (aa = 0; aa <= _maxAA; aa++)
129 f += (weighting[aa][pos] * matrix[aa][_gapPos1]);
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);
137 for (aa = 0; aa <= _maxAA; aa++)
139 f += (weighting[aa][pos] * matrix[aa][_gapPos2]);
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);