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

河田研のサイトに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)

drift(100, 0.5,1000,100)

感想

楽しい！
条件を変えて色々試すとめちゃくちゃ楽しいです。
数が多いほどワクワクしました。

みなさんもぜひ

宣伝

こういうのやってます。

学生のみなさん気軽に遊びに来てください。
ハイレベルな研究者（私は除く）がLTしてます。

Gatherにてアフタートークも設けているので、
他分野の研究者と交流を深めたり、研究の悩みを共有しやすいと思います。
「オンサイトの学会が激減した今こそ、オンラインで輪を広げよう」みたいな企画なので、沢山の参加者がいるとうれしいです。

たたた