binaries/src/muscle/color.cpp

   1 #include "muscle.h"\r
   2 #include "msa.h"\r
   3 \r
   4 static int Blosum62[23][23] =\r
   5         {\r
   6 //   A   B   C   D   E    F   G   H   I   K    L   M   N   P   Q    R   S   T   V   W    X   Y   Z \r
   7         +4, -2, +0, -2, -1,  -2, +0, -2, -1, -1,  -1, -1, -2, -1, -1,  -1, +1, +0, +0, -3,  -1, -2, -1,  // A\r
   8         -2, +6, -3, +6, +2,  -3, -1, -1, -3, -1,  -4, -3, +1, -1, +0,  -2, +0, -1, -3, -4,  -1, -3, +2,  // B\r
   9         +0, -3, +9, -3, -4,  -2, -3, -3, -1, -3,  -1, -1, -3, -3, -3,  -3, -1, -1, -1, -2,  -1, -2, -4,  // C\r
  10         -2, +6, -3, +6, +2,  -3, -1, -1, -3, -1,  -4, -3, +1, -1, +0,  -2, +0, -1, -3, -4,  -1, -3, +2,  // D\r
  11         -1, +2, -4, +2, +5,  -3, -2, +0, -3, +1,  -3, -2, +0, -1, +2,  +0, +0, -1, -2, -3,  -1, -2, +5,  // E\r
  12         \r
  13         -2, -3, -2, -3, -3,  +6, -3, -1, +0, -3,  +0, +0, -3, -4, -3,  -3, -2, -2, -1, +1,  -1, +3, -3,  // F\r
  14         +0, -1, -3, -1, -2,  -3, +6, -2, -4, -2,  -4, -3, +0, -2, -2,  -2, +0, -2, -3, -2,  -1, -3, -2,  // G\r
  15         -2, -1, -3, -1, +0,  -1, -2, +8, -3, -1,  -3, -2, +1, -2, +0,  +0, -1, -2, -3, -2,  -1, +2, +0,  // H\r
  16         -1, -3, -1, -3, -3,  +0, -4, -3, +4, -3,  +2, +1, -3, -3, -3,  -3, -2, -1, +3, -3,  -1, -1, -3,  // I\r
  17         -1, -1, -3, -1, +1,  -3, -2, -1, -3, +5,  -2, -1, +0, -1, +1,  +2, +0, -1, -2, -3,  -1, -2, +1,  // K\r
  18         \r
  19         -1, -4, -1, -4, -3,  +0, -4, -3, +2, -2,  +4, +2, -3, -3, -2,  -2, -2, -1, +1, -2,  -1, -1, -3,  // L\r
  20         -1, -3, -1, -3, -2,  +0, -3, -2, +1, -1,  +2, +5, -2, -2, +0,  -1, -1, -1, +1, -1,  -1, -1, -2,  // M\r
  21         -2, +1, -3, +1, +0,  -3, +0, +1, -3, +0,  -3, -2, +6, -2, +0,  +0, +1, +0, -3, -4,  -1, -2, +0,  // N\r
  22         -1, -1, -3, -1, -1,  -4, -2, -2, -3, -1,  -3, -2, -2, +7, -1,  -2, -1, -1, -2, -4,  -1, -3, -1,  // P\r
  23         -1, +0, -3, +0, +2,  -3, -2, +0, -3, +1,  -2, +0, +0, -1, +5,  +1, +0, -1, -2, -2,  -1, -1, +2,  // Q\r
  24         \r
  25         -1, -2, -3, -2, +0,  -3, -2, +0, -3, +2,  -2, -1, +0, -2, +1,  +5, -1, -1, -3, -3,  -1, -2, +0,  // R\r
  26         +1, +0, -1, +0, +0,  -2, +0, -1, -2, +0,  -2, -1, +1, -1, +0,  -1, +4, +1, -2, -3,  -1, -2, +0,  // S\r
  27         +0, -1, -1, -1, -1,  -2, -2, -2, -1, -1,  -1, -1, +0, -1, -1,  -1, +1, +5, +0, -2,  -1, -2, -1,  // T\r
  28         +0, -3, -1, -3, -2,  -1, -3, -3, +3, -2,  +1, +1, -3, -2, -2,  -3, -2, +0, +4, -3,  -1, -1, -2,  // V\r
  29         -3, -4, -2, -4, -3,  +1, -2, -2, -3, -3,  -2, -1, -4, -4, -2,  -3, -3, -2, -3,+11,  -1, +2, -3,  // W\r
  30         \r
  31         -1, -1, -1, -1, -1,  -1, -1, -1, -1, -1,  -1, -1, -1, -1, -1,  -1, -1, -1, -1, -1,  -1, -1, -1,  // X\r
  32         -2, -3, -2, -3, -2,  +3, -3, +2, -1, -2,  -1, -1, -2, -3, -1,  -2, -2, -2, -1, +2,  -1, +7, -2,  // Y\r
  33         -1, +2, -4, +2, +5,  -3, -2, +0, -3, +1,  -3, -2, +0, -1, +2,  +0, +0, -1, -2, -3,  -1, -2, +5,  // Z\r
  34         };\r
  35 \r
  36 static int toi_tab[26] =\r
  37         {\r
  38         0,      // A\r
  39         1,      // B\r
  40         2,      // C\r
  41         3,      // D\r
  42         4,      // E\r
  43         5,      // F\r
  44         6,      // G\r
  45         7,      // H\r
  46         8,      // I\r
  47         -1,     // J\r
  48         9,      // K\r
  49         10,     // L\r
  50         11,     // M\r
  51         12,     // N\r
  52         -1,     // O\r
  53         13,     // P\r
  54         14,     // Q\r
  55         15,     // R\r
  56         16,     // S\r
  57         17,     // T\r
  58         17,     // U\r
  59         18,     // V\r
  60         19,     // W\r
  61         20,     // X\r
  62         21,     // Y\r
  63         22,     // Z\r
  64         };\r
  65 \r
  66 static int toi(char c)\r
  67         {\r
  68         c = toupper(c);\r
  69         return toi_tab[c - 'A'];\r
  70         }\r
  71 \r
  72 static int BlosumScore(char c1, char c2)\r
  73         {\r
  74         int i1 = toi(c1);\r
  75         int i2 = toi(c2);\r
  76         return Blosum62[i1][i2];\r
  77         }\r
  78 \r
  79 /***\r
  80 Consider a column with 5 As and 3 Bs.\r
  81 There are:\r
  82         5x4 pairs of As.\r
  83         3x2 pairs of Bs.\r
  84         5x3x2 AB pairs\r
  85         8x7 = 5x4 + 3x2 + 5x3x2 pairs of letters\r
  86 ***/\r
  87 static double BlosumScoreCol(const MSA &a, unsigned uColIndex)\r
  88         {\r
  89         int iCounts[23];\r
  90         memset(iCounts, 0, sizeof(iCounts));\r
  91         const unsigned uSeqCount = a.GetSeqCount();\r
  92         unsigned uCharCount = 0;\r
  93         for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
  94                 {\r
  95                 char c = a.GetChar(uSeqIndex, uColIndex);\r
  96                 if (IsGapChar(c))\r
  97                         continue;\r
  98                 int iChar = toi(c);\r
  99                 ++iCounts[iChar];\r
 100                 ++uCharCount;\r
 101                 }\r
 102         if (uCharCount < 2)\r
 103                 return -9;\r
 104         int iTotalScore = 0;\r
 105         for (int i1 = 0; i1 < 23; ++i1)\r
 106                 {\r
 107                 int iCounts1 = iCounts[i1];\r
 108                 iTotalScore += iCounts1*(iCounts1 - 1)*Blosum62[i1][i1];\r
 109                 for (int i2 = i1 + 1; i2 < 23; ++i2)\r
 110                         iTotalScore += iCounts[i2]*iCounts1*2*Blosum62[i1][i2];\r
 111                 }\r
 112         int iPairCount = uCharCount*(uCharCount - 1);\r
 113         return (double) iTotalScore / (double) iPairCount;\r
 114         }\r
 115 \r
 116 /***\r
 117 Consider a column with 5 As and 3 Bs.\r
 118 A residue of type Q scores:\r
 119         5xAQ + 3xBQ\r
 120 ***/\r
 121 static void AssignColorsCol(const MSA &a, unsigned uColIndex, int **Colors)\r
 122         {\r
 123         int iCounts[23];\r
 124         memset(iCounts, 0, sizeof(iCounts));\r
 125         const unsigned uSeqCount = a.GetSeqCount();\r
 126         unsigned uCharCount = 0;\r
 127         for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
 128                 {\r
 129                 char c = a.GetChar(uSeqIndex, uColIndex);\r
 130                 if (IsGapChar(c))\r
 131                         continue;\r
 132                 int iChar = toi(c);\r
 133                 ++iCounts[iChar];\r
 134                 ++uCharCount;\r
 135                 }\r
 136         int iMostConservedType = -1;\r
 137         int iMostConservedCount = -1;\r
 138         for (unsigned i = 0; i < 23; ++i)\r
 139                 {\r
 140                 if (iCounts[i] > iMostConservedCount)\r
 141                         {\r
 142                         iMostConservedType = i;\r
 143                         iMostConservedCount = iCounts[i];\r
 144                         }\r
 145                 }\r
 146 \r
 147         double dColScore = BlosumScoreCol(a, uColIndex);\r
 148         int c;\r
 149         if (dColScore >= 3.0)\r
 150                 c = 3;\r
 151         //else if (dColScore >= 1.0)\r
 152         //      c = 2;\r
 153         else if (dColScore >= 0.2)\r
 154                 c = 1;\r
 155         else\r
 156                 c = 0;\r
 157 \r
 158         int Color[23];\r
 159         for (unsigned uLetter = 0; uLetter < 23; ++uLetter)\r
 160                 {\r
 161                 double dScore = Blosum62[uLetter][iMostConservedType];\r
 162                 if (dScore >= dColScore)\r
 163                         Color[uLetter] = c;\r
 164                 else\r
 165                         Color[uLetter] = 0;\r
 166                 }\r
 167 \r
 168         for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
 169                 {\r
 170                 char c = a.GetChar(uSeqIndex, uColIndex);\r
 171                 if (IsGapChar(c))\r
 172                         {\r
 173                         Colors[uSeqIndex][uColIndex] = 0;\r
 174                         continue;\r
 175                         }\r
 176                 int iLetter = toi(c);\r
 177                 if (iLetter >= 0 && iLetter < 23)\r
 178                         Colors[uSeqIndex][uColIndex] = Color[iLetter];\r
 179                 else\r
 180                         Colors[uSeqIndex][uColIndex] = 0;\r
 181                 }\r
 182         }\r
 183 \r
 184 void AssignColors(const MSA &a, int **Colors)\r
 185         {\r
 186         const unsigned uColCount = a.GetColCount();\r
 187         for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
 188                 AssignColorsCol(a, uColIndex, Colors);\r
 189         }\r