Rを使ったシミュレーション:「遺伝的浮動」を実感してみよう!

前回はプログラミングの基本について説明しました。もう皆さんは、繰り返しと代入と条件分岐については大体イメージがつかめたと思います。まずはとりあえず、ごく簡単なおさらいをやっておきましょう。

第9, 10回のレポート採点とコメント返却、まだ終わってません

すみません。先週から仕事が立て込んでしまって、まだ採点とコメント返却終わってません。今週中にはなんとかします。

前回のアンケートのうち、難解ポイントに関する意見:

繰り返し・代入・条件分岐のおさらい

次のそれぞれのプログラムの_部分に何を入れるべきかを選択肢から選びなさい。また、Rを使って実行してみましょう。

R editorを使う

Untitled-1.gif

Rでプログラミングを行う場合は、Rについているスクリプトエディター(R editor)を用いると、操作がずーっと楽になります。

Untitled-2.gif

繰り返し

for(i in 1:__){    #{と}の間を10回繰り返す。i の値は毎回1ずつ増える
  print(i)
}

選択肢: a) 1    b)  5     c) 10

代入

goukei=0           #goukeiという変数に初期値0を代入
for(i in 1:10) {   
  goukei=_____+i  #goukeiに前回までのgoukeiの値にiの値を足したもの代入
} 
print(goukei)

選択肢: a) i    b)  goukei   c) print

条件分岐

goukei=0           #goukeiという変数に初期値0を代入
 for(i in 1:10) {
   if(i %% 2 __ 0) { # %% は割り算をした余りを返す演算子。3 %% 2 の計算結果は 1 になる
     goukei=goukei+i # if文の内容は: 「もし i を 2 で割った余りが 0 ならば、{}内を実行
   }               #それ以外は何もしない
 }              #for 文のカッコ閉じるに注意
 print(goukei)

選択肢: a) =    b)  ==     c) <=

この3つの命令それぞれは単純なことなのですが、組み合わせれば、かなり複雑なことでも表現できます。そこで、今回は、上の命令を組み合わせて、シミュレーションに挑戦してみましょう。

遺伝的浮動のシミュレーション

drift.gif

右の図は見たことはありますか?これは、遺伝的浮動のシミュレーションの結果で、集団サイズN=100, 対立遺伝子が2つのとき、片方の対立遺伝子の頻度の変化を表したものです。。最初の世代の対立遺伝子の数は50:50。このグラフを使って、遺伝的浮動とはどういうものか、説明できる人はいますか?また、このグラフはどうやって作ったものですか?

プログラミングを知る前は、右のような図を作るのはとうてい無理と思ってたかもしれません。でも、Rを使って、簡単なプログラムを作れば簡単に描けてしまいます。

「生物学科の皆さんになじみが深い」
「Rをつかってやってみて、結果の見た目が面白い」
「思ったよりはプログラミングが簡単」

という3点がこのテーマを選んだ理由です。

 なお、シミュレーション(simulation, 「シュミレーションでは無いですよ!」)は、数学モデルなどを用いて現実に似た状況をコンピュータ上で再現することです。模擬実験とも呼ばれています。(大量の計算を必要とするので、Rのようなインタプリタ型のソフトでは、大規模なシミュレーションは難しいです)

0. シミュレーションの前段階:モデルの想定

 これからシミュレーションを行うのですが、つぎのような状況を想定してみましょう。受験勉強でやった集団遺伝のお話を、ちょっとだけ思い出してください。

mona2.gif
mona1.gif

まず、対象とする架空の生物(ここでは「モナー」を使いました。(注:モナーについて詳しく知りたい人は、ウィキペディアで「モナー」を検索してみてください。))は2倍体の生物で、雌雄の区別が無く、繁殖期になると集団内の個体が一斉に配偶子を非常にたくさん(数千億個!)ずつ放出した後死亡します。配偶子は生息域の空間内をランダムに漂った後、24時間後に、最も近くにある配偶子と接合し、新しい繁殖個体になります。集団の生息域は空間的に限られているため、親世代と同数の子供個体しか、成熟個体になれません。どの子供が成熟個体になるかは、ランダムに(偶然に)決まるものとします。つまり、このモデルだと、世代間で集団サイズに変動は無いということです。モナーにはトンガリ耳と丸耳の2タイプが対立遺伝子(A, a)で決まることが知られています。トンガリ耳の遺伝子型は(AA, Aa)、丸耳の遺伝子型は(aa)です。

なお、上のモナーの説明は、私が今回のシミュレーションのために、適当に設定したものです。簡単に言うと、「任意交配をする有限サイズの生物集団から遺伝子をランダム抽出して、同サイズの次世代集団を作る」ということです

1. 初めてのシミュレーション:10個体(遺伝子は20個)からなる集団における、遺伝子頻度変化のシミュレーション

個体数=10(遺伝子の集団サイズ=20)、最初の世代における対立遺伝子aの数10、観察世代数=100

 モナー10個体からなる、1つの集団を考えてみましょう。いまこの集団で、対立遺伝子Aの頻度が0.5, 対立遺伝子aの頻度が0.5のとき(10個体だから遺伝子の総数は、20個です。そのうち半分の10個がAで、10個がaです)。集団中の配偶子の接合はランダムに起こると仮定して、100世代の間に、集団中の対立遺伝子の頻度はどのように変化するでしょうか?

