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

Last change on this file was 10390, checked in by aboeckma, 11 years ago

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