gibbs sampling

3月に入って netwalker を購入しました。

電車通勤の時間を利用してコードを書きたいなと思っていたのですが、なかなかよろしい感じです。
といってもキーボードは使いやすいとは言い難く、長いコードは書きたくないですが...

ということで、netwalker で作成したプログラム第一号を公開します。
第一号は gibbs sampling を使ったモチーフ抽出アルゴリズムの実装です。

前々から gibbs sampling について調べてはいたのですが、適度な練習問題が無く実装はしていませんでした。
今月になって購入したバイオインフォマティクスの数理とアルゴリズムにちょうど良い例があったので、それを実装してみました。

そもそもモチーフって何という話ですが、一言で言うとタンパク質の配列パターンの事を指すそうです。
私は専門ではないので詳しい事は分かりませんが、ここでは与えられた複数のタンパク質配列から、特徴的な配列パターンを見つけてくる事をモチーフ抽出と言ってよさそうです。

具体的にどうモチーフを抽出するかですが、載っている例では、各タンパク質A、C、T、Gの出現確率をP(A)、P(C)、P(T)、P(G)、部分配列中の位置 j に各タンパク質が現れる確率を f_j(A)、f_j(C)、f_j(T)、f_j(G)とした時に、相対エントロピーを最大化するように部分配列を選んでくる事でモチーフを抽出します。

これをローカルマルチプルアラインメントと呼ぶそうですが、この問題はまともにやると NP 困難なので、gibbs sampling を使って効率的に解くよ、というお話しになっています。

という事で、アルゴリズムは以下の通り。

1. 各入力配列 s_k から, 長さ L の部分配列 t_k を一様ランダムに選択 (k = 1,...,n)
2. 入力配列から, s_i を一様ランダムに選択
3. f'_j(a) を t_i 意外の t_k から計算したモチーフ領域の j 番目の列に置ける文字 a の頻度とする ( f'_j(a) = (|{k|t_k[j] = a, k != i}|) / (n-1))
4. s_i から長さ L の部分配列 t'_i を, Πj f'_j(t'_i[j])/p(t'_i[j]) に比例する確率でランダムに選択
5. t_i を t'_i で置き換える
6. 2-5 を十分な回数繰り返す

ステップ 3. の頻度は確率の間違えかな？と思うのですが、一応そのまま書いておきます。
ステップ 4. では求めた事後分布に従って新しくモチーフ候補を選択し直すわけですが、ここはランダムで選択しないと初期値から動かなくなるので、まず計算したスコアの配列を確率分布として使えるように総和が1となるように正規化してから、PRMLの下巻に紹介されている棄却サンプリングを使いました。
提案分布には一様分布を使っています。

コードと実行例は以下。

#!/usr/local/bin/perl

use warnings;
use strict;

my $L = 5;
my @dataset;
my %wordset;
my $words = 0;