slide5.gif

今から実験をするのですから、前もって予想しておいた方が面白そうです
この集団の遺伝子頻度が1か0に固定するまでに、何世代ぐらいかかると思いますか?
0-25  26-50  51-75  76-100  100<

それでは、次のプログラムを走らせて、20個の遺伝子中対立遺伝子aが10個ある任意交配集団の、100世代の間のaの遺伝子頻度の変化をみてみましょう。

プログラムを実行

Rでは変数の名前が日本語でも動作します。今回は皆さんの理解を助けるために、変数の値を日本語で書いてみます。プログラムが見づらくなるので、変数名に英文字を使ったものも示しておきました。

まずは、下のプログラムをコピーして、__のところを埋めて実行してみてください。

集団サイズ = __             #1集団内の遺伝子の総数を集団サイズと呼ぶ。今回は20。
対立遺伝子aの数 = __        #対立遺伝子aの集団中の個数を10個がaであるとする
繰り返し世代数 = ___        #今回は100世代観察することにする
for(i in 1:繰り返し世代数){    #この行から対応する}までの間を100世代分(「繰り返し世代数」の値)繰り返す
 カウンタ=0                      #次世代にランダムに残るaの値を カウンタ という変数に入れる。初期化して0にする。
   for(j in 1:集団サイズ){          #配偶子プールから1個の配偶子を取り出す試行を集団サイズの数繰り返す
     if ( runif(1) < 対立遺伝子aの数/集団サイズ  ){    #乱数を一つ発生させ、
                                                       #それが対立遺伝子aの頻度よりも小さければ
       カウンタ = カウンタ + 1              #aの個数を入れるカウンタという変数の値を1つ増やす
     }
   }
 結果=c(i, 対立遺伝子aの数/集団サイズ)  #何世代目かを表す数字(i)と対立遺伝子aの頻度を結果というオブジェクトに入れる
 print(結果)                      #outの内容を画面に表示する
 対立遺伝子aの数=カウンタ         # 対立遺伝子aの数を新しい世代の対立遺伝子aの個数(カウンタの値)で置き換える
}
質問しますので挙手して答えてください。
100世代目が0だった人、 1だった人、 それ以外の数値だった人
固定した世代数 0-10  10 - 20  20-50  51-75  76-100  100<]

遺伝的浮動のプログラムをステップ・バイ・ステップで理解する

上のプログラム、ここまで説明してきた

繰り返し
代入
条件分岐

だけからできていますね。目新しいものと言えば、

runif()  乱数を発生させる。runif(1)とした場合、0から1までの間の値を1つ発生。

だけです。

では、プログラミングの手順を1つずつ、順を追って見てゆきましょう。この過程さえ理解できれば、シミュレーションのプログラムは自分で書けてしまうのです。

なお、プログラミングの作業自体は、「数学の能力」はさほど重要ではないと思います。
私も「数学」は得意じゃないのです。でもプログラミングはある程度はできます。
プログラミングというのは、石頭のコンピュータでも理解できるように手順を分かりやすく論理的に説明するという点で、
むしろ、「国語能力」とか「説明能力」が必要なのだろうと感じています

1. 「モナーの生活環」に従って、シミュレーションの条件を決める

最初は小さい数で自分の想像が及ぶ範囲の数にしておくと良いと思います。今回は、簡単に、

20個の遺伝子からなる集団
0世代目の対立遺伝子aの数が10(つまり対立遺伝子aの頻度は0.5)
この集団を100世代観察する

という条件にしておきます

2. 上に書いた数字に自分で考えた変数名をつける。(下の3の作業と同時に行われる)

名前を付けておくと、あとでコンピュータに命令するときに便利だからです。名前は、Rの規則で使えないもの以外なら、何でも良い。自分に分かりやすいようにつける。上で決めた条件で数値が決まっているものは代入しておく。

*自分で考えた変数名と定数
対立遺伝子aの数 = 10   #最初の世代の対立遺伝子aの数。次世代以降は1つ前の世代で抽出されたaの数が代入される
カウンタ = 0    #配偶子プールから、次世代の集団に取り出された対立遺伝子aの数。1世代が終わったあとで、
               数字が入るので、初期値として0を入れておく
集団サイズ = 20  #集団の大きさ(つまり、集団に含まれる遺伝子の数)
繰り返し世代数 = 100  #繰り返しの世代数

slide10-1.gif slide10-2.gif

3.遺伝的浮動という現象を「繰り返し、代入、条件分岐」で説明する方法を考える。(上の2の作業と同時に行われる)

紙に絵をかいて、「モナー」の生活環を「繰り返し、代入、条件分岐」で説明できるかを考え、言葉にしてみる。

*自分で考えたプログラムの構造。コマンドを使って書かなくても、
  とりあえず、構造が理解できるように書いてあればよい。

Untitled-5.gif Untitled-4.gif

  1. 何世代分か、以下の操作を繰り返す
  2.  1つの世代の中では、配偶子プールから配偶子を1つ選ぶという操作を20回(遺伝子の数)繰り返す
  3.    取り出した配偶子が対立遺伝子aだったら、変数countaの値を一つ増やす。そうでなかったら何もしない
  4.  集団内で、20回配偶子を取り出す繰り返しが終わったら、対立遺伝子aの数を numaに代入する
  5. 次世代の集団における対立遺伝子aの頻度(num/20)をプリントする

