Perl を利用してアミノ酸と塩基の出現頻度を計算する

アミノ酸・塩基の出現頻度

Perl スクリプトを利用して、塩基あるいはアミノ酸の出現頻度を計算するとき、ハッシュを利用すると便利。おおよその処理は次のように行う。

  1. $i を 0 とする。

  2. $i が文字列の最後に達するまで以下の処理を繰り返す。

    1. 入力された塩基あるいはアミノ酸配列 $seq に対して、位置 $i から長さ $leng だけ切り出して $str に保存する。

    2. $str はハッシュリファレンス $count のキーとして存在するかどうかを調べる。存在しなければ $count->{$str} = 0 と初期化する。

    3. $count->{$str} に 1 を加える。

    4. $i に 1 を加える。

  3. 最後に合計値で割り、割合を計算する。

use strict;
 
sub countString {
  my ($seq, $leng) = @_;

  my $count = {};

  $seq =~ tr/A-Z/a-z/;
  my $seq_len = length($seq);
 
  for (my $i = 0; $i < $seq_len; $i++) {
    my $str = substr($seq, $i, $leng);
    if (length($str) == $leng) {
      $count->{absolute}->{$str}++;
    }
  }

  my $count_sum = $seq_len - $leng + 1;
  foreach my $key (sort keys %{$count->{absolute}}) {
    $count->{relative}->{$key} = sprintf("%.5f", $count->{absolute}->{$key} / $count_sum);
  }
  return $count;
}
 

my $nt_seq = 'AGTAGTA';
my $nt_result = &countString($nt_seq, 3);
foreach my $key (keys %{$nt_result->{relative}}) {
  print "$key :  $nt_result->{relative}->{$key}\n";
}
## tag :  0.20000
## gta :  0.40000
## agt :  0.40000


my $aa_seq = 'TPEMPVLENK';
my $aa_result = &countString($aa_seq, 1);
foreach my $key (keys %{$aa_result->{relative}}) {
  print "$key :  $aa_result->{relative}->{$key}\n";
}
## e :  0.20000
## n :  0.10000
## v :  0.10000
## m :  0.10000
## l :  0.10000
## p :  0.20000
## k :  0.10000
## t :  0.10000