*Rを使ったプログラミング演習 [#maa24e16]
#contents
前回の授業で使ったRという計算ソフト。アンケートの結果を見ると「感動」してくれた人もいるぐらい、面白かったよう。コマンドを入力するのはちょと難しいけど、エクセルみたいにマウスを使って面倒な操作をしなくても、一瞬にしてグラフが描けてしまう。これなら、良いレポートを仕上げるために、何度も解析をやりなおして、グラフを描き直すっていう作業も面倒ではなさそう。。。
初めて使ったソフトなのに、皆さん、けっこうRの操作には慣れているし（ほとんどみんな、課題ができていた）、Rに魅力も感じているよう。なら、もうちょっと、Rと親しくなろうということで、今回は、Rを使ってプログラミングを勉強します。
*Rのおさらい [#j465be3c]
前回やったRを使った計算のおさらいを、サンプルデータをつかってやってみましょう。
**サンプルデータ： 仲岡先生の実習から [#ca1098d1]
下のデータは仲岡先生の実習でやった、調査コドラートの面積と、出現種数です
コドラート1 10*10 2種
コドラート2 15*15 4種
コドラート3 20*20 7種
コドラート4 25*25 7種
コドラート5 30*30 7種
コドラート6 35*35 9種
コドラート7 40*40 8種
コドラート8 45*45 12種
コドラート9 50*50 10種
コドラート10 55*55 10種
コドラート11 60*60 11種
コドラート12 65*65 11種
コドラート13 70*70 9種
コドラート14 80*80 10種
コドラート15 100*100 12種
（データ提供、してくれた方に感謝！）
Rを使ってこのデータの、面積と種数の相関関係を知りたいとき、必要なのは、コドラートの１辺の長さと種数のベクトルです。
&size(16){今から何をするか、もうお分かりですよね？};
予想通り、K2Editorを立ち上げて、要らないところを消し、
10 2
15 4
. .
というデータに整形して、エクセルにコピペして、１列ごとのデータを scan() 関数を使ってオブジェクトに読み込みます。では、やってみましょう。（そろそろミミタコだろうから、説明は簡単に書いておきます）
***1. K2Editorでの操作 [#p44e93f7]
++上の囲みのデータをコピペ
++正規表現置換
検索文字列： ^.*\* （説明：行頭から最初のアスタリスクまで） 置換文字列：＜無し＞（つまり削除ということ）
検索文字列： 種 置換文字列：＜無し＞
検索文字列： + （説明：左に書かれているのは半角のスペースと半角のプラス） 置換文字列： \t （説明：タブ）
++終わったら全体を選択してコピー。 エクセルを起動
***2. エクセルでの操作 [#t0e23894]
++K2Editorで整形したデータをペースト
++２列に分かれてデータが表示される
10 2
15 4
. .
++1列のデータをコピーして、Rに移動
***3. Rでの操作 [#h2054842]
++scan()関数を使って、コドラートの一辺の長さを qに、種数をyに代入
> q=scan()
[1] ← ここでペースト。yについても同様
> y=scan()
++（確認のため、それぞれ表示しておきます）
> q
[1] 10 15 20 25 30 35 40 45 50 55 60 65 70 80 100
> y
[1] 2 4 7 7 7 9 8 12 10 10 11 11 9 10 12
++欲しいのは面積なので、xにqの2乗を入れる
x=q^2
++それぞれの常用対数を使ってグラフにプロット（対数グラフ用紙必要無し！⌣)
> plot(log10(x), log10(y))
++回帰分析
> result=lm(log10(y) ~ log10(x))
> result
Coefficients:
(Intercept) log10(x)
-0.1375 0.3250
ここに表示された結果は、log(y)=0.325*log10(x)-0.1375 という回帰直線の式
++回帰直線の描画
> abline(result)
++解析のサマリー表示
> summary(result)
Call:
lm(formula = log10(y) ~ log10(x))
Residuals:
Min 1Q Median 3Q Max
-0.211429 -0.054045 0.006334 0.053616 0.142160
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.13750 0.15254 -0.901 0.384
log10(x) 0.32498 0.04713 6.896 1.09e-05 *** ←有意水準0.1%で有意
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ←これ、有意水準
Residual standard error: 0.09927 on 13 degrees of freedom
Multiple R-Squared: 0.7853, Adjusted R-squared: 0.7688
F-statistic: 47.55 on 1 and 13 DF, p-value: 1.092e-05
++相関の検定
Pearson's product-moment correlation
data: log10(x) and log10(y)
t = 6.8957, df = 13, p-value = 1.092e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.6847590 0.9618161
sample estimates:
cor
0.886173 ← 相関係数
この授業では、それぞれの分析の意味については説明しません。仲岡先生に聞いて下さい。（グラフ用紙で提出しなくてもいいと仰っていたような、いなかったような。。。）
＊大事なこと＊ 上のような解析の方法、私だってRでできるとは知りませんでした。では「回帰分析行いなさい」なんて言われたら、どうするか？もしかしてRでできるかなと思ったら、
マニュアルを読みましょう
--http://cse.naro.affrc.go.jp/takezawa/r-tips.pdf には竹澤さんという方が書かれた詳しい説明があります。私がやったことは、Adobe Acrobat Readerで開いて、検索機能を使って「回帰分析」が載っているページを調べただけです。あとは、例を見ながら、上の解析をやってみました。
--http://cse.naro.affrc.go.jp/takezawa/r-tips.pdf には船尾さんという方が書かれた詳しい説明があります。私がやったことは、Adobe Acrobat Readerで開いて、検索機能を使って「回帰分析」が載っているページを調べただけです。あとは、例を見ながら、上の解析をやってみました。
*プログラミングって何？ [#a076c86e]
ところで、プログラミングって何をすることでしょうか？IT用語辞典で調べてみると、、、
載っていません
唯一近いところでは、「プログラミング言語」というのがあります。説明文からプログラミングに関わるところを抜き出して、「プログラミング」の説明文を作ってみると、
-''プログラミング'': 人間に理解できるように英語などを元に作られたプログラミング言語を用いて、そのままではコンピュータが実行させる命令を作ること。できた命令はソースコード、プログラム、スクリプトなどの名前で呼ばれる。命令文をコンピュータに実行させるためには、コンピュータが理解可能な機械語の羅列(オブジェクトコード)に翻訳する必要がある。プログラミング言語によっては、コンパイルという翻訳作業を別にしなければならないコンパイラ型言語と、リアルタイムに翻訳作業をしてコンピュータに命令を渡すインタプリタ型言語がある。
今日の授業でコンピュータにあれこれ命令を与えるために使うRは、計算ソフトですが、インタプリタ型のプログラミング言語ということもできます。前回インストールしたPHPも同じくインタプリタ型のプログラミング言語です。
つまり、プログラミングとは、簡単に言うと、「%%%言葉を使ってコンピュータが理解できる命令を作ること%%%」 です。
**コンピュータにどんな命令ができるのか？ or どんな命令がしたいか？ [#hb1af07f]
「Wordを立ち上げて文章を作成し、印刷する」とか、「Excelで家計簿を管理する」とかは、市販のソフトウェアを使った作業ですが、いずれも、非常に高度で複雑な命令が、使用者からコンピュータに向かって発せられています。マウスを動かしたらポインタが動くなんていうのは、けっこう複雑な命令ですよ。でも、ほとんどの場合、一般の人がコンピュータでやりたいと思うほとんどの作業は、わざわざ「命令」なんていうものを意識する必要は無いはずです。利用者が理解しやすく・使いやすいように、専用のソフトウェアが開発され、販売されているでしょう。
生物学の研究者でも、プログラミングなどを経験したことの無い人は大勢います。だって、今や、DNAシーケンスの決定や、整列や、系統樹作成にだって、専用のソフトウェアがあるのですから。。。それぞれのソフトウェアの使い方を覚えれば、それですむわけです。私たちも新しい解析方法が開発されるたびに、そのソフトウェアの使い方を覚えるのに苦労しています。
でも、それだけで、本当にいいのでしょうか？良くないという点が、まず、２つ思いつきます。
1. 自分の目的にあったソフトウェアがいつも存在するというわけではない
2. ソフトが存在する場合、目的の数だけ、ソフトの使い方を覚えなければならない
ということです。特に１番目は致命的ですね。ソフトウェアが開発されてなかったら、目的が達成できないのですから。例えば、集団遺伝学で学ぶ、遺伝的浮動の解析ソフトなんていうのは、パソコンショップに行っても売ってないです（でも、じつは、このくらいメジャーなものなら、ネットで探せばフリーのものがあります）。
そこで、この授業でプログラミングを勉強する理由は次の３点です。
+市販のソフトではカバーできないような目的を、コンピュータに命令を与えることで達成する
--あなたのコンピュータは実は、すごい性能をもっています。その性能はできる限りつかわないと損
+プログラミングを行うことで、目的とする生物学的現象の理解が深まる
--コンピュータは石頭なので、物事を本当に論理立てて説明してやらないと、命令を聞いてくれません
&br;
プログラミングの過程で、生物学的現象の意味を論理立てて考える必要に迫られます
+プログラミングを行って、命令を実行させることで、初めて、コンピュータに''作業をさせている''ことが実感できる
--ワープロや表計算とかは、どちらかというと、道具を使って作業をさせられている感が強いですよね
*Rを使ったシミュレーション：集団遺伝学的解析を試してみよう！ [#if40831d]
#ref(./drift.gif,around,right,80%)
そこで、Rを使ってプログラミングを行うのですが、まず、目的が必要です。たぶん、皆さんが知っている生物学的現象で、Rをつかってやってみる一番ビジュアル的に面白いのは、集団遺伝学におけるいろんなシミュレーションでしょう。 （ただ、Rを使うデメリットは、計算速度が遅いことです。大規模なシミュレーションをやるには、インタプリタ型のソフトではつらいでしょう。）
シミュレーション（simulation, 「シュミレーションでは無いですよ！」）は、数学モデルなどを用いて現実に似た状況をコンピュータ上で再現することです。模擬実験とも呼ばれています。
きっと、右の図は見たことがあるでしょう？これは、遺伝的浮動の効果を説明するときによく使われるグラフです。このシミュレーションでは、集団サイズN=100, ２つの対立遺伝子のうち、片方の対立遺伝子の頻度を示す。最初の世代の対立遺伝子の数は50:50。このグラフを使って、遺伝的浮動とはどういうものか、説明できる人はいますか？また、このグラフはどうやって作ったものですか？
あてます
今日の授業では、Rを使って、「繰り返し」と「条件分岐」による最も単純なプログラムを、どんどん成長させていって、最終的には、上のグラフのようなものが作れるところまでを学びます。この過程で、遺伝的浮動とか、突然変異率とか、固定までの待ち時間とか、綿野さんの授業で勉強したようなことが、でてきて、たぶん、理解が深まると思います。
**Rの説明書 [#b6e9b891]
授業では一つ一つ解説しますが、復習課題をやるときとか、自分で新たな解析に挑戦するときは、説明書が必要になるでしょう。
-http://cse.naro.affrc.go.jp/takezawa/r-tips.pdf には竹澤さんという方が書かれた、日本語のRの解説があります。知りたい内容は、pdfを開いて、Acrobat Readerの検索機能でサーチし下さい。
-http://cse.naro.affrc.go.jp/takezawa/r-tips.pdf には船尾さんという方が書かれた、日本語のRの解説があります。知りたい内容は、pdfを開いて、Acrobat Readerの検索機能でサーチし下さい。
また、Rを使っている最中は
> ?plot
> example(plot)
という風に、?の後に関数名を打てば、説明や使用例が英語で表示されます。説明文の一番最後にはExample（使用例）が載っていますので、けっこう参考になります。
**プログラミングの基礎: 繰り返しと条件分岐 [#bd49f99e]
前回授業の課題で、プログラミングの基本は繰り返しと条件分岐だという話しをしました。
***繰り返し [#l4531379]
繰り返しは次のように for() を使って表現します。
for(i in 1:10) { print(i) }
この命令は、
-iを1から10まで変化させなさい。そして、それぞれのiの値の時に、{ } で囲まれた部分の命令を実行しなさい
-{ } の中の意味は、iの値を画面にプリントしなさい
というものです。 { } のところは複数行に分かれていてもよく
for(i in 1:10) {
print(i)
}
と書いても結果は同じです。
-縦に長くプリントされるのは格好悪いので、オブジェクトにベクトルとして代入してから表示させてみましょう。
a=c() #最初にNULLベクトルという空のベクトルを一つ作っておきます。
for(i in 1:10) {
a=append(a,i) #aにiの値をベクトル要素として追加します。
}
print(a) #最後に結果の内容を表示します
-for文を入れ子にすると、指数関数的に繰り返しの数が増えます
#＊＊次のプログラムの実行結果を予測してください。＊＊
for(j in 1:10) {
for(i in 1:10){
print(i*j)
}
}
-予想通りでしたか？ 縦に長く表示されるのも格好が悪いので、今日の課題で使う関数を一つかって、行列形式で出力してみましょう
#＊＊次のプログラムの実行結果を予測してください。＊＊
result=c()
for(j in 1:9) {
for(i in 1:9){
result=append(result, i*j)
}
}
matrix(result, nrow=i, ncol=j) #nrowとncolはそれぞれ行の数、列の数を示します。
-小学校のとき散々お世話になって表がでてきましたか？
-ここでresultの中身を確認して、行列表示と比較しておきましょう
> result
-２つの結果を比較してみると分かるように、matrix()という関数は、resultに入っているベクトルを、nrowとncolの数に従って、行列に変換したものです。
***条件分岐 [#r32ace1b]
条件分岐はif()を使って表現します
i=3
if(i==3){print("san")}
この例では、
-iに3を代入し
-もしiが3なら "san" とプリントしなさい
という命令です。{}の中は複数行にわかれていてもよく、また、else を使えば、それ以外の条件を指定できます
i=4
if(i==3){
print("san")
} else {
print("san deha nai")
}
***繰り返しと条件分岐 [#a28dddc9]
では、上の繰り返しと条件分岐を組み合わせて見ましょう
1から10までのうち、5以外を画面に表示しなさい
というプログラムなら
for(i in 1:10){
if(i != 5) { # != は「等しくない」ということを示します
print(i)
}
}
**初めてのプログラミング：繰り返しと条件分岐を使って [#ob71bc71]
では、次のプログラムを作ってみてください
1から10までの数字を表示させなさい、5だけは"five"と英単語で表示させなさい
**0. シミュレーションの前段階：モデルの想定 [#g314d5f0]
これからシミュレーションを行うのですが、つぎのような状況を想定してみましょう。
#ref(./mona2.gif,around,right)
#ref(./mona1.gif,around,right)
まず、対象とする架空の生物（仮に、「モナー」と呼びます&size(11){（注：モナーについて詳しく知りたい人は、ウィキペディアで「モナー」を検索してみてください。）};）は2倍体の生物で、雌雄の区別が無く、繁殖期になると集団内の個体が一斉に配偶子を数千億個ずつ放出して死亡します。配偶子は生息域の空間内をランダムに漂った後、24時間後に、最も近くにある配偶子と接合し、新しい繁殖個体になります。集団の生息域は空間的に限られているため、親世代と同数の子供個体しか、成熟個体になれません。どの子供が成熟個体になるかは、ランダムに決まるものとします。つまり、このモデルだと、世代間で集団サイズに変動は無いということです。モナーにはトンガリ耳と丸耳の２タイプが1対立遺伝子(A, a)で決まることが知られています。トンガリ耳の遺伝子型は(AA, Aa)、丸耳の遺伝子型は(aa)です。
&size(11){なお、上のモナーの説明は、私が今回のシミュレーションのために、適当に設定したものです。だって、''「任意交配をする有限サイズの生物集団から遺伝子をランダム抽出して、同サイズの次世代集団を作った」''なんて言い方をするよりは、生物を扱ってる感じがするでしょ？};
**1. 初めてのシミュレーション：10個体からなる集団における、遺伝子頻度変化のシミュレーション [#zeb9f52b]
モナー10個体からなる、１つの集団を考えてみましょう。いまこの集団で、対立遺伝子Aの頻度が0.5, 対立遺伝子aの頻度が0.5のとき（10個体だから遺伝子の総数は、20個です。そのうち半分の10個がAで、10個がaです）。集団中の配偶子の接合はランダムに起こると仮定して、100世代の間に、集団中の対立遺伝子の頻度はどのように変化するでしょうか？
&ref(./slide5.gif,30%); &ref(./slide6.gif,40%);
今から実験をするのですから、前もって予想しておいた方が面白そうです
この集団の遺伝子頻度が1か0に固定するまでに、何世代ぐらいかかると思いますか？
a=c(10) #対立遺伝子aの集団中の個数
for(i in 1:30){ #この行から対応する}までの間を 100世代分繰り返す
counta=0 #次世代にランダムに残るaの値を入れる変数を初期化
for(j in 1:20){ #数千億からなる配偶子プールから1個の配偶子を取り出す試行を20回行う
if ( runif(1) < ??? ){ # ＊＊＊ 取り出した遺伝子がaだったならば ＊＊＊
counta=counta+1 # ＊＊＊ aの個数を入れる変数の値を1増やす ＊＊＊
}
}
out=counta/20) #何世代目かという数字(i)と対立遺伝子aの頻度(a/20)をoutというオブジェクトに入れる
print(out) #outの内容を画面に表示する
a=counta # aの値を新しい世代の対立遺伝子aの個数(counta)で置き換える
}
上から順番に意味を考えてゆきましょう。
&ref(./slide1.gif,40%); &ref(./slide2.gif,40%);
&ref(./slide3.gif,40%);&ref(./slide4.gif,40%);
それでは、次のプログラムを走らせて、20個の遺伝子中対立遺伝子aが10個ある任意交配集団の、100世代の間のaの遺伝子頻度の変化をみてみましょう。
a=10 #対立遺伝子aの集団中の個数
for(i in 1:100){ #この行から対応する}までの間を 100世代分繰り返す
counta=0 #次世代にランダムに残るaの値を入れる変数を初期化
for(j in 1:20){ #数千億からなる配偶子プールから1個の配偶子を取り出す試行を20回行う
if ( runif(1) < a/20 ){ #乱数を一つ発生させ、それが対立遺伝子aの頻度(a/20)よりも小さければ
counta=counta+1 #aの個数を入れるcountaという変数の値を１つ増やす
}
}
out=c(i, counta/20) #何世代目かという数字(i)と対立遺伝子aの頻度(a/20)をoutというオブジェクトに入れる
print(out) #outの内容を画面に表示する
a=counta # aの値を新しい世代の対立遺伝子aの個数(counta)で置き換える
}
さて、皆さんの計算結果はどうでしたか？何世代目で固定された？
**講義についての連絡： [#n4cdcdd4]
※私が行う講義はあと４回です。7月20には、をPower Pointを使って班ごとにプロジェクトの成果を発表してもらいます。また、7月20日の残りの時間は、皆さんからの要望に応えて、もう一度解説して欲しい事項などを説明します。7月13日のアンケートで、希望をとりますので、優先順位を考えておいてください。
*第10回授業の課題 [#ic6d1e37]
全ての課題は、&size(16){http://bean.bio.chiba-u.jp/joho18/ に、「自分のID/10」という新しいページを作成し、これまでの提出例にならって、分かりやすく書き込むこと。あまりに分かりにくい回答は減点します。};
-提出期限：6月28日水曜正午（下記課題全て）
--提出期限を過ぎたものでも、点数を半分にするなどで評価しています。
**課題１．意見調査 [#u259e7fa]
下の囲みの中にあるアンケートをコピー・ペーストして、「回答：」の後に答えを書き込むこと。
*第10回授業アンケート
**氏名：
**課題への回答
-今日（6月22日）の授業の進み方は？（はやい、丁度いい、おそい）
--回答：
-今日の授業の難しさはどう感じましたか（簡単すぎ 簡単 丁度いい 難しい 難しすぎ）：
--回答：
-難しいと答えた人は、特にどの点が難しかったですか？：
--回答：
-今日の授業は（よく分かった 分かった 分からなかった）：
--回答：
-分からないと答えた人は、特にどの点が分からなかったですか？：
--回答：
-今日の講義で理解できなかった用語があったら挙げてください：
--回答：
-Ｒで何ができるのか、だんだん分かってきたと思います。今後自分の勉強にＲを使うつもりですか？
--回答：
-Ｒでも、Ｒ以外の言語でも、何かやってみたいプログラミングはありますか？
--回答：
-今回やったプログラミングはおもしろかったですか？
--回答：
**復習課題：Rによるプログラミング：繰り返しと条件分岐 [#yaa55602]
次の問題文の指定にあるプログラムを作成しなさい。作成したプログラムは課題提出ページに貼り付けなさい。
行頭には必ず半角の空白を入れること。
+ 1から20までの総和を計算するプログラム（?のこと）（注：繰り返し文(for())を使ってプログラムを書いてください。sum()関数は使わないこと）
+ 乱数を20個発生させて、%%%0.5未満のものだけ%%%を画面に表示させるプログラム
-評価：
--細かい指定はしないので、プログラムをどれだけ修飾してもらってもいいです
--この問題は7点で評価します
--繰り返し、条件分岐など、プログラミングについて理解できているかどうか
--Rを用いてプログラムを正しく実行できているかどうか
--プログラムの実行内容を理解できているかどうか
**予習課題：HTMLとウェブ＆プログラミング [#t8686d3c]
以下の設問に答えなさい
-問1：あなたが普段みているウェブページの本体はテキストファイルです。でも、ウェブページの画面には写真などが表示されています。文字情報だけのはずのテキストファイルで写真が表示できるのはどうしてでしょうか？
--a. ウェブページは写真をファイル内に保存できる、特殊な形式のテキストファイルである。
--b. ウェブブラウザがテキストファイルの指定に従って、別の場所に保存された写真を画面上に表示している
--c. ウェブサーバというものがテキストファイルと写真の入った画像ファイルを合成して発信している
-問2：皆さんの課題の採点結果閲覧ページは、PostgreSQLというデータベースシステムと、PHPというプログラミング言語が使われています。自分の成績ページを表示させるときに、どのような事が起きているでしょうか？
--a. インターネットにアクセスしているコンピュータに暗号化された全員分のデータがダウンロードされ、ブラウザの持っている機能で自分のデータだけ解読され、表示されている。
--b. アクセス先のウェブサーバーがデータベースサーバーからあなたのデータを取り出し、１回１回、あなた専用のウェブページを作成している。
--c. 教員が毎回、課題の採点後に、全員分の返却用ページを作成し、ウェブサーバーにアップロードしている。アカウントとパスワードを入力すると、成績データのページが納められたディレクトリに入った自分専用のページが表示される。
-問3：あなたが普段見ているウェブページは、HTMLという書式で作られており、タグと呼ばれるマークで、下線や太字を指定したり、他のページへのリンクを指定したりできます。HTML書類の中で、次のようなタグで指定された部分がありました。この部分は、ブラウザではどのように表示されるでしょうか？
<font size="16"><i><b>HTMLの世界にようこそ！</b></i></font>
-評価：
--この問題は短いので3点で評価します
--HTMLについて自分で調べて予習したかどうか
--どの処理がサーバー側で行われ、どの処理がローカル（自分が操作しているコンピュータ）で行われているかを理解しているかどうか
------
*Rを使ったプログラミング演習（２） [#p1b85cff]
＊＊素朴な疑問＊＊
「こういうシミュレーションのプログラミングって、これから生物学を学ぶのに、どのくらい必要なんですか？」
＝＝おへんじ＝
画一的な答えは難しいです。シミュレーションそのものが直接必要になるような分野もあれば、そうでない分野もあります。
ただ、私は、プログラミングは、生物学における発展的な情報処理をする上での基本技だと考えています。
千葉大の研究室に限って言えば、プログラミング技術が無いと卒業研究が終わらないというところは
現時点では（たぶん）無いでしょう。但し、プログラミングができると、卒業研究や修士の研究でより
進んだことができたり、理解が深まったり、研究が楽になるようなところはいくつもあります。
数理生物学、理論集団遺伝学、進化生態学、バイオインフォマティクスを志す人には、
敢えて、不可欠な技術だと言っておきます。
-&color(,"#eef"){&size(12){でもね、実は授業でこれを取り上げる理由は、「こんなことができる」と「体験」してもらうことです。こういう考え方や技術を今のうちに知っておけば、もしかしたら、皆さんの中から、素晴らしい解析プログラムを作っちゃう人が現れるかも知れませんからね。};};
*前回のおさらい [#id79d5cb]
それでは、次のプログラムを走らせて、20個の遺伝子中対立遺伝子aが10個ある任意交配集団の、100世代の間のaの遺伝子頻度の変化をみてみましょう。
numa=10 #対立遺伝子aの集団中の個数
for(i in 1:100){ #この行から対応する}までの間を 100世代分繰り返す
counta=0 #次世代にランダムに残るaの値を入れる変数を初期化
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)で置き換える
}
前回、全角の空白が混じっていたため、うまく動かなかったプログラムがこれです。まず、これを始めに動かしてみて、
もう一度、プログラムの内容を理解してゆきましょう。
**遺伝的浮動のプログラムをステップ・バイ・ステップで理解する [#s1b6e76c]
-&ref(slide5.gif,around,right,40%);
前回のレポートで、コピー・ペーストすれば動くのはわかるけど、どうやったらこういうのが作れるのか分からないという意見がありました。私としては、クラスのほとんどの人に、プログラミングとはどういうことをするものかを分かって欲しいので、作成過程を１つずつ、順を追って説明します。
+まず、前回説明した「モナーの生活環」をみて、シミュレーションの条件を自分で考える。最初のシミュレーションなので、簡単に、
20個の遺伝子からなる集団
0世代目の対立遺伝子aの数が10（つまり対立遺伝子aの頻度は0.5）
この集団を100世代観察する
という条件にしておく
+今、上に書いた数字に自分で考えた変数名をつける。名前を付けておくと、あとでコンピュータに命令すときに便利だから。数値が決まっているものは代入しておく。名前は、Rの内部で使われているもの以外なら、何でも良い。自分に分かりやすいようにつける。
＊自分で考えた変数名と定数
numa=10 #0世代目の対立遺伝子aの数。2世代目以降は、シミュレーション実験で抽出されたaの数が代入される
counta=0 #配偶子プールから、次世代の集団に取り出された対立遺伝子aの数。１世代が終わったあとで、
数字が入るので、とりあえず今は0を入れておく
20 #集団の大きさ（つまり、集団に含まれる遺伝子の数）
100 #繰り返しの世代数
&ref(授業/H18/情報処理/11/slide10-1.gif,40%);
&ref(授業/H18/情報処理/11/slide10-2.gif,40%);
+変数名が決まったら、プログラムの構造を考える。紙に絵をかいて、どういうプログラムなら、「モナー」の生活環を模倣できるか考えたら（ここはホワイトボードにお絵かき）、言葉にしてみる。
＊自分で考えたプログラムの構造。コマンドを使って書かなくても、
とりあえず、構造が理解できるように書いてあればよい。
&ref(授業/H18/情報処理/11/slide10-3.gif,40%);
++何世代分か、以下の操作を繰り返す
++ １つの世代の中では、配偶子プールから配偶子を１つ選ぶという操作を20回（遺伝子の数）繰り返す
++ 取り出した配偶子が対立遺伝子aだったら、変数countaの値を一つ増やす。そうでなかったら何もしない
++ 集団内で、20回配偶子を取り出す繰り返しが終わったら、対立遺伝子aの数を numaに代入する
++次世代の集団における対立遺伝子aの頻度(num/20)をプリントする
+プログラムの構造ができたのだから、次ぎに、一つ一つのステップを、コンピュータが理解できる文で表現してみる。
++まず、100世代この集団を観察するのだから、繰り返し文で100回繰り返す構造が必要になる。
＊100回繰り返す繰り返し文の構造
for(i in 0:100) {
＜※ここに繰り返される内容が入る＞
}
++１世代で行われる操作（上の※の部分）は、１つの集団から配偶子を親の集団と同数（つまり、20個）取り出すという繰り返し操作。なので、上の100回繰り返すfor文の中には、さらに、集団中の遺伝子の数だけ繰り返すfor文が入る
※集団中の遺伝子数分、つまり、20回繰り返す文
for(j in 0:20) {
＜※※ここに繰り返される内容が入る＞
}
++上のfor文の中では、配偶子プールからpopsize個取り出した配偶子の中に、対立遺伝子aがいくつ含まれているかを数えること。そのとき、
※※の操作：
「配偶子を一つ取り出すとき、遺伝子頻度の確率で、対立遺伝子aをひく。もし、対立遺伝子aを引いたら、
aの数を入れておく変数(counta)の値を１つ増やす」
&size(24){＝}; &ref(授業/H18/情報処理/11/slide10-4.gif,40%); &ref(授業/H18/情報処理/11/slide10-5.gif,40%);
上の『対立遺伝子aを引いたら』（このことは対立遺伝子頻度aの確率で生じる）というのは、
『乱数を一つ発生させたとき、それが対立遺伝子頻度aより低かったら』（これも対立遺伝子頻度aの確率で生じる）と
同じことと考える。そこで、
「乱数を一つ発生させたとき、それが対立遺伝子頻度aより低かったら、
aの数を入れておく変数(counta)の値を１つ増やす」
---ここが理解されなかったら、Rを使ってrunif(20)を表示させる。0.5より小さい値の数を数えて、次世代の集団の対立遺伝子aの頻度とし、もう一度Rを使って乱数発生。この操作を3回ぐらいやる
++上の操作を、20回繰り返したら、countaという変数には、取り出された対立遺伝子aの数（countaという名前にしてある）が入っているはず。この値は、次世代の集団における対立遺伝子aの数なので、次世代（つまり、外側のfor文の次の繰り返し）に渡さなければならない。そこで、numaという変数にcountaの値を代入して、次の繰り返しに渡す。countaという変数は、次の繰り返しでまた、0からaを数えなければならないから、ここで0を代入して、初期化しておく
・内側のfor文の繰り返しが１回終わったら
numa = counta
counta = 0
++次の世代でも上と同じことを繰り返すが、集団内の対立遺伝子頻度は、前の世代の対立遺伝子頻度で決まる。
**2. 対立遺伝子aの初期値を自由に変更できるように、関数にしてしまう。 [#kd03c7e0]
これでは使いにくいので、aの値を自由に変更できる関数を定義してみましょう。「関数の定義」なんていうと難しそうですが、ようするに、作ったプログラムに名前をつけて、いろいろと数値を変えて解析できるようにすることです。
drift1=function(a){ #関数定義の始まり
#上にあった、a=10 という命令は削除。この関数では、aの値を自由に変えられる
for(i in 1:30){
counta=0
for(j in 1:20){
if ( runif(1) < a/20 ){
counta=counta+1
}
}
out=counta/20
print(out)
a=counta
}
} #関数定義の終わり
さっきのプログラムの上下を関数を定義する関数であるfunction(){}で囲んで、オブジェクトに代入するだけですから、簡単です。この関数を実行するには、最初aの値を()の中に入力して、
drift1(10)
とします。
↑を使って何回も実行してみると、結果がいろいろ変わるのがわかりますよね。
**３．結果をグラフ表示。前につかったplot()関数 [#xdaae7e5]
では、せっかく結果が出たのだから、グラフにしてみましょう。でも、その前に、以前使った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()
です。
results=append(results, a/20)
と書けば、resultsというベクトルにa/20の値が追加されます。但し、ひとつだけ重要な注意があります。
resultsというベクトルはあらかじめ存在していなければ、append()は使えない
つまり、for文を使って繰り返しを行う前に、
results=c()
という命令で、空ベクトルを作っておく必要があるのです。
では、まず、結果をresults()に入れてベクトルで返すような変更を、上のプログラムに加えて見ましょう
drift2=function(a){
results=c() #空ベクトルを作る
for(i in 1:30){
counta=0
for(j in 1:20){
if ( runif(1) < a/20 ){
counta=counta+1
}
}
results=append(results, a/20) #ここにあった out=counta/20 の代わりに、append()関数でベクトル作成
a=counta
}
print(results) #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){
results=c() #空ベクトルを作る
for(i in 1:30){
counta=0
for(j in 1:20){
if ( runif(1) < a/20 ){
counta=counta+1
}
}
results=append(results, a/20)
a=counta
}
plot(results) #ここにあった print文の代わりに、plot()でグラフ表示
}
グラフがすぐに表示できました！でも、点々では見にくいので、線で結んでみましょう。plot()のところを、
plot(results, type="l")
とかえるだけです。
drift4=function(a){
results=c() #空ベクトルを作る
for(i in 1:30){
counta=0
for(j in 1:20){
if ( runif(1) < a/20 ){
counta=counta+1
}
}
results=append(results, a/20)
a=counta
}
plot(results, type="l") #ここにあった print文の代わりに、plot()でグラフ表示
}
**４．何回もの計算結果を、一つのグラフにまとめて表示したい [#qa15e0db]
上で作ったプログラムを何度も走らせてみると、グラフの形がどんどん変わります。これはこれで面白いのですが、グラフの軸が一定していないし、１つのグラフに結果が重ねて表示されないので、ちょっと感じがつかめませんね。
そこで、いっそのこと、この実験を20回ぐらい繰り返してみて、結果を１つのグラフに表示させることにしましょう。同じことを何度も繰り返すのですから、for文をもう一つ、外側に作ればいいって、容易に想像がつきますよね？
drift5=function(a){
for(i in 1:20){ #これまでのプログラムの外側に新しく作ったfor文。20回繰り返す。
results=c()
for(j in 1:30){
counta=0
for(k in 1:20){
if ( runif(1) < a/20 ){
counta=counta+1
}
}
results=append(results, a/20)
a=counta
}
plot(results, type="l")
} #新しく作ったfor文の終わり。ここまで繰り返し。
}
さあ、走らせてみましょう。
drift5(10)
あれ？なんか変ですね。最後は0.0で線が一本でるだけです。
なぜだか分かりますか？
こたえは、aの値です。最初、aは10で与えられましたが、１集団分の繰り返しが終わった後、次の繰り返しの前に、下の数字に戻すのを忘れていました。
drift5=function(a){
firsta=a #関数に渡された（つまり()の中に書かれた）aの値をfirstaというオブジェクトに保存
for(i in 1:20){ #これまでのプログラムの外側に新しく作ったfor文。20回繰り返す。
a=firsta #繰り返しの最初に、aの値を初期値(firsta)に戻す
results=c()
for(j in 1:30){
counta=0
for(k in 1:20){
if ( runif(1) < a/20 ){
counta=counta+1
}
}
results=append(results, a/20)
a=counta
}
plot(results, type="l")
} #新しく作ったfor文の終わり。ここまで繰り返し。
}
できました。実行してみましょう。ウーン。。。うまくいってるようなのですが、１回ずつ、グラフが変わっちゃいます。グラフを重ねて描くには、par(new=T) というコマンドをplot()の前に使います。
drift5=function(a){
firsta=a #関数に渡された（つまり()の中に書かれた）aの値をfirstaというオブジェクトに保存
for(i in 1:20){ #これまでのプログラムの外側に新しく作ったfor文。20回繰り返す。
a=firsta #繰り返しの最初に、aの値を初期値(firsta)に戻す
results=c()
for(j in 1:30){
counta=0
for(k in 1:20){
if ( runif(1) < a/20 ){
counta=counta+1
}
}
results=append(results, a/20)
a=counta
}
par(new=T) #複数のグラフを重ね合わせる
plot(results, type="l")
} #新しく作ったfor文の終わり。ここまで繰り返し。
}
あれれ？そうか。始点の軸が固定されていないから、ずれちゃってますね。plot()関数にはylimというオプションがあって、y軸の目盛りを0から1に固定するには、次のように書きます。
drift5=function(a){
firsta=a #関数に渡された（つまり()の中に書かれた）aの値をfirstaというオブジェクトに保存
for(i in 1:20){ #これまでのプログラムの外側に新しく作ったfor文。20回繰り返す。
a=firsta #繰り返しの最初に、aの値を初期値(firsta)に戻す
results=c()
for(j in 1:30){
counta=0
for(k in 1:20){
if ( runif(1) < a/20 ){
counta=counta+1
}
}
results=append(results, a/20)
a=counta
}
par(new=T) #複数のグラフを重ね合わせる
plot(results, type="l", ylim=c(0,1)) #y軸の目盛りを0から1に固定する
} #新しく作ったfor文の終わり。ここまで繰り返し。
}
いよいよ完成にちかづいてきました。図が重なって見にくくなったら、
>frame()
か
>plot.new()
で図を新しく書き直して下さい。
**５．もっと大きな集団サイズで解析してみたい。 [#u93d3778]
何回かグラフを描かしてみると、わりとすぐに固定してしまうことがわかりますね。じゃあ、もっと、集団サイズがもっと大きかったらどうなるか気になってくるでしょ？
どうすればよいか？プログラムの中で、集団サイズを20に指定しているところを100とかに変更してもいいのですが、どうせなら、関数に値を渡すときに、指定できれば、つまり、
drift(繰り返し数, 世代数, 集団サイズ, 対立遺伝子の数）
とか指定して実行できれば、好きな大きさの集団で、好きな数だけ、シミュレーションを繰り返せますよね。さっきのプログラムで、数値で指定してあったところを変数に置き換えるだけですから、そんなに難しくなさそうです。
では、やってみましょう
drift6=function(x, y, z, a){
plot.new() #ついでに、図を新しく書き直すという関数をここに追加しました
firsta=a
for(i in 1:x){ #上のプログラムでは20でしたが、ここでは、x回繰り返す。
a=firsta
results=c()
for(j in 1:y){ #上のプログラムでは30だったけど、ここではy世代繰り返す。
counta=0
for(k in 1:z){ #上のプログラムでは20だったけど、ここでは集団サイズはz。
if ( runif(1) < a/z ){ #上のプログラムでは20だったけど、ここでは集団サイズはz。
counta=counta+1
}
}
results=append(results, a/z) #上のプログラムでは20だったけど、ここでは集団サイズはz。
a=counta
}
par(new=T)
plot(results, type="l", ylim=c(0,1))
}
}
できちゃいました。
あんまり大きい数字を入れると計算にすごく時間がかかります。途中で終了したいときは、 ESC キーを押してください。
**演習： [#v54b4118]
下のプログラムは遺伝的浮動をシミュレートする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
実行するための命令文と、結果のグラフを、1, 2のいずれの場合についても、課題提出ページに貼り付けなさい（注：1番目のシミュレーションの実行には、数分かかるかもしれません）
-問3: 上の２つのグラフを比較して分かったことを述べなさい。
-------------------------------
*時間があったら説明： いろんなシミュレーション [#he54915e]
**７．丸耳の表現型(aaで示されるもの）の集団内での表現型頻度をグラフで表示 [#o15e8af8]
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)を集団から排除してみよう [#zea292e6]
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の頻度を青線で表示する
}
}
**９．固定までにかかった世代数を求める関数 [#r6b6ba61]
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)
}
**１０.有限集団におけるヘテロ接合体頻度の変化 [#g5da17c3]
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")
}
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")