4. プログラムの構造ができたのだから、次ぎに、一つ一つのステップを、コンピュータが理解できる文で表現してみる。

  1. まず、100世代この集団を観察するのだから、繰り返し文で100回繰り返す構造が必要になる。
    *100回繰り返す繰り返し文の構造
    for(i in 0:100) {
          <※ここに繰り返される内容が入る>
    }
  2. 1世代で行われる操作(上の※の部分)は、1つの集団から配偶子を親の集団と同数(つまり、20個)取り出すという繰り返し操作。なので、上の100回繰り返すfor文の中には、さらに、集団中の遺伝子の数だけ繰り返すfor文が入る
    ※集団中の遺伝子数分、つまり、20回繰り返す文
       for(j in 0:20) {
         <※※ここに繰り返される内容が入る>
       }
  3. 上のfor文の中では、配偶子プールか集団サイズの数だけ取り出した配偶子の中に、対立遺伝子aがいくつ含まれているかを数えること。そのとき、
     ※※の操作:
      「配偶子を一つ取り出すとき、遺伝子頻度の確率で、対立遺伝子aをひく。もし、対立遺伝子aを引いたら、
        aの数を入れておく変数(カウンタ)の値を1つ増やす」
        slide10-4.gif   Untitled-3.gif
      上の『対立遺伝子aを引いたら』(このことは対立遺伝子頻度aの確率で生じる)というのは、
       『乱数を一つ発生させたとき、それが対立遺伝子頻度aより低かったら』(これも対立遺伝子頻度aの確率で生じる)と
      同じことと考える。そこで、
     「乱数を一つ発生させたとき、それが対立遺伝子頻度aより低かったら、
        aの数を入れておく変数(counta)の値を1つ増やす」
    • ここが理解されなかったら、Rを使ってrunif(20)を表示させる。0.5より小さい値の数を数えて、次世代の集団の対立遺伝子aの頻度とし、もう一度Rを使って乱数発生。この操作を3回ぐらいやる
  4. 上の操作を、集団サイズの数だけ繰り返したら、カウンタという変数には、取り出された対立遺伝子aの数が入っているはず。この値は、次世代の集団における対立遺伝子aの数なので、次の世代に渡さなければならない。そこで、対立遺伝子aの数という変数にカウンタの値を代入して、次の繰り返しに渡す。%%%カウンタという変数は、次の繰り返しでaの数を数えるのに使われるから、内側のfor文が実施される前に0を代入して、初期化しておく
     ・内側のfor文の繰り返しが1回終わったら
      対立遺伝子aの数 = カウンタ  
      カウンタ= 0                 #内側のfor文の直前に書いても同じこと
  5. 次の世代でも上と同じことを繰り返すが、集団内の対立遺伝子頻度は、前の世代の対立遺伝子頻度(対立遺伝子aの数/集団サイズ)で決まる。

<授業後に追加>赤玉・青玉を取り出す実験でプログラミングを説明

授業では「遺伝子頻度」や「集団」という聞き慣れない言葉が出てきて、理解しにくかったのかもしれません。そこで、赤玉と青玉が同じ数だけたくさん入った筒から、20個の玉を取り出して、その中に含まれる赤玉の割合に従って筒に含まれる赤玉の個数を変え、また20個取り出す。。。こういう操作を何度も繰り返すというモデルで、プログラミングの説明を追加しておきます。

Untitled-9.gif

「赤玉と青玉の同じ数入っている筒から無作為に20個の玉を取り出したら、赤玉が7個入っていた」
  -----------------------------------------------------
上でやった操作には、「20個取り出して赤玉かどうか判断する」というところに「繰り返し」が含まれている。
「繰り返し」がはっきりと分かるように説明すると、
 1. 筒からボールを一つ取り出す
 2. そのボールが赤玉だったら数える(1つ目だったら「一つ」、2つめは「2つ」..)。青玉だったら数えない
 3.  1 と 2 の操作を20回繰り返す

この部分を大雑把でいいので、プログラミングしてみる。繰り返しはfor文で表現できる

 カウンタの値は最初は0にしておく
 for ( 20回繰り返す  ) {
   if (取り出したのが赤玉だったら) {
     カウンタの値を1増やす
   }
 }

カウンタのところは代入文で表現できる。for文も、Rで実行できる形式に直す

 カウンタ = 0
 for ( i in 1:20 ) {
   if (取り出したのが赤玉だったら) {
     カウンタ = カウンタ + 1
   }
 }

「もし、取り出したのが赤玉だったら」というif文をどうプログラミングすれば良いかが、今回のシミュレーションの最重要ポイントだろう。この部分は乱数を使って表現できる。

Untitled-10.gif
つまり、0-1の範囲の乱数を1つ発生させ、それが0.5よりも小さい値だったら、赤玉が取り出されたとみなせばいい。 乱数を1つ発生させるには runif(1) を使う。

 カウンタ = 0
 for ( i in 1:20 ) {
   if (runif(1) < 0.5 ) {
     カウンタ = カウンタ + 1
   }
 }

これで20個取り出したうちの赤玉の個数を数えるという部分は完成です

