-
Notifications
You must be signed in to change notification settings - Fork 30
Expand file tree
/
Copy pathblosum.cpp
More file actions
92 lines (80 loc) · 6.11 KB
/
blosum.cpp
File metadata and controls
92 lines (80 loc) · 6.11 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
#include "muscle.h"
#include "m3alnparams.h"
#include "alpha.h"
// src/blosum/matrices/blosumNN.sij
// Converted to C syntax by blosum_qij_reformat.py
float Blosum62_sij[20][20] =
{
// A C D E F G H I K L M N P Q R S T V W Y
/* A */ { 1.9646f, -0.2043f, -0.8767f, -0.4319f, -1.1050f, 0.0798f, -0.8126f, -0.6609f, -0.3670f, -0.7323f, -0.4676f, -0.7654f, -0.4071f, -0.4020f, -0.7068f, 0.5579f, -0.0227f, -0.0947f, -1.2634f, -0.8820f, }, // A
/* C */ { -0.2043f, 4.2911f, -1.7300f, -1.8062f, -1.1877f, -1.2502f, -1.4939f, -0.6138f, -1.5182f, -0.6387f, -0.7099f, -1.3299f, -1.3976f, -1.4509f, -1.6946f, -0.4375f, -0.4333f, -0.4038f, -1.1521f, -1.2036f, }, // C
/* D */ { -0.8767f, -1.7300f, 2.8871f, 0.7552f, -1.7419f, -0.6568f, -0.5595f, -1.5606f, -0.3509f, -1.8028f, -1.5293f, 0.6358f, -0.7401f, -0.1567f, -0.8029f, -0.1305f, -0.5254f, -1.5713f, -2.1072f, -1.5325f, }, // D
/* E */ { -0.4319f, -1.8062f, 0.7552f, 2.4514f, -1.5962f, -1.0551f, -0.0588f, -1.5972f, 0.3877f, -1.4232f, -0.9990f, -0.1340f, -0.5581f, 0.9273f, -0.0577f, -0.0735f, -0.4316f, -1.2211f, -1.4177f, -1.0102f, }, // E
/* F */ { -1.1050f, -1.1877f, -1.7419f, -1.5962f, 3.0230f, -1.5537f, -0.6171f, -0.0804f, -1.5393f, 0.2074f, 0.0063f, -1.4970f, -1.7986f, -1.5822f, -1.3932f, -1.1845f, -1.0538f, -0.4245f, 0.4588f, 1.4696f, }, // F
/* G */ { 0.0798f, -1.2502f, -0.6568f, -1.0551f, -1.5537f, 2.7816f, -1.0204f, -1.8624f, -0.7640f, -1.8135f, -1.3383f, -0.2114f, -1.0668f, -0.8926f, -1.1521f, -0.1462f, -0.7877f, -1.5694f, -1.2457f, -1.5199f, }, // G
/* H */ { -0.8126f, -1.4939f, -0.5595f, -0.0588f, -0.6171f, -1.0204f, 3.7555f, -1.6158f, -0.3605f, -1.3934f, -0.7756f, 0.2892f, -1.0805f, 0.2240f, -0.1249f, -0.4408f, -0.8429f, -1.5587f, -1.1711f, 0.8463f, }, // H
/* I */ { -0.6609f, -0.6138f, -1.5606f, -1.5972f, -0.0804f, -1.8624f, -1.6158f, 1.9993f, -1.3351f, 0.7608f, 0.5634f, -1.6085f, -1.3783f, -1.3848f, -1.4951f, -1.1741f, -0.3588f, 1.2735f, -1.2903f, -0.6657f, }, // I
/* K */ { -0.3670f, -1.5182f, -0.3509f, 0.3877f, -1.5393f, -0.7640f, -0.3605f, -1.3351f, 2.2523f, -1.2234f, -0.6774f, -0.0895f, -0.5068f, 0.6363f, 1.0544f, -0.1017f, -0.3348f, -1.1312f, -1.4782f, -0.9100f, }, // K
/* L */ { -0.7323f, -0.6387f, -1.8028f, -1.4232f, 0.2074f, -1.8135f, -1.3934f, 0.7608f, -1.2234f, 1.9247f, 0.9959f, -1.6895f, -1.4300f, -1.0670f, -1.0773f, -1.2213f, -0.5987f, 0.3942f, -0.8159f, -0.5310f, }, // L
/* M */ { -0.4676f, -0.7099f, -1.5293f, -0.9990f, 0.0063f, -1.3383f, -0.7756f, 0.5634f, -0.6774f, 0.9959f, 2.6963f, -1.0754f, -1.2382f, -0.2105f, -0.6836f, -0.7404f, -0.3331f, 0.3436f, -0.7124f, -0.4974f, }, // M
/* N */ { -0.7654f, -1.3299f, 0.6358f, -0.1340f, -1.4970f, -0.2114f, 0.2892f, -1.6085f, -0.0895f, -1.6895f, -1.0754f, 2.8266f, -1.0002f, 0.0008f, -0.2199f, 0.3005f, -0.0230f, -1.4382f, -1.8480f, -1.0409f, }, // N
/* P */ { -0.4071f, -1.3976f, -0.7401f, -0.5581f, -1.7986f, -1.0668f, -1.0805f, -1.3783f, -0.5068f, -1.4300f, -1.2382f, -1.0002f, 3.6823f, -0.6410f, -1.0543f, -0.4045f, -0.5376f, -1.1744f, -1.8271f, -1.4599f, }, // P
/* Q */ { -0.4020f, -1.4509f, -0.1567f, 0.9273f, -1.5822f, -0.8926f, 0.2240f, -1.3848f, 0.6363f, -1.0670f, -0.2105f, 0.0008f, -0.6410f, 2.6426f, 0.4914f, -0.0506f, -0.3377f, -1.0992f, -0.9732f, -0.7105f, }, // Q
/* R */ { -0.7068f, -1.6946f, -0.8029f, -0.0577f, -1.3932f, -1.1521f, -0.1249f, -1.4951f, 1.0544f, -1.0773f, -0.6836f, -0.2199f, -1.0543f, 0.4914f, 2.7367f, -0.3824f, -0.5612f, -1.2513f, -1.3397f, -0.8469f, }, // R
/* S */ { 0.5579f, -0.4375f, -0.1305f, -0.0735f, -1.1845f, -0.1462f, -0.4408f, -1.1741f, -0.1017f, -1.2213f, -0.7404f, 0.3005f, -0.4045f, -0.0506f, -0.3824f, 1.9422f, 0.6906f, -0.8231f, -1.3759f, -0.8429f, }, // S
/* T */ { -0.0227f, -0.4333f, -0.5254f, -0.4316f, -1.0538f, -0.7877f, -0.8429f, -0.3588f, -0.3348f, -0.5987f, -0.3331f, -0.0230f, -0.5376f, -0.3377f, -0.5612f, 0.6906f, 2.2727f, -0.0278f, -1.2145f, -0.8030f, }, // T
/* V */ { -0.0947f, -0.4038f, -1.5713f, -1.2211f, -0.4245f, -1.5694f, -1.5587f, 1.2735f, -1.1312f, 0.3942f, 0.3436f, -1.4382f, -1.1744f, -1.0992f, -1.2513f, -0.8231f, -0.0278f, 1.8845f, -1.4171f, -0.6038f, }, // V
/* W */ { -1.2634f, -1.1521f, -2.1072f, -1.4177f, 0.4588f, -1.2457f, -1.1711f, -1.2903f, -1.4782f, -0.8159f, -0.7124f, -1.8480f, -1.8271f, -0.9732f, -1.3397f, -1.3759f, -1.2145f, -1.4171f, 5.2520f, 1.0771f, }, // W
/* Y */ { -0.8820f, -1.2036f, -1.5325f, -1.0102f, 1.4696f, -1.5199f, 0.8463f, -0.6657f, -0.9100f, -0.5310f, -0.4974f, -1.0409f, -1.4599f, -0.7105f, -0.8469f, -0.8429f, -0.8030f, -0.6038f, 1.0771f, 3.2975f, }, // Y
};
void GetSubstMx_Letter_Blosum(uint PctId, float Mx[20][20])
{
#define X(x) \
case x: \
for (unsigned i = 0; i < 20; ++i) \
for (unsigned j = 0; j < 20; ++j) \
Mx[i][j] = Blosum##x##_sij[i][j]; \
return;
switch (PctId)
{
X(62);
}
#undef X
Die("SetSubstMx_Letter_Blosum(%u)", PctId);
}
void GetGapParams_Blosum(uint PctId, uint n, float *ptrGapOpen, float *ptrCenter)
{
#define X(p, i, G, C) if (PctId == p && i == n) { *ptrGapOpen = float(G); *ptrCenter = float(C); return; }
X(90, 0, -7.3333335, 1.2)
X(90, 1, -8.1662216, 1.0788642)
X(90, 2, -6.7398319, 1.0459337)
X(90, 3, -7.0647068, 1.2546233)
X(80, 0, -6.6666665, 0.99999994)
X(80, 1, -7.2274466, 0.91091353)
X(80, 2, -7.6157303, 0.86217165)
X(80, 3, -7.1673636, 0.85966408)
X(70, 0, -6.2208495, 0.88161403)
X(70, 1, -7.3177958, 0.70952064)
X(70, 2, -7.1693735, 0.93325645)
X(70, 3, -6.7926803, 0.71609467)
X(62, 0, -6, 0.79999995)
X(62, 1, -5.6413326, 0.71837389)
X(62, 2, -6.6825562, 0.59377569)
X(62, 3, -5.574501, 0.66151822)
#undef X
Die("GetGapParams_Blosum(%u, %u)", PctId, n);
}
float GetBlosumScoreChars(byte a, byte b)
{
uint ia = g_CharToLetterAmino[a];
uint ib = g_CharToLetterAmino[b];
if (ia >= 20 || ib >= 20)
return 0;
return Blosum62_sij[ia][ib];
}
float GetBlosumScoreLetters(byte ia, byte ib)
{
if (ia >= 20 || ib >= 20)
return 0;
return Blosum62_sij[ia][ib];
}