!!!課題1-3（塩基使用頻度その3）

!!!0.1　課題1-3

二連続塩基(dinucleotide)の頻度解析を行うようにプログラムを変更する。マイコプラズマ全ゲノムDNA配列の二連続塩基の頻度パターンを調べ、1-2の結果からの期待値と比較、考察する。

この問題の中核は「期待値をどうやって求めるか」と「実際の値と期待値をどう比較するか」ということになると思います。観測値(observation)を期待値(expectation)で割ったものをO/E値と言います。これを使うと例えば「AAは実際には期待値の半分しかなかった。しかしCTは2倍以上ある。何か理由があるに違いない」などという議論ができますので、今回は比較の方法としてこれを採用しましょう。

!!!1　　各dinucleotideの頻度を数える

!!1.1　データ変数の用意

 my %dinuc;   #dinucleotide の数を数える
 my %diexp;   #dinucleotide の期待値
 my %diobs;   #dinucleotide の実頻度
 my %oe;       #dinucleotide のO/E値

こういう変数が必要になることが予想されます。ここで%が先頭についたものはハッシュ（連想配列）というものです。ハッシュとは、基本的には先頭に@がつく配列と似たような機能を持っていますが、配列と違う点は引数に文字列が使える点です。ですから、例えば配列では

 @array = (1, 2, 3);

としたら$array[1]は2となりますが、ハッシュの場合は

 %hash = {'string'=>1, 'message'=>2, 'line'=>3}

としたら$hash{'message'}が2となるように、任意の文字列を引数に使えます。このハッシュを使えば引数としてdinucleotide自身を使えるようになりますので、プログラムが非常に楽になります。Perlではこのハッシュを使う機会が多いので、是非しっかり覚えてください。

!!1.2　カウント部分の改良

基本は変わらないのですが、二連続塩基ですので1文字置換をするtr///を使って数えることはできません。複数置換のs///では置換された数を返さないのでこれを使うこともできません。しかたがないので、ここではfor文を使って2文字の並びを少しずつずらして読んでいきましょう。ここでは塩基配列の最初から2文字の並びを一文字ずつずらして読むのですから、例えば"atgcggctg"という並びなら一つ目が"at"、二つ目が"tg"というようになります。つまりこれは[一文字目]から[配列の長さ-1文字目]までfor文で繰り返すことになります。ここで注意しなければいけないことは、Perlでは先頭の番号は1ではなく0から始まる点です。つまり、for文で繰り返す場合は[0文字目]から[配列の長さ-2文字目]までとなります。

文字列から一部分、例えばこの場合は2連続塩基を切り出すには、substr()関数を使います。substr()は

 $parts = substr($seq, 0, 2);

のように（文字列, 切り出す場所, 切り出す文字数）という引数を与えて使います。
　
さて、ここでハッシュが大変役に立ちます。ハッシュの引数は文字列でいいので、

 $diobs{$parts} ++;

としてやることで、$partsの中身が何であるかを気にせずにカウントすることができます。ここで、あくまでカウントしているのはハッシュの中身でありハッシュではないので、

 %diobs{$parts} ++;

ではなく

 $diobs{$parts} ++;

となっている点に注意してください。

!!!2　　O/E値を求める

!!2.1　期待値を求める：基本

期待値ってなんでしょう。文字通り「期待される値」です。たとえばdinucleotideの場合何も考えなければ、4×4=16通りあるわけですから各dinucleotideの期待値は一律1/16です。しかし、1-2で見たように生物はa,t,g,cをそれぞれ均等に使っているわけではありません。したがって、もともと使用頻度が低い塩基（仮にCとします）を含む組合せの期待値も当然低くならなければいけません。

したがって、単に一律1/16とするのではなく、各生物ごとの期待値というもので議論する必要があるわけです。

この時、期待値を求める式は以下のようになります。

1文字目の塩基の頻度 × 2文字目の塩基の頻度

したがってまず、a,t,g,cを頻度に直すことから考えましょう。頻度は小数で表したいですね。うまく変換して下さい。

!!2.2　期待値を求める：応用

この部分の計算は

 $diexp{'aa'} = $percent_a * $percent_a;
 $diexp{'at'} = $percent_a * $percent_t;
                     .
                     .
                     .

のように、普通に書くと式が16本できてしまいます。これをどうにかしてまとめてみましょう。PerlにはCやJavaにはないforeach文というものがありますが、これを使えばどうにかなりそうです。

foreach 文は、配列やハッシュの全ての要素を順番に計算します。例えば

 foreach $content (@array){
 	print $content;
 }

とすれば、配列の各要素が順次 $content の中に入れられるので、次々とその中身が出力されます。ハッシュの場合も同じように

 foreach $key (sort keys %hash){
 	print "$key = " , $hash{$key};
 }

とすることで毎回ハッシュのキー（引数となる文字列）が$keyに入れられます。sort keysというのは、foreachで計算する順序を、キーをASCII順に並べた順番にする、という意味です。

それでは、foreach を使って、%diobsの全ての要素を順次とりだしながら、毎回$keyに入力されるキー（ニ連続塩基）の期待値とO/E値も同時に計算されるようにプログラムを改良してみましょう。

!!2.3　O/E値を求める

O/E値はその名の通り、実測値/期待値です。同じforeach文の中であらかじめ計算してある期待値と実測値をもとにO/E値を計算し、これもハッシュに代入しましょう。

!!!3　　出力

ここまで来れば%oeというハッシュにaaから順番に求めるデータが格納されているはずです。後は出力させるだけですね。

「組合せ、O/E値、実際値、期待値」と1行に出せばいいですね。最後の2つはおまけです。これもforeach文の中でやってしまいましょう。

［出力例］

      O/E: Observation/Expectation
 aa  *.22:        0.*5/0.*2
 at  *.77:        *.*9/*.*2
                  .
                  .
                  .

前回同様printfを使って整形します。

 print ("     O/E: Observation/Expectation\n");
 printf ("%s  %.2f:       %.2f/%.2f\n", $key, $oe{$key}, $diobs{$key}, $diexp{$key});

ちなみに出力例を出したプログラムは上のようになっています。こういうやり方もあるということを一応知っておいて損はないでしょう。