&ref(): File not found: "Untitled-11" at page "授業/H19/情報処理/11";
上の図で、最初の矢印の部分の操作を行い、筒から取り出した赤玉の数を数えるシミュレーションのプログラミングは終わりました。次に、取り出した赤玉の割合になるように筒の中の赤玉の数を変更して、同様の操作を100回繰り返すという部分のプログラミングを考えます。繰り返しが含まれているのでfor文を使って大雑把に書くと、

 for ( j in 1:100 ) {   #さっきiという変数を使ったので、ここではjを使う
   ※20個取り出して赤玉の数を数える
 }

たったこれだけです。上で作った赤玉の数を数えるプログラムを※のところに入れてみると

 for ( j in 1:100 ) {   #さっきiという変数を使ったので、ここではjを使う
   カウンタ = 0
   for ( i in 1:20 ) {
     if (runif(1) < 0.5 ) {
       カウンタ = カウンタ + 1
     }
   }
 }

ん??これでは最初に20個取り出した中に7個含まれていたので、筒の中の赤玉の割合を0.35に変更したことが反映されていませんね。取り出した赤玉の数は、内側のfor文が終わった後に、カウンタの値を20で割ることで求められますから、「赤玉の割合」という変数を一つ作って、毎回、新しい値をこの変数に入れるようにしましょう。

 赤玉の割合=0.5         #初期値は0.5
 for ( j in 1:100 ) {   #さっきiという変数を使ったので、ここではjを使う
   カウンタ = 0
   for ( i in 1:20 ) {
     if (runif(1) < 赤玉の割合 ) {
       カウンタ = カウンタ + 1
     }
   }
   赤玉の割合 = カウンタ/20  #20個の中の赤玉の割合を算出
 }

2. 対立遺伝子aの初期値を自由に変更できるように、関数にしてしまう。

ここから先、サンプルプログラムの重要な部分がところどころ空欄になっています。
実行するまえに★★のついた行では、空欄を補ってください

これでは使いにくいので、aの値を自由に変更できる関数を定義してみましょう。「関数の定義」なんていうと難しそうですが、ようするに、作ったプログラムに名前をつけて、いろいろと数値を変えて解析できるようにすることです。

drift1=function(対立遺伝子aの数){       #関数定義の始まり
繰り返し世代数 = 100        #今回は100世代観察することにする
集団サイズ = 20             #1集団内の遺伝子の総数を集団サイズと呼ぶ。今回は20。
 for(i in 1:繰り返し世代数){             
   カウンタ=___               #★★カウンタの値を初期化
   for(j in 1:集団サイズ){           
     if ( runif(1) < 対立遺伝子aの数/集団サイズ  ){  
      カウンタ=カウンタ+1         
     }
   }
   結果=____/集団サイズ   #★★結果は対立遺伝子aの頻度       
   print(結果)
   対立遺伝子aの数=カウンタ
  }
}                      #関数定義の終わり

さっきのプログラムの上下を関数を定義する関数であるfunction(){}で囲んで、オブジェクトに代入するだけですから、簡単です。この関数を実行するには、最初aの値を()の中に入力して、

drift1(10)

とします。

↑を使って何回も実行してみると、結果がいろいろ変わるのがわかりますよね。

3.結果をグラフ表示。前につかったplot()関数

では、せっかく結果が出たのだから、グラフにしてみましょう。でも、その前に、以前使ったplot()という関数を思い出してみましょう。

x=c(1,2,3,4,5)
plot(x)

plot()という関数は、ベクトルの内容をプロットしてくれるのでした。そうすると、今のdrift1関数の実行結果はどうやれば表示できるでしょうか?

関数を実行した結果が全て、一つのオブジェクトにベクトルとして格納されていれば良い。例えば:
0.6, 0.55, 0.65, 0.8, 0.95, 1, 1, 1, ...

これをやるには、上のout=count/20としているところを、次から次へとベクトルの要素を付け加えるようにすれば良さそうです。ベクトル要素を付け加える関数は前回使った

append()

です。

結果=append(結果, 対立遺伝子aの数/20)

と書けば、結果というオブジェクトのベクトルにa/20の値が追加されます。但し、

結果というベクトルはあらかじめ存在していなければ、append()は使えないのでfor文を使って繰り返しを行う前に、
results=c()

という命令で、空ベクトルを作っておく必要があります。

では、まず、結果をresults()に入れてベクトルで返すような変更を、上のプログラムに加えて見ましょう
drift2=function(対立遺伝子aの数){       #関数定義の始まり
繰り返し世代数 = 100        #今回は100世代観察することにする
集団サイズ = 20             #1集団内の遺伝子の総数を集団サイズと呼ぶ。今回は20。
結果=c()       #結果を入れる空ベクトルを作る
 for(i in 1:繰り返し世代数){             
   カウンタ=0                  
   for(j in 1:集団サイズ){           
     if ( runif(1) < 対立遺伝子aの数/集団サイズ  ){  
      カウンタ=カウンタ+1         
     }
   }
   結果=_____(____, カウンタ/集団サイズ)   #★★ここにあった 結果=カウンタ/集団サイズ  の代わりに、
                                            #append()関数でベクトル作成 
   対立遺伝子aの数=カウンタ
  }
   print(結果)                 #printの文はfor文の外に出したことに注意!理由は分かりますよね
}                      #関数定義の終わり

