DNA配列の統計的性質 第7回 〜DNA walking〜

前回までの記事で、GenBankからファイルをとってきて、そのDNA配列データのヌクレオチド頻度を計算することができるようになりました。さらにスクリプトに手を加えることで、あるゲノムDNA配列に対してGC%の分布をも描けるようになりました。

GC%の分布のほかにもDNA配列の並びのパターンといいますか性質を調べるのに様々な測定量があります。たとえば、GC-skewですが、これは一本鎖DNAにおけるG含有量とC含有量のバイアスを表す指標で、(C-G)/(C+G) の式で表される測定量で、バクテリアの種によってはバイアス傾向が変わる個所が複製開始・終結地点と一致することが多いという指摘がなされていまして、積極的に生物学的意義付けもされています[2]。さて今回は、別の視点・方向からDNA配列の性質に迫ってみたいと思います。例によって前回のスクリプトを叩き台にして発展させます。

この連載の最初の方でPerlの導入に伴い、また世界物理年を祝ってｗ、ブラウン運動を取り上げました。その際にランダムウォークモデルを使って数値計算実験をして、拡散現象と正規分布の関係すなわちアインシュタインの関係式について論じました。要約しますと、ブラウン運動の実験事実の一つ「平均ニ乗変位σ^2が時間tに比例する」がアインシュタインの関係式であり、ブラウン粒子が原点からどれくらい離れたかという目安であるσは√tに比例します。

それでは、前回のヌクレオチド頻度を計算するスクリプトの以下の箇所を見てください。

foreach $base (@DNA){
if ($base eq 'a'){
++$count_a;
}elsif ( $base eq 't' ){
++$count_t;
}elsif ( $base eq 'g' ){
++$count_g;
}elsif ( $base eq 'c' ){
++$count_c;
}else{
print "Error! Can't recognize: $base\n";
++$errors;
}
}

という箇所です。これはDNA配列データを一文字ずつ読み、各塩基ごとにその数をカウントしています。このカウント作業を少しいじって、次の測定量を計算します。

一次元の仮想的な粒子を考えて、これをはじめ原点におく。DNA配列データを順に読み込み、ピリミジン (塩基T or C)が出現したら粒子を+1だけ動かし、プリン (塩基A or G)が出現したら粒子を-1だけ動かす。この操作をDNA配列の端から端まで行い、一次元の仮想粒子の軌道としてマッピングする。

すなわちゲノムのポジションを各時刻tとみなしたときに、DNA配列の並びによって仮想粒子の位置x(t)を順次決めていくという作業です。この手法はDNA walkとも呼ばれ、いかに各塩基組成 (この場合A+G or T+C)が局所的に変化するかを調べることに相当しています。DNA walkを計算するにはさきほどのカウントサブルーチンに手を加えることで実現できます。


