授業/H19/情報処理/11

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

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

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

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

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

• 括弧記号の使い方。
• for( i in 1:20)は繰り返す回数を示すと書いてあるのに、iの値が一ずつ増えていくのがなぜだかわからない

  説明不足だったかもしれませんね。for(i in 1:20){  }は「iを１つずつ増やす操作を20回行いながら{}の中を実行するという意味です。

• ベクトル

  (1, 2, 3, 5)など、複数の数値の集まりのことです。配列とも言います。

• 定義するという意味が理解できない

  Rがもともと持っていない関数を新たに作ることを定義すると言っています。

• 数学が苦手なので、関数を理解するのが大変でした。
• プログラミングのあたり
• プロミラングのところ（特にappendの使い方。）

  以上3件は今日の授業で出てくるでしょう

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

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

R editorを使う †

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

• ファイルメニューから「新しいスクリプト」をクリック

  R Editor というエディタのウィンドウが開く

  このウィンドウの中に授業ページのサンプルプログラムをペーストして編集
  
  

• 編集が終わったら、実行したい部分をマウスで選択して右クリック

  選択した部分のプログラムが実行され、結果がRコンソールに表示される。

繰り返し †

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) <=

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

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

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

• 参考資料：http://bean.bio.chiba-u.jp/joho/index.php?joho19%2FPROJECT#n4f5442c

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

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

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

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

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

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

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

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

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

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

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

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

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

プログラムを実行 †

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

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

集団サイズ = __             #１集団内の遺伝子の総数を集団サイズと呼ぶ。今回は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の個数を入れるカウンタという変数の値を１つ増やす
     }
   }
 結果=c(i, 対立遺伝子aの数/集団サイズ)  #何世代目かを表す数字(i)と対立遺伝子aの頻度を結果というオブジェクトに入れる
 print(結果)                      #outの内容を画面に表示する
 対立遺伝子aの数=カウンタ         # 対立遺伝子aの数を新しい世代の対立遺伝子aの個数(カウンタの値)で置き換える
}

• 変数名に日本語を使わない場合：

  numa = 10                   #対立遺伝子aの集団中の個数を numa という変数で表す。今回は10個がaであるとする
  numgen = 100                #何世代について観察するかを numgen で表す。今回は100世代観察することにする
  for(i in 1:numgen){            #この行から対応する}までの間を numgen で指定された100世代分繰り返す
   counta=0                        #次世代にランダムに残るaの値を counta という変数に入れる。初期化して0にする。
     for(j in 1:20){               #数千億からなる配偶子プールから1個の配偶子を取り出す試行を20回行う
       if ( runif(1) < numa/20  ){    #乱数を一つ発生させ、それが対立遺伝子aの頻度(a/20)よりも小さければ
         counta=counta+1              #aの個数を入れるcountaという変数の値を１つ増やす
       }
     }
   out=c(i, counta/20)             #何世代目かを表す数字(i)と対立遺伝子aの頻度(a/20)をoutというオブジェクトに入れる
   print(out)                      #outの内容を画面に表示する
   numa=counta                     # aの値を新しい世代の対立遺伝子aの個数(counta)で置き換える
  }

質問しますので挙手して答えてください。
100世代目が0だった人、　1だった人、　それ以外の数値だった人
固定した世代数　0-10  10 - 20  20-50  51-75  76-100  100<]

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

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

繰り返し
代入
条件分岐

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

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

だけです。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1. まず、100世代この集団を観察するのだから、繰り返し文で100回繰り返す構造が必要になる。

   ＊100回繰り返す繰り返し文の構造
   for(i in 0:100) {
         ＜※ここに繰り返される内容が入る＞
   }

2. １世代で行われる操作（上の※の部分）は、１つの集団から配偶子を親の集団と同数（つまり、20個）取り出すという繰り返し操作。なので、上の100回繰り返すfor文の中には、さらに、集団中の遺伝子の数だけ繰り返すfor文が入る

   ※集団中の遺伝子数分、つまり、20回繰り返す文
      for(j in 0:20) {
      　　＜※※ここに繰り返される内容が入る＞
      }