これで

[1] 0.50 0.60 0.80 0.70 0.55 0.50 0.35 0.40 0.20 0.30 0.1....

というベクトルが表示できました。これをplot()でグラフにするには、

plot(drift2(10))

とやるだけです。

しかし。。。ちょっと待てよ?どうせなら、関数の中にplot()を入れてしまえばすぐにグラフが表示されるのでは?

drift3=function(対立遺伝子aの数){       #関数定義の始まり
繰り返し世代数 = 100        #今回は100世代観察することにする
集団サイズ = 20             #1集団内の遺伝子の総数を集団サイズと呼ぶ。今回は20。
結果=c()       #結果を入れる空ベクトルを作る
 for(i in 1:繰り返し世代数){             
   カウンタ=0                  
   for(j in 1:集団サイズ){           
     if ( runif(1) < 対立遺伝子aの数/集団サイズ  ){  
      カウンタ=カウンタ+1         
     }
   }
   結果=append(結果, カウンタ/集団サイズ)   #ここにあった 結果=カウンタ/集団サイズ  の代わりに、
                                            #append()関数でベクトル作成 
   対立遺伝子aの数=カウンタ
  }
    ____(____)     #★★ここにあった print文の代わりに、plot()でグラフ表示
}                      #関数定義の終わり

グラフがすぐに表示できました!でも、点々では見にくいので、線で結んでみましょう。plot()を、

 plot(結果, type="l")

に置き換えるだけです。

drift4=function(対立遺伝子aの数){       #関数定義の始まり
繰り返し世代数 = 100        #今回は100世代観察することにする
集団サイズ = 20             #1集団内の遺伝子の総数を集団サイズと呼ぶ。今回は20。
結果=c()       #結果を入れる空ベクトルを作る
 for(i in 1:繰り返し世代数){             
   カウンタ=0                  
   for(j in 1:集団サイズ){           
     if ( runif(1) < 対立遺伝子aの数/集団サイズ  ){  
      カウンタ=カウンタ+1         
     }
   }
   結果=append(結果, カウンタ/集団サイズ)   #ここにあった 結果=カウンタ/集団サイズ  の代わりに、
                                            #append()関数でベクトル作成 
   対立遺伝子aの数=カウンタ
  }
    ____(____, type="l")     #★★ここにあった print文の代わりに、plot()でグラフ表示
}                      #関数定義の終わり

4.何回もの計算結果を、一つのグラフにまとめて表示したい

 上で作ったプログラムを何度も走らせてみると、グラフの形がどんどん変わります。これはこれで面白いのですが、グラフの軸が一定していないし、1つのグラフに結果が重ねて表示されないので、ちょっと感じがつかめませんね。

そこで、いっそのこと、この実験を20回ぐらい繰り返してみましょう。つまり、20個の集団について別々に観察を行って、結果を1つのグラフに表示させるわけです。同じことを何度も繰り返すのですから、for文をもう一つ、外側に作ればいいって、容易に想像がつきますよね?

drift5=function(対立遺伝子aの数){       #関数定義の始まり
観察集団数=20
繰り返し世代数 = 100        #今回は100世代観察することにする
集団サイズ = 20             #1集団内の遺伝子の総数を集団サイズと呼ぶ。今回は20。
 for(k in 1:_____){    #★★これまでのプログラムの外側に新しく作ったfor文。観察集団数繰り返す。
  結果=c()       #結果を入れる空ベクトルを作る。集団をいくつも観察するので、1集団ごとに初期化
  for(i in 1:繰り返し世代数){             
    カウンタ=0                  
    for(j in 1:集団サイズ){           
      if ( runif(1) < 対立遺伝子aの数/集団サイズ  ){  
       カウンタ=カウンタ+1         
      }
    }
    結果=append(結果, カウンタ/集団サイズ)   #ここにあった 結果=カウンタ/集団サイズ  の代わりに、
                                             #append()関数でベクトル作成 
    対立遺伝子aの数=カウンタ
   }
     plot(結果, type="l")     #ここにあった print文の代わりに、plot()でグラフ表示
 }
}                      #関数定義の終わり

さあ、走らせてみましょう。

>drift5(10)

あれ?なんか変ですね。最後は0.0で線が一本でるだけです。

なぜだか分かりますか?

こたえは、対立遺伝子aの数の値です。最初、対立遺伝子aの数は10で与えられましたが、1集団分の繰り返しが終わった後、次の繰り返しの前に、最初に与えたの数字に戻すのを忘れていました。

