blusom 矩阵计算过程

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

模块替换矩阵BLOSUM(BLOcks SUbstitution Matrix),用于计算氨基酸替换频率。
attachments-2019-09-AvBH8rOa5d770fffeb913.pngYouTube上看到一个关于讲解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";

    }

}

结果可输出类似下图矩阵


image.png


  • 发表于 2019-09-10 11:13
  • 阅读 ( 63 )
  • 分类:编程语言

你可能感兴趣的文章

相关问题

0 条评论

请先 登录 后评论
不写代码的码农
yami

1 篇文章

作家榜 »

  1. 合肥国肽生物 113 文章
  2. 祝让飞 104 文章
  3. 刘永鑫 64 文章
  4. 生信分析流 47 文章
  5. SXR 44 文章
  6. 调研图 38 文章
  7. 张海伦 31 文章
  8. 爽儿 25 文章