while(<>)
{
   chomp;
   my $line = $_;
   my @sequence = split(//,$line);
   push @dataset,\@sequence;
}

# init words and probs
$words = &setwords(\%wordset, \@dataset);
&initwordprobs($words,\%wordset);
my $datasize = @dataset;

print "dataset\n";
foreach my $ar (@dataset)
{
   print "@$ar\n";
}

# init motif set
my @pos;
foreach my $ar (@dataset)
{
   my $p = int rand(@$ar-$L);
   push @pos, $p;
   #print $p,"\n";
}

# output first
print "init\n";
for (my $i = 0; $i < @dataset; $i++)
{
   my $ar = $dataset[$i];
   my $p = $pos[$i];
   for (my $j = $p; $j < $p+$L; $j++)
   {
      print "$ar->[$j] ";
   }
   print "\n";
}

# iterration
for (my $iter = 0; $iter < 50; $iter++)
{
   #print "epoch: $iter\n";
   my $id = int rand($datasize);
   my %wordsfreq;
# wordset * position
# +---+-------------+
# |   |a|b|c|d|e|f|g|
# +---+-------------+
# | 0 |             |
# | 1 |             |
# | 2 |  frequency  |
# | . |             |
# | L |             |
# +---+-------------+
   &initwordsfreq(\%wordset,$L,\%wordsfreq);
# count word frequency
   &countwordfreq($id,\@dataset,\%wordsfreq,\@pos);
# pick up new motif i
   my $ar = $dataset[$id];
   $pos[$id] = &pickup($id,$ar,$L,\%wordsfreq, \%wordset);
}

# output
print "result\n";
for (my $i = 0; $i < @dataset; $i++)
{
   my $ar = $dataset[$i];
   my $p = $pos[$i];
   for (my $j = $p; $j < $p+$L; $j++)
   {
      print "$ar->[$j] ";
   }
   print "\n";
}

# subroutine
sub pickup
{
   my $id = shift;
   my $ar = shift;
   my $L = shift;
   my $wordsfreqhr = shift;
   my $wordsethr = shift;

   my @scores;
   my $total = 0;
   for (my $p = 0; $p <= @$ar-$L; $p++)
   {
      my $score = 1;
      my $index = 0;
      for (my $j = $p; $j < $p+$L; $j++)
      {
         my $w = $ar->[$j];
         $score *= $wordsfreqhr->{$index}->{$w}/$wordsethr->{$w};
         $index++;
      }
      push @scores, $score;
      $total += $score;
   }
   return int rand(@scores) if $total == 0;
   for(my $i = 0; $i < @scores; $i++)
   {
      $scores[$i] /= $total;
   }
   #print "@scores\n";
   my $flg = 1;
   while($flg)
   {
      my $npos = int rand(@scores);
      my $s = rand(1);
      if ($s <= $scores[$npos])
      {
         $flg = 0;
         return $npos;
      }
   }
}
sub countwordfreq
{
   my $id = shift;
   my $datasetar = shift;
   my $wordsfreqhr = shift;
   my $posar = shift;

   for (my $i = 0; $i < @$datasetar; $i++)
   {
      next if ($i == $id);
      my $ar = $datasetar->[$i];
      for (my $j = $posar->[$i]; $j < $L+$posar->[$i]; $j++)
      {
         my $p = $j-$posar->[$i];
         my $w = $ar->[$j];
         $wordsfreqhr->{$p}->{$w}++;
      }
   }
   for (my $l = 0; $l < $L; $l++)
   {
      foreach my $k (keys %{$wordsfreqhr->{$l}})
      {
         $wordsfreqhr->{$l}->{$k} /= (@$datasetar-1);
      }
   }
}

sub initwordsfreq
{
   my $wordsethr = shift;
   my $L = shift;
   my $wordfhr = shift;
   for (my $i = 0; $i < $L; $i++)
   {
      my %freqonpos;
      foreach my $k (keys %$wordsethr)
      {
         $freqonpos{$k} = 0;
      }
      $wordfhr->{$i} = \%freqonpos;
   }
}
sub initwordprobs
{
   my $words = shift;
   my $wordsethr = shift;

   foreach my $k (keys %$wordsethr)
   {
      $wordsethr->{$k} = 1/$words;
   }
}

sub setwords
{
   my $wordshr = shift;
   my $dataar = shift;
   my $words = 0;

   foreach my $ar (@$dataar)
   {
      foreach my $w (@$ar)
      {
         unless (defined $wordshr->{$w})
         {
            $wordshr->{$w} = defined;
            $words++;
         }
      }
   }
   return $words;
}

サンプルデータ1

c a t c g a a a a a a
a a a a c a t c g a a
a a c a t c g a a a a
a a a a a a c a t c g

サンプルデータ2

t c a g a g a t t c g t a g
t t c a t t c g g g c g t
c g a t t c g c g a c t c
t g a a t t c g g a a

実行例1

# perl gibbs.pl < sample1.txt
dataset
c a t c g a a a a a a
a a a a c a t c g a a
a a c a t c g a a a a
a a a a a a c a t c g
init
g a a a a 
a a a c a 
a t c g a 
a a a a a 
result
c a t c g 
c a t c g 
c a t c g 
c a t c g 

実行例2

# perl gibbs.pl < sample2.txt
dataset
t c a g a g a t t c g t a g
t t c a t t c g g g c g t
c g a t t c g c g a c t c
t g a a t t c g g a a
init
t c a g a 
c a t t c 
t c g c g 
a a t t c 
result
a t t c g 
a t t c g 
a t t c g 
a t t c g 

上手い事特徴的な文字列を抽出できる事が確認できました。