binaries/src/clustalw/src/multipleAlign/ProfileStandard.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 "ProfileStandard.h"
  10
  11 namespace clustalw
  12 {
  13
  14 /**
  15  *
  16  * @param prfLen
  17  * @param firstS
  18  * @param lastS
  19  * @return
  20  */
  21 ProfileStandard::ProfileStandard(int prfLen, int firstS, int lastS)
  22  : ProfileBase(prfLen, firstS, lastS)
  23 {
  24
  25 }
  26
  27 /**
  28  *
  29  */
  30 void ProfileStandard::resetPrf2()
  31 {
  32     profile.clear();
  33 }
  34
  35
  36 /**
  37  *
  38  * @param seqArray
  39  * @param seqWeight
  40  */
  41 void ProfileStandard::calcStandardProfile(SeqArray* seqArray, vector<int>* seqWeight)
  42 {
  43     /** DONT FORGET TO CHECK THE SIZES ARE CORRECT */
  44
  45     int sum1, sum2;
  46     int i, d;
  47     int r;
  48     int _maxAA = userParameters->getMaxAA();
  49     profile.resize(prfLength + 2, vector<int>(LENCOL + 2));
  50     int _gapPos1 = userParameters->getGapPos1();
  51     int _gapPos2 = userParameters->getGapPos2();
  52
  53     for (r = 0; r < prfLength; r++)
  54     {
  55         /*
  56          * calculate sum2 = number of residues found in this column
  57          */
  58         sum2 = 0;
  59         for (i = firstSeq; i < lastSeq; i++)
  60         {
  61             sum2 += (*seqWeight)[i];
  62         }
  63         /*
  64          * only include matrix comparison scores for those residue types found in this
  65          * column
  66          */
  67         if (sum2 == 0)
  68         {
  69             for (d = 0; d <= _maxAA; d++)
  70             {
  71                 profile[r + 1][d] = 0;
  72             }
  73
  74             profile[r + 1][_gapPos1] = 0;
  75             profile[r + 1][_gapPos2] = 0;
  76         }
  77         else
  78         {
  79             for (d = 0; d <= _maxAA; d++)
  80             {
  81                 sum1 = 0;
  82                 for (i = firstSeq; i < lastSeq; i++)
  83                 {
  84                     if (d == (*seqArray)[i][r])
  85                     {
  86                         sum1 += (*seqWeight)[i];
  87                     }
  88                 }
  89                 profile[r + 1][d] = (int)(10 *(float)sum1 / (float)sum2);
  90             }
  91             sum1 = 0;
  92
  93             for (i = firstSeq; i < lastSeq; i++)
  94             {
  95                 if (_gapPos1 == (*seqArray)[i][r])
  96                 {
  97                     sum1 += (*seqWeight)[i];
  98                 }
  99             }
 100             profile[r + 1][_gapPos1] = (int)(10 *(float)sum1 / (float)sum2);
 101             sum1 = 0;
 102
 103             for (i = firstSeq; i < lastSeq; i++)
 104             {
 105                 if (_gapPos2 == (*seqArray)[i][r])
 106                 {
 107                     sum1 += (*seqWeight)[i];
 108                 }
 109             }
 110             profile[r + 1][_gapPos2] = (int)(10 *(float)sum1 / (float)sum2);
 111         }
 112     }
 113
 114 }
 115
 116
 117 }