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
93
94
95
96
97
98
99
|
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "clustalw.h"
/*
* Prototypes
*/
/*
* Global variables
*/
extern sint max_aa,gap_pos1,gap_pos2;
void calc_prf1(sint **profile, char **alignment, sint *gaps,
sint matrix[NUMRES][NUMRES],
sint *seq_weight, sint prf_length, sint first_seq, sint last_seq)
{
sint **weighting, sum2, d, i, res;
sint numseq;
sint r, pos;
int f;
float scale;
weighting = (sint **) ckalloc( (NUMRES+2) * sizeof (sint *) );
for (i=0;i<NUMRES+2;i++)
weighting[i] = (sint *) ckalloc( (prf_length+2) * sizeof (sint) );
numseq = last_seq-first_seq;
sum2 = 0;
for (i=first_seq; i<last_seq; i++)
sum2 += seq_weight[i];
for (r=0; r<prf_length; r++)
{
for (d=0; d<=max_aa; d++)
{
weighting[d][r] = 0;
for (i=first_seq; i<last_seq; i++)
if (d == alignment[i][r]) weighting[d][r] += seq_weight[i];
}
weighting[gap_pos1][r] = 0;
for (i=first_seq; i<last_seq; i++)
if (gap_pos1 == alignment[i][r]) weighting[gap_pos1][r] += seq_weight[i];
weighting[gap_pos2][r] = 0;
for (i=first_seq; i<last_seq; i++)
if (gap_pos2 == alignment[i][r]) weighting[gap_pos2][r] += seq_weight[i];
}
for (pos=0; pos< prf_length; pos++)
{
if (gaps[pos] == numseq)
{
for (res=0; res<=max_aa; res++)
{
profile[pos+1][res] = matrix[res][gap_pos1];
}
profile[pos+1][gap_pos1] = matrix[gap_pos1][gap_pos1];
profile[pos+1][gap_pos2] = matrix[gap_pos2][gap_pos1];
}
else
{
scale = (float)(numseq-gaps[pos]) / (float)numseq;
for (res=0; res<=max_aa; res++)
{
f = 0;
for (d=0; d<=max_aa; d++)
f += (weighting[d][pos] * matrix[d][res]);
f += (weighting[gap_pos1][pos] * matrix[gap_pos1][res]);
f += (weighting[gap_pos2][pos] * matrix[gap_pos2][res]);
profile[pos+1][res] = (sint )(((float)f / (float)sum2)*scale);
}
f = 0;
for (d=0; d<=max_aa; d++)
f += (weighting[d][pos] * matrix[d][gap_pos1]);
f += (weighting[gap_pos1][pos] * matrix[gap_pos1][gap_pos1]);
f += (weighting[gap_pos2][pos] * matrix[gap_pos2][gap_pos1]);
profile[pos+1][gap_pos1] = (sint )(((float)f / (float)sum2)*scale);
f = 0;
for (d=0; d<=max_aa; d++)
f += (weighting[d][pos] * matrix[d][gap_pos2]);
f += (weighting[gap_pos1][pos] * matrix[gap_pos1][gap_pos2]);
f += (weighting[gap_pos2][pos] * matrix[gap_pos2][gap_pos2]);
profile[pos+1][gap_pos2] = (sint )(((float)f / (float)sum2)*scale);
}
}
for (i=0;i<NUMRES+2;i++)
weighting[i]=ckfree((void *)weighting[i]);
weighting=ckfree((void *)weighting);
}
|
