河田研のサイトにRで行う、遺伝的浮動のシミュレーションが載ってたのでやってみました。
###ちなみに
私は、Rstudioで動かしてます。
別にどっちでも問題ないと思いますが、個人的には楽なのでRstudioが好きです。
なんでもかんでもAnacondaから起動させちゃうマン
##スクリプトファイルの作成
まず、シミュレーションを実行するためのスクリプトファイルを作ります。
サイトにはコピペしろって書いてありますが、コピペしてもエラーでます。
ので、以下に書き換えました。
drift <- function(N,iFreq,Gen,Rep){
Frequency<-numeric(Gen+1)
Num<-numeric(Gen+1)
Generation<-numeric(Gen)
par(new=F)
for(j in 1:Rep){
Frequency[1]<-iFreq
for(i in 1:Gen){
Generation[i]<-i
Num[i]<-rbinom(1,N,Frequency[i])
Frequency [i+1]<-Num[i]/N
}
Generation[i+1]<-i+1
plot(Generation,Frequency,ylim=c(0,1),xlim=c(1,Gen),type="l",lty=1)
par(new=T)
}
}
スクリプトはこれでおっけーです。
Source on Save で保存します。
名前はなんでもいいですが、元サイトに則り'drift.R'とします。
##ファイルの実行
Rstudioの作業環境です。
右上のSource→を押します。
もしくは
source("drift.R")
でファイル内のコマンドを実行します。
##シミュレーション
シミュレーションの実行です。
こうやります。
drift(10, 0.5,50,1) #個体数(10)遺伝子頻度(0.5)世代(50)実行回数(1)
上記の結果
実行回数(100)の場合の結果
個体数(10)遺伝子頻度(0.5)世代(100)実行回数(1000)の場合
##更新:色分け
@StrawBerryMoon様に色分けを提案していただきました!
楽しさ倍増ですね。
drift <- function(N,iFreq,Gen,Rep){
Frequency<-numeric(Gen+1)
Num<-numeric(Gen+1)
Generation<-numeric(Gen)
par(new=F)
for(j in 1:Rep){
Frequency[1]<-iFreq
for(i in 1:Gen){
Generation[i]<-i
Num[i]<-rbinom(1,N,Frequency[i])
Frequency [i+1]<-Num[i]/N
}
Generation[i+1]<-i+1
plot(Generation, Frequency, ylim=c(0, 1), xlim=c(1, Gen), type="l", lty=1, col=j)
par(new=T)
}
}
```
```
drift(10, 0.5,50,100)
```
![スクリーンショット 2021-09-22 0.49.49.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1572428/0165af9a-3ce0-50f4-4b8c-8c01ea9cd6c5.png)
```
drift(100, 0.5,1000,100)
```
![スクリーンショット 2021-09-22 0.50.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1572428/08a3d489-3ee3-1f3a-027d-7ba3ef36f8f4.png)
##感想
楽しい!
条件を変えて色々試すとめちゃくちゃ楽しいです。
数が多いほどワクワクしました。
みなさんもぜひ
##宣伝
こういうのやってます。
https://flcl-browser.studio.site/page
学生のみなさん気軽に遊びに来てください。
ハイレベルな研究者(私は除く)がLTしてます。
Gatherにてアフタートークも設けているので、
他分野の研究者と交流を深めたり、研究の悩みを共有しやすいと思います。
「オンサイトの学会が激減した今こそ、オンラインで輪を広げよう」みたいな企画なので、沢山の参加者がいるとうれしいです。
たたた