3. 上のfor文の中では、配偶子プールか集団サイズの数だけ取り出した配偶子の中に、対立遺伝子aがいくつ含まれているかを数えること。そのとき、

   　※※の操作：
     「配偶子を一つ取り出すとき、遺伝子頻度の確率で、対立遺伝子aをひく。もし、対立遺伝子aを引いたら、
       aの数を入れておく変数(カウンタ)の値を１つ増やす」

   　 ＝　　　　　

     上の『対立遺伝子aを引いたら』（このことは対立遺伝子頻度aの確率で生じる）というのは、
      『乱数を一つ発生させたとき、それが対立遺伝子頻度aより低かったら』（これも対立遺伝子頻度aの確率で生じる）と
   　　同じことと考える。そこで、
   　「乱数を一つ発生させたとき、それが対立遺伝子頻度aより低かったら、
       aの数を入れておく変数(counta)の値を１つ増やす」

   • ここが理解されなかったら、Rを使ってrunif(20)を表示させる。0.5より小さい値の数を数えて、次世代の集団の対立遺伝子aの頻度とし、もう一度Rを使って乱数発生。この操作を3回ぐらいやる
4. 上の操作を、集団サイズの数だけ繰り返したら、カウンタという変数には、取り出された対立遺伝子aの数が入っているはず。この値は、次世代の集団における対立遺伝子aの数なので、次の世代に渡さなければならない。そこで、対立遺伝子aの数という変数にカウンタの値を代入して、次の繰り返しに渡す。%%%カウンタという変数は、次の繰り返しでaの数を数えるのに使われるから、内側のfor文が実施される前に0を代入して、初期化しておく

    ・内側のfor文の繰り返しが１回終わったら
     対立遺伝子aの数 = カウンタ  
     カウンタ= 0                 #内側のfor文の直前に書いても同じこと

5. 次の世代でも上と同じことを繰り返すが、集団内の対立遺伝子頻度は、前の世代の対立遺伝子頻度（対立遺伝子aの数/集団サイズ）で決まる。

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

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

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

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

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

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

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

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

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

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

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

上の図で、最初の矢印の部分の操作を行い、筒から取り出した赤玉の数を数えるシミュレーションのプログラミングは終わりました。次に、取り出した赤玉の割合になるように筒の中の赤玉の数を変更して、同様の操作を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個取り出した中に含まれていた赤玉の数の割合に、筒の中の赤玉の割合を変更したことが反映されていませんね。取り出した赤玉の数は、内側の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個の中の赤玉の割合を算出
 }

あとは結果をプリントする文を入れれば完成です。何回目に固定したかが分かるように、j の値も一緒にプリントするようにしましょう。

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

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

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

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

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

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

drift1(10)

とします。

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

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

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

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

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