drift5=function(対立遺伝子aの数){       #関数定義の始まり
 対立遺伝子aの初期値=_____  #★★関数に渡された「対立遺伝子aの数」値を対立遺伝子aの初期値に保存
 観察集団数=20
 繰り返し世代数 = 100        #今回は100世代観察することにする
 集団サイズ = 20             #1集団内の遺伝子の総数を集団サイズと呼ぶ。今回は20。
 for(k in 1:観察集団数){    #これまでのプログラムの外側に新しく作ったfor文。観察集団数繰り返す。
  対立遺伝子aの数=_____     #★★繰り返しの最初に、「対立遺伝子aの数」を「対立遺伝子aの初期値」に戻す
  結果=c()       #結果を入れる空ベクトルを作る。集団をいくつも観察するので、1集団ごとに初期化
  for(i in 1:繰り返し世代数){             
    カウンタ=0                  
    for(j in 1:集団サイズ){           
      if ( runif(1) < 対立遺伝子aの数/集団サイズ  ){  
       カウンタ=カウンタ+1         
      }
    }
    結果=append(結果, カウンタ/集団サイズ)   #ここにあった 結果=カウンタ/集団サイズ  の代わりに、
                                             #append()関数でベクトル作成 
    対立遺伝子aの数=カウンタ
   }
     plot(結果, type="l")     #ここにあった print文の代わりに、plot()でグラフ表示
 }                     #新しく作ったfor文の終わり
}                      #関数定義の終わり

できました。実行してみましょう。ウーン。。。うまくいってるようなのですが、1回ずつ、グラフが変わっちゃいます。グラフを重ねて描くには、par(new=T) というコマンドをplot()の前に使います。

drift5=function(対立遺伝子aの数){       #関数定義の始まり
 対立遺伝子aの初期値=対立遺伝子aの数  #関数に渡された(()の中に書かれた)aの値を対立遺伝子aの初期値に保存
 観察集団数=20
 繰り返し世代数 = 100        #今回は100世代観察することにする
 集団サイズ = 20             #1集団内の遺伝子の総数を集団サイズと呼ぶ。今回は20。
 for(k in 1:観察集団数){    #これまでのプログラムの外側に新しく作ったfor文。観察集団数繰り返す。
  対立遺伝子aの数=対立遺伝子aの初期値     #繰り返しの最初に、対立遺伝子aの数の値を初期値に戻す
  結果=c()       #結果を入れる空ベクトルを作る。集団をいくつも観察するので、1集団ごとに初期化
  for(i in 1:繰り返し世代数){             
    カウンタ=0                  
    for(j in 1:集団サイズ){           
      if ( runif(1) < 対立遺伝子aの数/集団サイズ  ){  
       カウンタ=カウンタ+1         
      }
    }
    結果=append(結果, カウンタ/集団サイズ)   #ここにあった 結果=カウンタ/集団サイズ  の代わりに、
                                             #append()関数でベクトル作成 
    対立遺伝子aの数=カウンタ
   }
     ___(new=T)     #★★複数のグラフを重ね合わせる
     plot(結果, type="l")     #ここにあった print文の代わりに、plot()でグラフ表示
 }                     #新しく作ったfor文の終わり
}                      #関数定義の終わり

あれれ?そうか。始点の軸が固定されていないから、ずれちゃってますね。しかも、何回も繰り返すとグラフがどんどん1つの画面に重ねて表示されています。次のようにして解決しましょう。

drift5=function(対立遺伝子aの数){       #関数定義の始まり
 _______()     #★★画面にこれまでに表示されているグラフをクリアする
 対立遺伝子aの初期値=対立遺伝子aの数  #関数に渡された(()の中に書かれた)aの値を対立遺伝子aの初期値に保存
 観察集団数=20
 繰り返し世代数 = 100        #今回は100世代観察することにする
 集団サイズ = 20             #1集団内の遺伝子の総数を集団サイズと呼ぶ。今回は20。
 for(k in 1:観察集団数){    #これまでのプログラムの外側に新しく作ったfor文。観察集団数繰り返す。
  対立遺伝子aの数=対立遺伝子aの初期値     #繰り返しの最初に、対立遺伝子aの数の値を初期値に戻す
  結果=c()       #結果を入れる空ベクトルを作る。集団をいくつも観察するので、1集団ごとに初期化
  for(i in 1:繰り返し世代数){             
    カウンタ=0                  
    for(j in 1:集団サイズ){           
      if ( runif(1) < 対立遺伝子aの数/集団サイズ  ){  
       カウンタ=カウンタ+1         
      }
    }
    結果=append(結果, カウンタ/集団サイズ)   #ここにあった 結果=カウンタ/集団サイズ  の代わりに、
                                             #append()関数でベクトル作成 
    対立遺伝子aの数=カウンタ
   }
      par(new=T)     #複数のグラフを重ね合わせる
      plot(結果, type="l", ____=c(0,1))   #★★ylimを使ってy軸の目盛りを0から1に固定する
 }                     #新しく作ったfor文の終わり
}                      #関数定義の終わり

5.もっと大きな集団サイズで解析してみたい。

うえのプログラムで、集団サイズが20の場合の遺伝的浮動のシミュレーションができました。何度もグラフを描かせてみると、わりとすぐに固定してしまうことがわかります。では、もっと、集団サイズがもっと大きかったらどうなるか気になってくるでしょ?

 集団サイズを大きくするにはどうすれば良いでしょうか?プログラムの中で、集団サイズを20に指定しているところを100とかに変更してもいいのですが、どうせなら、関数の中に集団サイズも指定できればいいですよね。ついでに、繰り返しの数や世代数、対立遺伝子の数まで指定できればいうことありません。つまり、

drift(繰り返し数, 世代数, 集団サイズ, 対立遺伝子の数)

とか指定して実行できれば、好きな大きさの集団で、好きな数だけ、シミュレーションを繰り返せるわけです。さっきのプログラムで、数値で指定してあったところを変数に置き換えるだけですから、そんなに難しくなさそうです。

