Bioinformatics_Basics/BLOSUM.m at main · zandigohar/Bioinformatics_Basics · GitHub

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
clc;
close all;
clear;

S1='ALIAA';
S2='ALLLA';
S3='AILIL';

S=[S1;S2;S3];%;S4;S5];

[m,n]=size(S);

SS='';
for i=1:m
    SS=strcat(SS,S(i,:));
end
elements=unique(SS);
p=length(elements);

nij=zeros(1,p);
cijk=zeros(p,p);
for k=1:n
    %step1
    for j=1:p
        temp=S(:,k)';
        nij(k,j)=length(find(temp==elements(j)));
    end
    %step2
    for i=1:p
        for j=1:p
            if i==j
                cij(i,j)=nij(k,i)*(nij(k,i)-1)/2;
            else
                cij(i,j)=nij(k,i)*nij(k,j);
            end
        end
    end
    cij
    cijk=cijk+cij;
end

T=n*(m*(m-1)/2);
qij=cijk/T;

for i=1:p
    pr(i)=length(find(S==elements(i)))/(m*n);
end

for i=1:p
    for j=1:p
        if i==j
            eij(i,j)=pr(i)^2;
        else
            eij(i,j)=2*pr(i)*pr(j);
        end
    end
end

SSA=round(2*log2(qij./eij));
disp(triu(SSA));


clc;


fprintf('\t');
for i=1:p
    fprintf("%s\t",elements(i));
end
fprintf('\n');
for i=1:p
    fprintf("%s \t",elements(i))
    for j=1:i-1
        fprintf("\t")
    end
    for j=i:p
        if SSA(i,j)==-Inf
            fprintf("-\t");
        else
            fprintf("%d\t",SSA(i,j));
        end
    end
    fprintf('\n');
end