drift5=function(対立遺伝子aの数){       #関数定義の始まり
観察集団数=20
繰り返し世代数 = 100        #今回は100世代観察することにする
集団サイズ = 20             #１集団内の遺伝子の総数を集団サイズと呼ぶ。今回は20。
 for(k in 1:_____){    #★★これまでのプログラムの外側に新しく作ったfor文。観察集団数繰り返す。
  結果=c()       #結果を入れる空ベクトルを作る。集団をいくつも観察するので、１集団ごとに初期化
  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で与えられましたが、１集団分の繰り返しが終わった後、次の繰り返しの前に、最初に与えたの数字に戻すのを忘れていました。

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

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

drift5=function(対立遺伝子aの数){       #関数定義の始まり
 対立遺伝子aの初期値=対立遺伝子aの数  #関数に渡された（()の中に書かれた）aの値を対立遺伝子aの初期値に保存
 観察集団数=20
 繰り返し世代数 = 100        #今回は100世代観察することにする
 集団サイズ = 20             #１集団内の遺伝子の総数を集団サイズと呼ぶ。今回は20。
 for(k in 1:観察集団数){    #これまでのプログラムの外側に新しく作ったfor文。観察集団数繰り返す。
  対立遺伝子aの数=対立遺伝子aの初期値     #繰り返しの最初に、対立遺伝子aの数の値を初期値に戻す
  結果=c()       #結果を入れる空ベクトルを作る。集団をいくつも観察するので、１集団ごとに初期化
  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文の終わり
}                      #関数定義の終わり

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

• y軸の目盛りを0から1に固定するには、plot()関数のylimというオプションを使います
• グラフがどんどん重なって見にくくなったら、

  >frame()

  　か

  >plot.new()

  というコマンドを使います。これもプログラムの中に組み込んでしまいましょう。

drift5=function(対立遺伝子aの数){       #関数定義の始まり
 _______()     #★★画面にこれまでに表示されているグラフをクリアする
 対立遺伝子aの初期値=対立遺伝子aの数  #関数に渡された（()の中に書かれた）aの値を対立遺伝子aの初期値に保存
 観察集団数=20
 繰り返し世代数 = 100        #今回は100世代観察することにする
 集団サイズ = 20             #１集団内の遺伝子の総数を集団サイズと呼ぶ。今回は20。
 for(k in 1:観察集団数){    #これまでのプログラムの外側に新しく作ったfor文。観察集団数繰り返す。
  対立遺伝子aの数=対立遺伝子aの初期値     #繰り返しの最初に、対立遺伝子aの数の値を初期値に戻す
  結果=c()       #結果を入れる空ベクトルを作る。集団をいくつも観察するので、１集団ごとに初期化
  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文の終わり
}                      #関数定義の終わり

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

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

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

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

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

では、やってみましょう

drift6=function(対立遺伝子aの数,  集団サイズ, 繰り返し世代数, 観察集団数){       #関数定義の始まり
 plot.new()     #画面にこれまでに表示されているグラフをクリアする
 対立遺伝子aの初期値=_____  #★★関数に渡された「対立遺伝子aの数」を対立遺伝子aの初期値に保存
 for(k in 1:_____){    #★★これまでのプログラムの外側に新しく作ったfor文。観察集団数繰り返す。
  対立遺伝子aの数=対立遺伝子aの初期値     #繰り返しの最初に、対立遺伝子aの数の値を初期値に戻す
  結果=c()       #結果を入れる空ベクトルを作る。集団をいくつも観察するので、１集団ごとに初期化
  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回授業の課題 †

• 提出期限：7月3日水曜正午（下記課題全て）
  • 提出期限を過ぎたものでも、点数を半分にするなどで評価しています。
• http://bean.bio.chiba-u.jp/joho19/ に、「自分のID」/11　という新しいページを作成し、下の囲みの中にあるアンケートをコピー・ペーストして、「回答：」の後に答えを書き込むこと。
  • 手順
  1. 個人ページのトップページ（上の方のページタイトルに、「joho19/自分のID」のみが書かれているページ）、画面の上の方にある〔 新規 〕をクリック
  2. ページ名を尋ねる入力スペースが表示されるので、半角英数字で、ドット・スラッシュ・1・1を下のように入力

     ./11

     • 注：課題提出ページが正しく作れていない場合、課題の点数から１点減点です

       良くある間違い： joho19/07s9999/06 というページを作るべきなのに joho19/07s9999/05/06 としてしまったとか

• 今回の課題はいつもより２点多い９点満点で採点します(予習6点、復習3点）
  • 復習課題の方は理解できないところがあったらメールでも、自分のページでもいいので、気軽に質問してください。

課題１．意見調査 †

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

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

課題2: 復習課題 †

下のプログラムは遺伝的浮動をシミュレートするRのプログラムです。このシミュレーションで観察する集団の数(num_repeats)、１集団あたり観察する世代数(num_generations)、集団サイズ（遺伝子数）(size_population)、０世代目における対立遺伝子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"))
}

• 問1: プログラム中【1】〜【3】に入る変数名を、下のリストから選んで答えなさい

  num_repeats, num_generations, size_population, num_a_allele

  • それぞれの変数は次のような意味を持っています
    • num_repeats : シミュレーションを行う集団数（１つの集団についてある世代数だけ観察するという実験を、何回（何個の集団について）おこなうか）、すなわち観察集団数。何本の折れ線グラフを表示させるかはこの値で決まる。
    • num_generations : 1つの集団について何世代続けて観察するか。x軸の長さになります。
    • size_population : 集団のサイズ、つまり、１つの集団に含まれる遺伝子の総数
    • num_a_allele : aアリルの数（つまり、この値をnum_populationで割ったものが、0世代目でのaアリルの遺伝子頻度になる）
• 問2: 上のプログラム中、【1】〜【3】を問1で答えた変数で置き換え、次の２つの場合を実行しなさい。

  1. 実験の繰り返し数10回、1集団での観察世代数500世代、集団サイズ（遺伝子数）200, 初期状態での対立遺伝子ａの数100
  2. 実験の繰り返し数10回、1集団での観察世代数100世代、集団サイズ（遺伝子数）20, 初期状態での対立遺伝子ａの数10
  いずれの場合についても、できたグラフをWinShotで切り取って、課題提出ページに貼り付けなさい（注：1番目のシミュ
  レーションの実行には、数分かかるかもしれません）（注：あまり大きい図（画面全面に広がる図とか）にはしないでく
  ださい）

• 問3: 上の２つのグラフを比較して分かることを述べなさい。
• 評価：
  • この課題は6点で評価します
  • Rによるプログラミングとシミュレーションを行うことができるかどうか

予習課題：HTMLとウェブ＆プログラミング †

以下の設問に答えなさい

• 問1：あなたが今みているウェブページの本体はテキストファイルです。でも、ウェブページの画面には写真などが表示されています。文字情報だけのはずのテキストファイルで写真が表示できるのはどうしてでしょうか？正しい答えを１つ以下の選択肢から選びなさい。
  • a. ウェブページは写真をファイル内に保存できる、特殊な形式のテキストファイルである。
  • b. ウェブブラウザがテキストファイルの指定に従って、別の場所に保存された写真を画面上の指定された位置に表示している
  • c. ウェブサーバというものがテキストと写真の入った画像ファイルを合成して発信している
• 問2：皆さんの課題の採点結果閲覧ページは、PostgreSQLというデータベースシステムと、PHPというプログラミング言語が使われています。自分の成績ページを表示させるときに、どのような事が起きているでしょうか？
  • a. インターネットにアクセスしているコンピュータに暗号化された全員分のデータがダウンロードされ、ブラウザの持っている機能で自分のデータだけ解読され、表示されている。
  • b. アクセス先のウェブサーバーがデータベースサーバーからあなたのデータを取り出し、１回１回、あなた専用のウェブページを作成している。
  • c. 教員が毎回、課題の採点後に、一人一人のための返却用ページを作成し、ウェブサーバーにアップロードしている。アカウントとパスワードを入力すると、成績データのページが納められたディレクトリに入った自分専用のページが表示される。
• 問3：あなたが普段見ているウェブページは、HTMLという書式で作られており、タグと呼ばれるマークで、下線や太字を指定したり、他のページへのリンクを指定したりできます。HTML書類の中で、次のようなタグで指定された部分がありました。この部分は、ブラウザではどのように表示されるでしょうか？

  <font size="3"><i><b>HTMLの世界にようこそ！</b></i></font>

• 評価：
  • この課題は3点で評価します
  • HTMLについて自分で調べて予習したかどうか
  • どの処理がサーバー側で行われ、どの処理がローカル（自分が操作しているコンピュータ）で行われているかを理解しているかどうか

---

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

７．丸耳の表現型(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  ){     #もし２つの遺伝子のうち最初のものがaならば
            count_p=count_p+1         #count_pの値を１増やす
          }
       } else {                     #jが偶数のとき
          if ( runif(1) < a/z  ){    #もし２つの遺伝子のうち最初のものがaならば
            count_p=count_p+1        #count_pの値を１増やす
          }
          if(count_p == 2){         #連続する奇数回・偶数回で合計したpの値が2ならば（つまり、aが2個だったら）
            count_p_type=count_p_type+1  #count_p_typeの値を1増やす（aaとう表現型が1つあると数える）
          }
          count_p=0                 #２個数え終わったら、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の頻度を青線で表示する
 }                        
}

８．選択係数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  ){     #もし２つの遺伝子のうち最初のものがaならば
            count_p=count_p+1         #count_pの値を１増やす
          }
       } else {                     #jが偶数のとき
          if ( runif(1) < a/z  ){    #もし２つの遺伝子のうち２番目のものがaならば
            count_p=count_p+1        #count_pの値を１増やす
          }
          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の頻度を青線で表示する
 }                        
}

９．固定までにかかった世代数を求める関数 †

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)                      
}

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

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")
}