模块替换矩阵BLOSUM(BLOcks SUbstitution Matrix),用于计算氨基酸替换频率。
YouTube上看到一个关于讲解blusom 矩阵计算过程的视频，给出了一个例子：
seq1 B A B A 

seq2 A A A C

seq3 A A C C

seq4 A A B A

seq5 A A C C 

seq6 A A B C 

简单而言就是指氨基酸对的出现频率除以期望频率后取对数。

（1）计算blocks 中氨基酸对的数目：等于每一列的频率值之和

AA=C52+C62+0+C21=26
AB=C51C11+0+C11C31+0=8
AC=0+0+C11C21+C21C41=10
BB=0+0+C32=3
BC=0+0+C31C21=6
CC=0+0+C22+C42=7
(2) 计算氨基酸对出现频率：
Total=26+8+10+3+6+7=60;
p_AA=26/60;
...
p_CC=7/60;
(3）计算氨基酸对的期望频率
a)统计各个氨基酸的总数A=5+6+1+2=14；B=1+0+3+0=4；C=0+0+2+4=6
b)氨基酸总数=A+B+C=14+4+6=24
c)计算期望频率e
e_AA=(14/2414/24)2；
e_AB=14/244/24;
...
e_CC=6/246/242

(4) 计算blusom score
AA_score=2log2(p_AA/e_AA)
AB_score=2log2(p_AB/e_AB)
...
CC_score=2*log2(p_CC/e_CC)

理解了计算过程之后，对于具体对一个group的蛋白，可通过序列比对得到.aln 格式文件(比对时候参数选择Identity matrix).

下面是编写的perl脚本，可用来计算score matrix (输入文件不能含有gap "-"，需要从比对结果选取blocks)

open M,"$ARGV[0]"; ##input clustw file *.aln

while (my $line=<M>){

       chomp $line;

       next if ($line=~/CLUSTAL.*/);

       next if ($line!~/\w.*/);

       $num++

       @m=split(/\s+/,$line);

       $seq=$m[1];

       $hash{$m[0]}{$num}=$seq;

       $len=length ($seq);

    }

close M;

foreach $k (keys %hash){

    $seq_all=();

    foreach $k1 (keys %{$hash{$k}}){

         $seq_all.=$hash{$k}{$k1};#connect into one line for each sequence

    }

    @amino=split(//,$seq_all);#split seq into each aminofor($i=0;$i<scalar @amino; $i++){

        $hash_str{$i}{$k}=$amino[$i];

    } 

}

                        

@amino_acid=(A,C,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y);#20种氨基酸简写

%freq_AB=();

%hash_aa_all=(); ## amino num for all column

$amino_num=(); 

foreach $col (keys %hash_str){ 

    @pro=();

    %hash_aa=();# amino num for per column 

    foreach $col_seq( keys %{$hash_str{$col}}){          

       push @pro, $hash_str{$col}{$col_seq}; ## put all amino each column into an array

    }  

    $amino_num+=scalar @pro;#stat number for all amino

    foreach $aa (@pro){

      $hash_aa{$aa}++; ##stat num of each amino in one column

      $hash_aa_all{$aa}++;#stat num of each amino for all column

    }

     for(my $i=0;$i<20;$i++){

        $A=$amino_acid[$i];

        for(my $j=$i;$j<20;$j++){

          $B=$amino_acid[$j];         

          $AB=$amino_acid[$i].$amino_acid[$j];  #two amino if($i==$j){

              if($hash_aa{$A}>1){

                 $freq=($hash_aa{$A}*($hash_aa{$A}-1))/2;

              }

             else{

             $freq=0;

             }

           }

         elsif($i!=$j){

            if(($hash_aa{$A}>0)&($hash_aa{$B}>0)){           

                $freq=$hash_aa{$A}*$hash_aa{$B};

            }

           else{$freq=0;}

         }

       $freq_all+=$freq;# freq for two amino permutation

       $freq_AB{$AB}+=$freq; # freq for AB amino

    }

  }

}

print join("\t",@amino_acid);

print "\n";

foreach my $k (sort keys %freq_AB){

    my @str =split(//,$k);

    my $A=$str[0];

    my $B=$str[1];

    $p_AB=$freq_AB{$k}/$freq_all; ##observed probability for AB

   if ($p_AB>0){

      if(($hash_aa_all{$A}>0)&&($hash_aa_all{$B}>0)){

        if($A eq $B){

          $e_AB=($hash_aa_all{$A}/$amino_num)*($hash_aa_all{$B}/$amino_num);

        }       

        elsif($A ne $B){

          $e_AB=($hash_aa_all{$A}/$amino_num)*($hash_aa_all{$B}/$amino_num)*2;

        }

      $log_odd=(log($p_AB/$e_AB)/log(2))*2;

      }

    }

    else{$log_odd=0;}

    if ($B eq "Y" ){

      print "$log_odd\t$A\n";

       $n++;#number of "\t"

      $num_t="\t"x"$n"; #matrix align

      print"$num_t";

    }

    else{

      print "$log_odd\t";

    }

}

结果可输出类似下图矩阵