では、やってみましょう

drift6=function(対立遺伝子aの数,  集団サイズ, 繰り返し世代数, 観察集団数){       #関数定義の始まり
 plot.new()     #画面にこれまでに表示されているグラフをクリアする
 対立遺伝子aの初期値=_____  #★★関数に渡された「対立遺伝子aの数」を対立遺伝子aの初期値に保存
 for(k in 1:_____){    #★★これまでのプログラムの外側に新しく作ったfor文。観察集団数繰り返す。
  対立遺伝子aの数=対立遺伝子aの初期値     #繰り返しの最初に、対立遺伝子aの数の値を初期値に戻す
  結果=c()       #結果を入れる空ベクトルを作る。集団をいくつも観察するので、1集団ごとに初期化
  for(i in 1:______){        #★★繰り返し世代数だけ繰り返す(「forループを回す」とも言います)
    カウンタ=0                  
    for(j in 1:______){      #★★内側のループでは集団サイズの数だけ繰り返す  
      if ( runif(1) < 対立遺伝子aの数/集団サイズ  ){  
       カウンタ=カウンタ+1         
      }
    }
    結果=append(結果, カウンタ/集団サイズ)   #ここにあった 結果=カウンタ/集団サイズ  の代わりに、
                                             #append()関数でベクトル作成 
    対立遺伝子aの数=カウンタ
   }
      par(new=T)     #複数のグラフを重ね合わせる
      plot(結果, type="l", ylim=c(0,1))   #y軸の目盛りを0から1に固定する
 }                     #新しく作ったfor文の終わり
}                      #関数定義の終わり

できちゃいました。 あんまり大きい数字を入れると計算にすごく時間がかかります。途中で終了したいときは、 ESC キーを押してください。

最後に

昨年、同様の授業を行ったところ、受講生から次のような質問がありました。

「こういうシミュレーションのプログラミングって、これから生物学を学ぶのに、どのくらい必要なんですか?」

これまでプログラミングなどしたことの無い人にとって、とても素朴な疑問だろうと思います。そこで、次のように答えておきました。

画一的な答えは難しいです。シミュレーションそのものが直接必要になるような分野もあれば、そうでない分野もあります。
ただ、私は、プログラミングは、生物学における発展的な情報処理をする上での基本技だと考えています。
千葉大の研究室に限って言えば、プログラミング技術が無いと卒業研究が終わらないというところは
現時点では(たぶん)無いでしょう。但し、プログラミングができると、卒業研究や修士の研究でより
進んだことができたり、理解が深まったり、研究が楽になるようなところはいくつもあります。
将来、数理生物学、理論集団遺伝学、進化生態学、バイオインフォマティクスを志す人には、
敢えて、不可欠な技術だと言っておきます。

第11回授業の課題

課題1.意見調査

 下の囲みの中にあるアンケートをコピー・ペーストして、「回答:」の後に答えを書き込むこと。

*第11回授業アンケート
**氏名:
**課題への回答
-今日(6月29日)の授業の進み方は?(はやい、丁度いい、おそい)
--回答:
-今日の授業の難しさはどう感じましたか(簡単すぎ 簡単 丁度いい 難しい 難しすぎ):
--回答:
-難しいと答えた人は、特にどの点が難しかったですか?:
--回答:
-今日の授業は(よく分かった 分かった 分からなかった):
--回答:
-分からないと答えた人は、特にどの点が分からなかったですか?:
--回答:
-今日の講義で理解できなかった用語があったら挙げてください:
--回答:

課題2: 復習課題

下のプログラムは遺伝的浮動をシミュレートするRのプログラムです。このシミュレーションで観察する集団の数(num_repeats)、1集団あたり観察する世代数(num_generations)、集団サイズ(遺伝子数)(size_population)、0世代目における対立遺伝子aの数(num_a_allele)を与えて実行すれば、、ランダムに選ばれる遺伝子が、次世代の集団の遺伝子頻度にどのように影響するかを、シミュレートする集団数(観察集団数)の数だけ色つきの折れ線グラフで表示することができます。下のプログラムを見て以下の問に答えなさい。

drift= function(num_repeats,num_generations,size_population, num_a_allele){
	 results=c()
        a=num_a_allele
	 for(i in 1:【1】){
		  for(j in 1:【2】){
		      count_a=0
		       for(k in 1:size_population){
			    if ( runif(1) < a/【3】 ){
			     count_a=count_a+1
			    }
		       }
		       a=count_a
		       results=append(results, a/size_population)
		  }
		  a=num_a_allele
	 }
	 rmatrix=matrix(results, nrow=num_generations, ncol=num_repeats)
	 return(matplot(rmatrix, type="l"))
}

予習課題:HTMLとウェブ&プログラミング

以下の設問に答えなさい


時間があったら解説: いろんなシミュレーション

7.丸耳の表現型(aaで示されるもの)の集団内での表現型頻度をグラフで表示