1 #!/usr/bin/perl -w
2 # DNA walk
3
4 use strict;
5 use warnings;
6
7 my $DNA = '';
8 $DNA = read_gbk();
9 DNA_walk($DNA);
10
11 exit;
12
13 ##########################
14 # gbkファイルからDNAの読み込み
15 ##########################
16 sub read_gbk{
17 my $str;
18 my $DNA = '';
19 while(<>){ # 標準入力からデータ読み込む
20 if($_ =~ /^ORIGIN/){
21 while(<>){
22 if($_ =~ /^\/\//){last;} # "//"がファイルの最後
23 $str = $_;
24 $str =~ s/[0-9]//g; # 数字は取り除く
25 $str =~ s/\s//g; # 空白文字を取り除く
26 $DNA = $DNA.$str; # $strを1つの文字列$DNAへと連結する
27 }
28 }
29 }
30
31 return $DNA
32 }
33
34 ##########################
35 # DNA walk の計算
36 ##########################
37 sub DNA_walk{
38
39 my $t = 0;
40 my $x = 0;
41
42 my $errors = 0; # 塩基以外の文字エラーの数
43 my $base;
44
45 my ($DNA) = @_;
46 my @DNA = split( '', $DNA);
47
48 foreach $base (@DNA){
49 if ($base eq 'a'){
50 --$x;
51 }elsif ( $base eq 't' ){
52 ++$x;
53 }elsif ( $base eq 'g' ){
54 --$x;
55 }elsif ( $base eq 'c' ){
56 ++$x;
57 }else{
58 print "Error! Can't recognize: $base\n";
59 ++$errors;
60 }
61 printf("%d\t%d\n", $t, $x);
62 $t++;
63 }
64 }


実行結果
$ perl DNA_walk.pl < HUMHBB.gbk > HUMHBB.dw
$ head HUMHBB.dw
0 -1
1 -2
2 -3
3 -2
4 -1
5 0
6 1
7 0
8 -1
9 0
...
...
...

これを一次元のランダムウォークt-x(t)図になぞらえてプロットしますと、以下のように
なります。



おお、なんじゃこりゃ！！なんだか意味深なグラフが得られました。GC-skewのように塩基組成のバイアスをよく表していそうです。なぜなら仮想粒子は、塩基A or Gが固まって出現したときには一次元x軸上でプラス方向へ動き続けるからです。反対に、T or Cが連続して出現すればマイナス方向へ動き続け、結果として局所的な塩基組成の変化がこのt-x(t)グラフから読み取れることになります。

また、上のスクリプトでは A or Gという塩基に着目していますが、これも任意に変えることが
できます。そのマッピングルールとしては以下のようなものも考えられます。

---------------
SW C or G
RY A or G
KM G or T
---------------

蛇足ですが、DNA配列などのデータにはヌクレオチドのシンボルとして
以下のようなコード表が用いられていることが多いです。(詳細)
----------------------------
A アデニン (Adenine)
T チミン (Thymine)
G グアニン (Guanine)
C シトシン (Cytosine)
U ウラシル (Uracil)
M A or C
R A or G
W A or T
S C or G
Y C or T
K G or T
N 任意のヌクレオチド
----------------------------




References

[1] バイオインフォマティクスのためのPerl入門
[2] J.R. Lobry (1999) "Genomic Landscapes", Microbiology Today 26: 164-165.
http://www.socgenmicrobiol.org.uk/pubs/micro_today/pdf/049906.pdf

DNA配列の統計的性質 第6回 〜DNA配列データ処理(その2)

前回までの過去記事

前回、簡単なPerlスクリプトを用いてGenBank形式ファイル中のDNA配列データを抜き出しました。さて今回は、このスクリプトを雛型にして統計解析をやっていきたいと思います。DNA配列データを目の前にして統計解析をやるぞ！と意気込んだところで、どのように何をやるかは非常にセンスを問われるところです。むしろ何か知りたいことがあって、はじめて統計的解析の意味があるといって間違いありません。うーむ、何か知りたいこと、、、


これといってありません


とは口が裂けても言えませんので、ひとまずDNA中のヌクレオチド頻度 (塩基A, T, G, Cの数)を計算してみます。これは前回gbkファイルから抜き出したDNA配列データを一文字ずつ読んで、Aが何個、Tが何個、、、と数えてやればよいです。またこのヌクレオチド頻度から平均のGC含有量を求めてみます。


1 #!/usr/bin/perl -w
2 # calc nucleotide frequency
3
4 use strict;
5 use warnings;
6
7 my $DNA = '';
8 $DNA = read_gbk();
9 freq($DNA);
10
11 exit;
12
13 ##########################
14 # gbkファイルからDNAの読み込み
15 ##########################
16 sub read_gbk{
17 my $str;
18 my $DNA = '';
19 while(<>){ # 標準入力からデータ読み込む
20 if($_ =~ /^ORIGIN/){
21 while(<>){
22 if($_ =~ /^\/\//){last;} # "//"がファイルの最後
23 $str = $_;
24 $str =~ s/[0-9]//g; # 数字は取り除く
25 $str =~ s/\s//g; # 空白文字を取り除く
26 $DNA = $DNA.$str; # $strを1つの文字列$DNAへと連結する
27 }
28 }
29 }
30
31 return $DNA
32 }
33
34 ##########################
35 # frequency の計算
36 ##########################
37 sub freq{
38
39 my $count_a = 0; # 塩基aの数
40 my $count_t = 0;
41 my $count_g = 0;
42 my $count_c = 0;
43 my $errors = 0; # 塩基以外の文字エラーの数
44 my $base; my $total;
45 my $gc_content = 0.0;
46 my $at_content = 0.0;
47
48 my ($DNA) = @_;
49 my @DNA = split( '', $DNA);
50
51 foreach $base (@DNA){
52 if ($base eq 'a'){
53 ++$count_a;
54 }elsif ( $base eq 't' ){
55 ++$count_t;
56 }elsif ( $base eq 'g' ){
57 ++$count_g;
58 }elsif ( $base eq 'c' ){
59 ++$count_c;
60 }else{
61 print "Error! Can't recognize: $base\n";
62 ++$errors;
63 }
64 }
65
66 # calc total
67 $total = $count_a + $count_t + $count_g + $count_c;
68
69 # calc GC% and AT%
70 $gc_content = ($count_g + $count_c) / $total * 100;
71 $at_content = ($count_a + $count_t) / $total * 100;
72
73 # output
74 print "a: ", $count_a, " t: ", $count_t," g: ", $count_g," c: ", $count_c, " errors: ", $errors, "\n";
75 print "total: $total bp\n";
76 printf("ave GC percent %.1f\n", $gc_content);
77 printf("ave AT percent %.1f\n", $at_content);
78
79 }


出力結果

$ perl nucleotide_freq.pl < HUMHBB.gbk
a: 22068 t: 22309 g: 14785 c: 14146 errors: 0
total: 73308 bp
ave GC percent 39.5
ave AT percent 60.5


非常にシンプルな計算結果ですね。クールですね。塩基a, t, g, cの数と全塩基数およびGC含有量・AT含有量が出力されています。一応、ちゃんと読めているかどうか塩基のトータル数を確認してみます。コマンドheadを叩いて、、、


LOCUS HUMHBB 73308 bp DNA linear PRI 13-DEC-2004

大丈夫そうです。

以下、余談になります。このスクリプトではGC%を求めていますが、あくまで入力したDNA配列の「平均」のGC%を求めていることになります。つまりこの量だけでは、配列全体でGC%が高いか低いかだけをいうことができます。しかしながら、予想されますようにゲノム上ではGC-richな箇所、-poorな箇所があります。そこで、このスクリプトにもう一工夫加えますと、配列に沿って局所的なGC%を求めることができます。あるウィンドウ幅を設けて、ゲノムDNA配列をスキャンしていき (移動ステップを定めて)、そのウィンドウ内のGC%を求めるということです。なお、このスクリプトは課題として残しておきます。

さてHUMHBB.gbkでウィンドウ幅1000 bpのステップ幅500 bpとしますと以下のようなグラフが描けます。

拡大図

これでもはや、大腸菌ゲノムであろうが、インフルエンザ菌ゲノムであろうがドンと来い！という感じですね。ヌクレオチドの頻度はわかるし、ゲノム GC%もわかります。さきほど計算しましたとおり、HUMHBBの配列は78 kb程度の長さです。現在利用可能なコンプリートゲノムなどを用いて、このようなGC%の分布を比較ゲノム的に見てみるのも興味深いと思います。特に原核生物と高等真核生物ゲノムの違いや遺伝子領域との関連、生物系統との関係など。またヒトゲノムにおいてはS期の複製タイミングとの関連も指摘されています。

おっと！勘違いされては困りますので、きついツッコミがある前に言っておきますが、「DNAの統計的性質」はこれで終りではありません！！当たり前ですね。このまま終わったら、詐欺ですね(注1)。あくまで徐々に段階を追っていってるだけです。誰がなんと言おうともネタ切れだとか、夜は眠くて作業できないよーっていう言い訳をしたいわけではありません。

次回ですが、しばらくオタマジャクシ本のパクリみたいな話しの流れでいこうかと思います。そんで、ブラウン運動のことをすっかり忘れた頃、新たなトリビアが生まれるはずです・・・。なんちって！！ｗ

#スクリプトに関してですが、このブログ上ではインデントをまったく行ってない風になっていますが、ホントはちゃんとしているんですよ！！空白が自動で詰められたのを直すのが面倒なのですわ。

すっかり余談の方が長くなってしまいました。(汗)



(注1) じつは以前のブログで、ホントにこの記事を最後に活動が停止してしまいました。心待ちにしていた方がおられたとしたら、この場を借りてお詫びいたします。

References

[1] バイオインフォマティクスのためのPerl入門
[2] Woodfine K, et al. (2004) "Replication timing of the human genome", Hum Mol Genet. 13:191-202.
[3] Complete Genomes in KEGG

DNA配列の統計的性質 第5回 〜DNA配列データ処理(その1)

今回からDNA塩基配列の統計的性質を調べる段階に入ることになります。そもそも表題にある通り、DNA配列を扱うのがメインですから。前回までは、Perlの基礎のおさらいと、この後に必要な背景としての確率的シミュレーション(ブラウン運動)を取り扱いました。アインシュタインが確立したブラウン運動の理論を用いてのペランらの実験によって原子の存在証明がなされたという奇跡のような偉業の雰囲気もちょっとだけ掴めたのではないでしょうか。(えっ、掴めてません？・・・大丈夫。世の中わかったような気になった人間の勝ちです、たぶん)

ちなみに、なぜPerlを扱うのか？とかC/C++じゃダメなのか？といった質問が各界から多数寄せられていますがｗ、多分に個人的趣味です！とだけお答えしておきます。なお、前回までのことは、私も激推しする統計Rというソフトならさらに簡単にできてしまいます。ハッハッハ。

さて、早速バリバリと進みたいところですが、ちょっと待ってください。まだ手元に解析対象であるDNA配列データがありません。ええ、わかっています。なので、データとってくるところから説明したいと思います。いやぁ、なんて親切なんでしょう。(別にネタに困って、お茶を濁してるわけではありませんYO！)

とりあえず、NCBIのサイトに行ってみましょう。ポチッとな。ｗ
TOPページにたどり着きましたら、Search欄からNucleotideを選び、HUMHBBというキーワードで検索してみます。



すると、しょっぱいいっぱい検索結果が返ってきます。この中からU01317というACCESSIONをクリックします。



DisplayをGenBankとしてGet Sequenceを押してファイルに保存します。1〜2コの配列をゲットするならばこれで良いです。大量にとってくる場合は、それなりのデータベースの扱い方をせねばなりませんが、ここでは割愛させていただきます。HUMHBBとは、赤血球中にある酸素を運搬するタンパク質であるヘモグロビンを構成するβグロビンをコードする遺伝子のことです。ちなみにHUMとなってますので、ヒトのβグロビンを表します。

このゲットしたGenBank形式のファイルは、前半部分に配列データのfeature(LOCUSとかACCESSIONや、参照すべき論文、遺伝子の位置などのアノテーション)の記載があり、後半にDNA配列が載せてあります。

目下、必要な情報は後半のDNA配列部分ですので、これを抜き出す操作をせねばなりません。これはPerlを使えば、簡単に(?)行うことができます。別に手で書き写したって構わないですよ、わたしゃ。←オイ！

ファイル名: gb2seq.pl
1 #!/usr/bin/perl -w
2 #GenBank file --> DNA seq
3
4 use strict;
5 use warnings;
6
7 my $str;
8
9 while(<>){ # 標準入力からデータ読み込む
10 if($_ =~ /^ORIGIN/){
11 while(<>){
12 if($_ =~ /^\/\//){last;} # "//"がファイルの最後
13 $str = $_;
14 $str =~ s/[0-9]//g; # 数字は取り除く
15 $str =~ s/\s//g; # 空白文字を取り除く
16 print "$str\n";
17 }
18 }
19 }
20
21 exit;

かなり省略した書き方をしてしまったので、可読性が低くてスンマセン。"< >"とか"$_"に注意してください。GenBank形式の場合、DNA配列データは必ず"ORIGIN"という行の直後、行番号とともに記載され、"//"行で終わります。したがって、"ORIGIN"行が現われるまではスルーして、DNA塩基を表す文字(a, t, g, c)だけを読み込みます。

これを下記のように実行します。これはGenBank形式のファイルを標準入力からリダイレクトしてやって、ファイルに出力しなくてはなりません。(あくまでUNIX系OSでの実行を想定しています)

$ perl gb2seq.pl < HUMHBB.gbk > HUMHBB.seq
$ less HUMHBB.seq
gaattctaatctccctctcaaccctacagtcacccatttggtatattaaagatgtgttgt
ctactgtctagtatccctcaagtagtgtcaggaattagtcatttaaatagtctgcaagcc
aggagtggtggctcatgtctgtaattccagcactggagaggtagaagtgggaggactgct
tgagctcaagagtttgatattatcctggacaacatagcaagacctcgtctctacttaaaa
aaaaaaaaattagccaggcatgtgatgtacacctgtagtcccagctactcaggaggccga
aatgggaggatcccttgagctcaggaggtcaaggctgcagtgagacatgatcttgccact
gcactccagcctggacagcagagtgaaaccttgcctcacgaaacagaatacaaaaacaaa
caaacaaaaaactgctccgcaatgcgcttccttgatgctctaccacataggtctgggtac
tttgtacacattatctcattgctgttcgtaattgttagattaattttgtaatattgatat
...

というように、GenBank形式のファイルにあったDNA配列データだけが抜き出せました。数字やスペースもちゃんと取り除かれています。このように、Perlではわずかな行のスクリプトで、テキストファイルの処理できることがウリなわけですが、最初に申しました通り、C/C++でやろうがJavaでやろうが、それは個人的趣味ということになります。

次回は、このスクリプトを雛型にして(今度こそ)統計処理を少しずつやっていきたいと思います。


References

[1] バイオインフォマティクスのためのPerl入門

DNA配列の統計的性質 第4回 〜ブラウン粒子と拡散(その2)

統計学の教科書に必ず出てくる正規分布に関して、ブラウン運動を通し違った面から理解が深まると思います。来年2005年はアインシュタインの奇跡の年(1905年)から100年を記念して国際物理年となっていることはご存知でしょうか。アインシュタインの偉業のひとつブラウン運動に関する仕事の一端が、この記事を通して少しでもご理解いただければ幸いです。(Dec 17, 2004)


前回は、多数のブラウン粒子がいかに振る舞うのかという問題を扱いました。そのとき、個々の粒子一つ一つの動きはランダムですが、その集団としての振る舞いは確率的・統計的に扱えると述べました。結局、ブラウン粒子が正規分布の形で広がることもシミュレーションしてわかったと思います(厳密に証明したわけではないですが)。

ここでは、正規分布の形で広がるということを前提に話しを進めます。当然、確率・統計学の教科書に出てくるような正規分布の式が活用されるべきところですが、数式をここに打ち出すのがメンドイので、省略させていただいて。まぁ、正規分布の式をアタマに思い浮かべるなり、教科書で確認するなりしていただければと思います。

正規分布を規定するのに、平均と標準偏差(or 分散)という変量があったはずです。ブラウン粒子の場合には、一次元x軸上で原点x = 0から出発して各ステップは確率1/2で左右に移動するわけですから、平均は0です。また拡散のシミュレーションで見ました通り、時間tとともに分布は広がっていきます。この分布の広がり具合を表すのが標準偏差(or 分散)でした。すなわち、時間の経過とともに標準偏差(or 分散)は増えていきます。

標準偏差を慣例に従いσで表しますと、正規分布ですので、ブラウン粒子はσ、2σ、3σの距離内に、それぞれ確率68.3%、95.5%、99.7%で見出すことができることになるわけです(この展開、いいですよね？)。そこで、ある時間経過したときにブラウン粒子が平均してどれくらい広がっているのかが標準偏差(or 分散)という量を測ることで知ることができます。

このことを調べるには、x^2(xの2乗)の平均値を計算して、平方根をとります(分散の平方根が標準偏差であったことを思い出してください)。x^2の平均は平均ニ乗変位とも呼ばれる量でσ^2です。もう一度言いますが、このσはブラウン粒子の原点からの距離の平均の目安を示します。したがって、時間tの間にブラウン粒子がどれだけ広がったかを調べるプログラムは以下のようになります。

1 #!/usr/bin/perl -w
2 #random walk 1D multi mean square
3
4 use strict;
5 use warnings;
6
7
8 #カウンタ変数、フラグ
9 my $i; #粒子のステップ数
10 my $j; #粒子数
11 my $t; #経過時間t
12 my $flag;
13
14 #ブラウン粒子の座標
15 my $x;
16 my $xnew;
17
18 #$xを格納する配列
19 my @distance;
20 @distance=();
21
22 #変数初期化
23 $x = $xnew = 0;
24
25 srand(time|$$);
26
27 for ($t = 1; $t < 7000; $t = $t * 2) {
28 for ($j = 0; $j < 1000; $j++) {
29 for ($i = 0; $i < $t; $i++) {
30 $flag = int(dice2());
31 if ($flag == 0){
32 $xnew = $x + 1;
33 }elsif ($flag == 1) {
34 $xnew = $x - 1;
35 }
36 $x = $xnew;
37 }
38 push (@distance, $x);
39 $x = $xnew = 0;
40 }
41 print $i, " ", mean_square(@distance), "\n";
42 @distance = ();
43 }
44
45 exit;
46
47
48 ###################################
49 # dice2
50 # call すると確率1/2で0 or 1を返す
51 ###################################
52 sub dice2{
53 my $i;
54 my $flag;
55 $flag = int(rand() * 2);
56 if ($flag == 1) {
57 $i = 1;
58 }else {
59 $i = 0;
60 }
61 return $i;
62 }
63
64 ###################################
65 # mean square (平均ニ乗変位)
66 # call すると平均ニ乗変位を返す
67 ###################################
68 sub mean_square{
69 my $i;
70 my $sum; #距離の合計
71 my $sum_sq; #距離の合計のニ乗
72 my $mean; #距離の平均
73 my $mean_sq; #ニ乗変位
74 my $n;
75
76 foreach $x (@_) {
77 $sum += $x;
78 $sum_sq += $x ** 2;
79 }
80
81 $n = @_; #配列の要素数
82 return (0, 0) if $n == 0; #データ無
83 return ($sum, 0) if $n == 1; #データ1コ
84
85
86 $mean = $sum / $n;
87 $mean_sq = ($sum_sq - $mean**2) / ($n - 1);
88
89 return $mean_sq;
90 }

カウンタ変数の意味について注意してください(27-43行目)。粒子を動かす経過時間の$tは最終的にグラフ化するときに両対数にしますので、インクリメントを指数関数的にしています。なお、標準偏差σではなく分散σ^2を求めています。87行目で分散を求める際に、個数n-1で割っていることにも注意してください。それから、解りやすさを追求したため、キレイなソースではないことにも注意してください(というか現時点での私の実力です)。(^^;

さて、この実行結果を図示したものが以下です。



結局、これは平均ニ乗変位σ^2を時間tの関数として描いたことになります(log-log plot)。直線は傾きが1であることを表しています。図を見ますと、明らかにσ^2は時間tに比例しています(一次関数)。したがって、ブラウン粒子の広がり具合、すなわち粒子が原点からどれくらい離れたかという目安であるσは√tに比例します。じつはこれ、ブラウン運動の実験事実の一つなのです。「平均ニ乗変位σ^2が時間tに比例する」は、アインシュタインの関係式とも呼ばれています。なお、実際にペランの弟子が顕微鏡を覗いて乳香粒子の位置を30秒ごとに測定し、出発点からの距離が√tの関数になっていることを確かめました。こうして原子の実在性が証明され、厳密なアヴォガドロ数の算定へと繋がります(大雑把でスンマセン！！)。

一応、今回までで準備運動(！)は終りで、次回からDNA配列データを使っていろいろ調べたりしてみたいと思います。じつは今日扱った事柄は極めて重要なことです。それが徐々に明らかになっていくはずです(なっていって欲しいです！！)。

DNAはどこに消えた？

当たり前過ぎて何も考えてなかったのですが、DNAの統計的性質を調べるに当たって、その肝心なDNA配列データはどこにあるの？という疑問もあるかと思います。ググっちゃえば(google検索すれば)、すぐに腐るほど見つかると思うのですが、基本ということで超超超超有名データベースを紹介しておきます。バイオインフォマティクスのテキストとか見ればそれこそクドイほどの説明があるのですが、軽く紹介しておきます。このブログでも次回以降使うと思いますし。(いや、次々回以降になる可能性大です)

National Center for Biotechnology Information

日本DNAデータバンク (DDBJ)


こういったデータベースにアクセスしてみて、何等かの遺伝子やらのDNA配列データとか眺めるのもオツなものですね！？なお、これだけ簡単に手に入ってしまう遺伝情報ですが、そのシーケンシング(sequencing: DNA配列を実験的に読み取ること)には多大なる苦労があったろうと思います。そのことに感謝して、使用しましょう。いや、個人的にはホントその苦労が全然わからないんですけどネ。

あ。あと、ちなみに関連事項として本日はこんな新聞記事が掲載されました！！

ヒト遺伝子のカタログ、世界に無償公開　国立遺伝研など（朝日新聞）

ヒトの体内で実際に働くことが実験的に確認された遺伝子を載せたデータベースが、日本を中心にした１２カ国の協力で完成し、１６日から世界に無償で公開される。遺伝子とその構造や機能情報などを結びつけたヒト遺伝子の詳細なカタログだ。病気と遺伝子の関連、薬の開発、生物学基礎研究など広く役立ちそうだ。(asahi.com: アサヒ・コムより)

H-Invitational データベース
H-Invitational Database CIB-DDBJ Flat File Server

DNA配列の統計的性質 第3回 〜ブラウン粒子と拡散(その1)

前回はPerlのスクリプトを用いて乱数列を生成し、ブラウン運動をシミュレートしてみました。今回はもう少し現実的にブラウン運動および拡散を扱ってみたいと思います。

コップの中に入った水に赤いインクを適量たらしてみたときの現象を想像します。水にもインクにも膨大な数の分子が存在していることは学校で習ったはずです。となりますと、当然ブラウン運動する粒子も多数でなければ、この現象を理解することができません。すなわち、互いに独立にあちらこちらへと動き回るブラウン粒子団全体の振る舞いを見なければならないのです。といいますか、微小な粒子の一つ一つの動きはランダムなのですが、その集団の振る舞いは規則的な扱いができるというのが統計物理のセントラルドグマです。ｗ　←大雑把ですいません

前回の一次元ランダムウォークのところで述べましたが、あの計算機実験はあくまで、ある一つの粒子に対するものでした。したがいまして、Perlのスクリプトにもう一つループを加えて1000個 (どこが現実的なんだ？)のブラウン粒子を発生させてみましょう。このとき、各ブラウン粒子の最終到達位置だけを記録することにします。コードは単純に以下のようになります。

1 #!/usr/bin/perl -w
2 #random walk 1D
3
4 use strict;
5 use warnings;
6
7
8 #カウンタ変数、フラグ
9 my $i; #粒子のステップ数
10 my $j; #粒子数
11 my $flag;
12
13 #ブラウン粒子の座標
14 my $x;
15 my $xnew;
16
17 #変数初期化
18 $x = $xnew = 0;
19
20 srand(time|$$);
21
22 for ($j = 0; $j < 1000; $j++) {
23 for ($i = 0; $i < 5000; $i++) {
24 $flag = int(dice2());
25 if ($flag == 0){
26 $xnew = $x + 1;
27 }elsif ($flag == 1) {
28 $xnew = $x - 1;
29 }
30 $x = $xnew;
31 }
32 print $x, "\n";
33 $x = $xnew = 0;
34 }
35
36 exit;
37
38
39 ###################################
40 # dice2
41 # call すると確率1/2で0 or 1を返す
42 ###################################
43 sub dice2{
44 my $i;
45 my $flag;
46 $flag = int(rand() * 2);
47 if ($flag == 1) {
48 $i = 1;
49 }else {
50 $i = 0;
51 }
52 return $i;
53 }
54

ループカウンタ変数の$i, $jが何を表しているかについて注意してください。さて、これを実行しますと、1000個分のブラウン粒子の最終到達位置が出力されます。この出力データをもとにExcelなり統計ソフトR [1]なりでヒストグラムを書いてみます。このヒストグラムは、横軸が最終到達位置となり、縦軸は各々の位置に至ったブラウン粒子の数を表します(下図)。ただし、縦軸は規格化しているので(積分値が1)、正確には原点から距離x離れた位置にブラウン粒子を見出す確率密度を表します。


各粒子が1000ステップ後のヒストグラム


各粒子が5000ステップ後のヒストグラム


各粒子が10000ステップ後のヒストグラム


あ、拡大して見て下さいね。ステップが進むごとに粒子集団が拡散しているのがわかると思います。さて、この分布はどこかで見たことがありませんか。そうです、実験の測定なんかで独立な誤差が積もり積もったときに現れるガウス分布(正規分布)なのです。この粒子系の実験を何回もやり、ヒストグラムの幅を無限に小さくとれば、ガウス分布に近づくことがわかります。このように、time ステップtを変えつつ、このヒストグラムを順次書いてみますと、これは1000個の粒子が初期状態t=0, x=0から拡散していく様子を表すことになります。つまり、ほとんどのブラウン粒子はその到達距離が平均のx=0のまわりに分布し、tが経過するほどブラウン粒子のいくつかは平均からそこそこ離れた位置まで到達することもあることがわかります。ここで確率・統計学で出てくる分布について思い出しますと、分布の幅つまり標準偏差(or分散)ですが、このモーメント量はブラウン粒子でいえば、まさにブラウン粒子系の到達距離に相当しているのです。

次回は、この時間とともに広がる粒子群の振る舞いを定量づける分散について議論します。



References

フリーの統計環境Rは、非常にオススメです[1]。日々新しい機能が追加されてますし、英語はもちろんのこと日本語訳のマニュアルもありますし、様々な分野の方々のご尽力によってtipsが書かれています。その中のひとつが[2]。よくできてます。スバラシイです！！

[1] The R Project for Statistical Computing
[2] 統計解析用フリーソフト・R-Tips

DNA配列の統計的性質 第2回 〜擬似乱数の適用例 (すぃすぃすぃーっと酔歩)

タイトルにあります「DNAの統計的性質」を考察するところまでなかなか到達できないでいますが(^^;、基本は牛歩で参りたいと思います。さて、今回はまず準備として以下のスクリプトを叩き台にします。これは出目が2つ (0か1) しかないサイコロです。前回の擬似乱数発生器を少しもじったものです。

1 #!/usr/bin/perl -w
2 #dice2
3
4 use strict;
5 use warnings;
6
7 my $i;
8
9 #擬似乱数のseed
10 srand(time|$$); #
11
12 for ($i = 0; $i < 10; $i++) {
13 print dice2(), "\n";
14 }
15
16 exit;
17
18 ###################################
19 # dice2
20 # call すると確率1/2で0 or 1を返す
21 ###################################
22 sub dice2{
23 my $i;
24 my $flag;
25 $flag = int(rand() * 2);
26 if ($flag == 1) {
27 $i = 1;
28 }else {
29 $i = 0;
30 }
31 return $i;
32 }

当然ですが、プログラミングはPerlに限らず機能ごとに関数を定義して利用するのが上等手段です。つまり、ひとまとめにできる処理はサブルーチンにまとめるってことなのです。22〜32行目のサブルーチンdice2は引数の受け渡しがありません。呼び出すと確率1/2で0か1を返します。ちょっとした工夫が25行目の処理です。rand()から返ってくる値は0以上1未満の実数ですので、これを2倍してから整数型に変換しますと、ちょうど0.5を境に0か1になります。

dice2.plの実行例
$ perl dice2.pl
0
1
0
1
1
1
1
0
0
0

これにしたがって、出目が4つのサイコロはすぐに作ることができますね。さて、出目が4のサイコロdice4ができたとして、それを利用した例を作ってみます。乱数をもっともダイレクトに応用した例がブラウン運動の計算機実験となります。(なんてスバラシイ展開なんだ！ｗ)

1827年、イギリスの植物学者ブラウン(Robert Brown)は、水中に入れた花粉から出た微粒子がフラフラと不規則に動くことを顕微鏡で観察しました。この微粒子の運動はブラウン運動とよばれています。ブラウン運動は、かのアインシュタインが数学的定式化を行い、その結果として後に原子・分子の存在を証明することとなりました(詳しくは文献[1]を参照)。このブラウン運動を離散化した数学モデルをランダムウォーク (酔歩)と呼びます。

ここでは最初に簡単のため一次元x軸上の粒子を考えます。のちに2次元に拡張するのはたやすいです。さて、その粒子が1 timeステップごとに1/2の確率で左右のどちらかに一歩進むという(1次元ランダムウォーク)シミュレーションをします。そのスクリプトは以下のようになります。サブルーチンのdice2()がすでにあるので簡単ですね。

1 #!/usr/bin/perl -w
2 #random walk 1D
3
4 use strict;
5 use warnings;
6
7
8 #カウンタ変数、フラグ
9 my $i;
10 my $flag;
11
12 #ブラウン粒子の座標
13 my $x;
14 my $xnew;
15
16 #変数初期化
17 $x = $xnew = 0;
18
19 srand(time|$$);
20
21 print $x, "\n";
22 for ($i = 0; $i < 10000; $i++) {
23 $flag = int(dice2());
24 if ($flag == 0){
25 $xnew = $x + 1;
26 }elsif ($flag == 1) {
27 $xnew = $x - 1;
28 }
29 $x = $xnew;
30 print $x, "\n";
31 }
32
33 exit;
34
35
36 ###################################
37 # dice2
38 # call すると確率1/2で0 or 1を返す
39 ###################################
40 sub dice2{
41 my $i;
42 my $flag;
43 $flag = int(rand() * 2);
44 if ($flag == 1) {
45 $i = 1;
46 }else {
47 $i = 0;
48 }
49 return $i;
50 }


22行目のfor文によって粒子を何ステップ動かすかを指定しています(例: 10000ステップ)。結果をファイルにリダイレクトし、gnuplotのようなグラフ描画ソフトで時間tに対する粒子の位置x(t)をプロットしたのが下図です。ただし、時刻t = 0において粒子の初期位置はx(0) = 0とします。


これは1次元のランダムウォークの時系列t-x(t)図 (10000ステップ)、横軸はt、縦軸がx(t)です。

このように一つの粒子は、不規則に上下をして10000ステップの後、ある位置に到達します。このような粒子がたくさん(たとえば、〜10^23個)あるときに各粒子の到達点はどうなるでしょうか (次回以降徐々に説明したいと思います)。じつはこのランダムウォーク、拡散現象に密接に関連しています。数学的にも厳密に理解がなされていて、経済学における金融理論なんかにも応用されて、ランダムな現象を理解する上で非常に強力な理論なのです。

次に2次元に拡張して、粒子の位置を(x,y)としてx-y図をプロットしたものが下図です。これはフラフラと動いた1個の粒子の軌跡を現しています。先ほど、拡散に関係していると言いましたが、だんだんそれっぽくなってきました。まさに水の入ったコップに赤いインクを落としたりしたときのアレですよ。




以上、計算機で作り出した乱数列を用いて、ブラウン運動のシミュレーションをしてみました。2次元の場合、ブラウン粒子は各ステップで上下左右のいずれかにジャンプしたわけですが、その移動位置の決定は粒子が以前どこに居たかによらず (数学的にいうと、互いに独立な事象といいます) 行われていたことに注意してください。また、悪い乱数 (出現に相関がある場合)を使うと、このような確率的シミュレーションでは非常に具合が悪いことも付け加えておきます。

今回はまったくbioinformaticsとは思えない内容でしたが、少しの間の辛抱です。といいますか、結局はどこまでいっても一般に描くようなbioinformaticsの範疇に入らないかもしれません。ここまで書いてきて、だんだん不安になってきました。。。


References
[1] 米沢 富美子 著, ブラウン運動 物理学One Point

DNA配列の統計的性質 第1回 〜はじめに〜

以下続く、DNA配列の統計的性質ネタは昔書き溜めたものでして、若干古臭い感じを受けなくもないですが、それなりに味わいのある点もあると思いますし、誰かの役に立つはずだ！とポジティブに考え公開します。Perlを勉強しつつ進行して行ったものなので、初心者の方にはある意味勇気を与える内容だと自負しています。ｗ


今回から数回にわたってperlを用いたDNAのモデル化について連載をしたいと思います。モデルといいましても立体構造がこうなって・・・という話ではありません。ゲノムプロジェクトで得られたテキストデータとしてのDNA配列に対する数理的考察を行います。

はじめに、モデル化に必要なperlの基礎事項を確認します。後々それをベースに改良を加えていき最終的な考察に必要なプログラムとなるように導きます。本来は、背景となる研究分野のreviewを最初に行うべきかもしれません。でもそんなの面倒くさいし、背景を知らないままperlのスクリプトを動かして、(テキストデータですが)DNA配列を扱って遊んでみるのも一興と考えました。なお、レベル的にはバイオインフォマティクスのためのPerl入門の範囲内でいけるはずです。というより、この本をパラパラとめくっていたら今回ネタを思いつきました。かなりマニアックな内容で狭〜いところ突いていきますので、背景が必要なときには随時説明をすることにします(私自身のためにも)。

ちなみに原書のページに行きますと、原書本文中に掲載されているperlサンプルコードがダウンロードできます。

さて、ではperlの復習からやります。擬似乱数列を発生させるスクリプトをいじっていじってイジり倒しながらテクニカルな確認をしてみます (とてもじゃありませんが、Perlのすべての文法を網羅できませんので、適宜Perl解説本を参照してください)。また実行例はUNIX系OSを想定していますので、$はシェルプロンプトを表します。

乱数列とは文字通り数が不規則に並んだ列であり、その出現順に相関があってはならないものです。プログラミングで使用するような計算機で生成する乱数 (擬似乱数と呼びます) は、厳密な意味で乱数とはいえません。そして考えうる統計的検定のすべてをクリアする擬似乱数は存在しないのです、、、などと、これ以上小難しいことを
考えず、デタラメっぽいニセの乱数を擬似乱数生成アルゴリズムで作るということです。

まずは次のスクリプト (ファイル名: random0.pl) を見てみます。


1 #!/usr/bin/perl -w # -wは警告表示
2 #random number generator
3
4 use strict;　　　　　 　　　 # おまじない
5 use warnings; # おまじない2
6
7 my $i; # カウンタループ変数$iの宣言
8
9
10 srand(time|$$); # 擬似乱数のseed
11
12
13 for ($i = 0; $i < 10; ++$i) { # ++は,$iを１増やす演算子
14 print rand(), "\n"; # "\n" は改行．
15 }
16
17 exit; # 終了の明示


これは、10個の乱数を生成するプログラムです。非常に単純ですね！ｗ
"#"はコメントです。ただし、1行目はコマンドインタプリタ行と呼ばれてまして、このテキストファイルがPerlのプログラムであることを明示しています。さらにいうと、/usr/bin/perlという部分は、Perlの実行ファイルがインストールされて
いるディレクトリです。なので、当然システムによって変わる場合があります。

「おまじない」の部分は、書いておく方が今後良いでしょう。ｗ
次にループカウンタ$iを宣言します。i, j, kといった文字は通常、プログラミングのときに整数型の変数に使用し、特にループカウンタとして好まれます。

さて、このプログラムの肝腎な部分の一つが10行目のsrand()関数です。この関数は後に出てくるrand()関数が乱数を発生させるための種（seed）を作ります。かなり大雑把な説明ですが。乱数発生器というのは、seedが同じであればまったく同じ乱数列を作ってしまいます。カッコ内のtime|$$は計算機の時刻とプロセスIDとの結合をseedとすることを表します。

13〜15行目がforループになっていて、10回分の乱数発生を行います。結果を出力するためのprint文の中で、rand()関数があります。これは先ほどのsrand()とともに乱数列を発生させます。なお、rand() は引数なしで呼び出すと、0以上1未満の数値を返します。

これでやっと実行ができます。以下は実行例です。

$ perl random0.pl
0.518842552629398
0.60127624412732
0.244999595956084
0.736491115155282
0.364612149202987
0.124185651759021
0.973730076706797
0.25492918122935
0.337549689449961
0.693960767996714


10個の乱数列が得られました。次回はこの擬似乱数を使って、簡単な確率論的シミュレーションをしてみたいと思います。説明、結構疲れました。。。ｗ


References
擬似乱数の詳細については、が詳しいです。ウェブサイトではオンラインで様々な数値計算アルゴリズムの解説が読めます(ただし英語ｗ)。C言語ならば日本語版もあります[2]。

[1] Numerical Recipes
[2] ニューメリカルレシピ・イン・シー―C言語による数値計算のレシピ