drift6=function(x, y, z, a){ 
 plot.new()         #ついでに、図を新しく書き直すという関数をここに追加しました
 firsta=a            
 for(i in 1:x){     #上のプログラムでは20でしたが、ここでは、x回繰り返す。
   a=firsta         
   results=c()
   results_p=c()       
   for(j in 1:y){              
     counta=0
     count_p=0
     count_p_type=0        
     for(k in 1:z){ 
       if ( runif(1) < a/z  ){  
        counta=counta+1
       }
       if (k%%2 != 0){              #jが奇数のとき
          if ( runif(1) < a/z  ){     #もし2つの遺伝子のうち最初のものがaならば
            count_p=count_p+1         #count_pの値を1増やす
          }
       } else {                     #jが偶数のとき
          if ( runif(1) < a/z  ){    #もし2つの遺伝子のうち最初のものがaならば
            count_p=count_p+1        #count_pの値を1増やす
          }
          if(count_p == 2){         #連続する奇数回・偶数回で合計したpの値が2ならば(つまり、aが2個だったら)
            count_p_type=count_p_type+1  #count_p_typeの値を1増やす(aaとう表現型が1つあると数える)
          }
          count_p=0                 #2個数え終わったら、count_pの値を初期化
       }
     }
     results=append(results, a/z)    
     results_p=append(results_p, count_p_type/(z/2))  #results_pに集団内のaaの頻度(count_p_type/z/2)を入れる
     count_p_type=0                                 #count_p_typeの値を初期化
     a=counta                   
    }
      par(new=T)     
      plot(results, type="l", ylim=c(0,1))  
      par(new=T)     
      plot(results_p, type="l", ylim=c(0,1), col="blue")  #aaの頻度を青線で表示する
 }                        
}

8.選択係数sで、丸耳の表現型(aa)を集団から排除してみよう

 drift6=function(x, y, z, a, s){ #選択係数sでaaが集団から排除される
 plot.new()         
 firsta=a            
 for(i in 1:x){     
   a=firsta         
   results=c()
   results_p=c()       
   for(j in 1:y){              
     counta=0
     count_p=0
     count_p_type=0        
     for(k in 1:z){ 
       if (k%%2 != 0){              #jが奇数のとき
          if ( runif(1) < a/z  ){     #もし2つの遺伝子のうち最初のものがaならば
            count_p=count_p+1         #count_pの値を1増やす
          }
       } else {                     #jが偶数のとき
          if ( runif(1) < a/z  ){    #もし2つの遺伝子のうち2番目のものがaならば
            count_p=count_p+1        #count_pの値を1増やす
          }
          if(count_p == 2){         #連続する奇数回・偶数回で合計したpの値が2ならば(つまり、aが2個だったら)
            if (runif(1) < s) {      #遺伝子型がaaのとき乱数を一つとって、それが選択係数より小さければ
              k=k-2                  #jの値から2を引いて
              count_p=0              #count_pを初期化し           
              next                   #強制的に次の繰り返し
            } else {
              count_p_type=count_p_type+1  #count_p_typeの値を1増やす(aaという表現型が1つあると数える)
            }
         }
          counta=counta+count_p
          count_p=0                 
       }
     }
     results=append(results, a/z)    
     results_p=append(results_p, count_p_type/(z/2))  #results_pに集団内のaaの頻度(count_p_type/z/2)を入れる
     count_p_type=0                                 #count_p_typeの値を初期化
     a=counta                   
    }
      par(new=T)     
      plot(results, type="l", ylim=c(0,1))  
      par(new=T)     
      plot(results_p, type="l", ylim=c(0,1), col="blue")  #aaの頻度を青線で表示する
 }                        
}

9.固定までにかかった世代数を求める関数

wtime=function(x, y, z, a){ 
firsta=a            
results_w=c()  
for(i in 1:x){     
  a=firsta         
  results=c()     
  w_time=0
  for(j in 1:y){              
    counta=0                  
    for(k in 1:z){ 
      if ( runif(1) < a/z  ){  
       counta=counta+1         
      }
    }
    a=counta                   
    if((w_time==0) && ((a==0) || (a==z))) { #もしw_timeが空で、aが0かzに固定しているとき、
       w_time=j                             #w_timeにそのときに世代数jを入れる
    }
  }
  results_w=append(results_w, w_time)       #その集団における固定までの世代時間をresults_wに入れる
}  
  return(results_w)                      
}

10.有限集団におけるヘテロ接合体頻度の変化

drift6=function(x, y, z, a){ 
 plot.new()        
 firsta=a
 resultp_m=c()            
 for(i in 1:x){    
   a=firsta         
   results=c()
   results_p=c()       
   for(j in 1:y){              
     counta=0
     count_p=0
     count_p_type=0        
     for(k in 1:z){ 
       if ( runif(1) < a/z  ){  
        counta=counta+1
       }
     }
     results=append(results, a/z)    
     results_p=append(results_p, 2*(a/z)*(1-a/z))  #results_pに集団内のaaの頻度(count_p_type/z/2)を入れる
     a=counta                   
    }
    resultp_m=append(resultp_m, results_p)
      par(new=T)     
      plot(results, type="l", ylim=c(0,1), col="gray")  
 }                        
    rmatrix=matrix(resultp_m, nrow=y, ncol=x)
    hetero=c()
    for(i in 1:y) {
      hetero=append(hetero,mean(rmatrix[i,]))
    }
    par(new=T)
    plot(hetero, type="l", ylim=c(0,1), col="blue")
    par(new=T)
    curve(0.5*(1-1/z)^x, xlim=c(0,y), ylim=c(0,1), col="red")
}