Introduction to population genetics simulation

ここでは集団遺伝学のイメージを掴むために簡単なシミュレーションを行います。以下ではpythonを使ってシミュレーションを行います。他のプログラミング言語でも可能ですので、すでに他の言語に親しんでいる人はお好きな言語で実装していただいて構いません。

遺伝的浮動の影響をシミュレートして、その結果を図示するところまでが目標です

おまじない

標準的なpythonに加えて、いくつか機能を追加します。今回は以下のライブラリを使います。

In [1]:
# time
import time

# algebra
import numpy as np
# plot
import matplotlib.pyplot as plt
# make your plot more attractive 
import seaborn

遺伝子頻度の確率的変動

単純化された状況を考えます。ゲノム上のある一箇所の遺伝子座に注目します。その遺伝子座に祖先型と変異型の2タイプが存在する場合を想定します(1-locus 2-allele model)。その場所には祖先型と変異型のどちらかしか存在しません。つまり「変異型のさらに変異型」などは考えません。この仮定のもとで変異型の個数が変化するプロセスをコンピューター上で再現します。

伝統的な集団遺伝学(教科書など)では、二倍体生物の集団の大きさを$N$として、集団全体では$2N$のゲノムが存在する状態を考えるケースが多いです。半数体生物を想定している場合は集団サイズが$2N$と考えれば、結果をそのまま使えます。ここでは二倍体生物特有の現象(優性劣性のある自然選択など)は考えないこととします。以下では集団遺伝の慣習に従って$N$を個体数、集団中のゲノムの数を$2N$で表現します

いま集団に存在する$2N$個のゲノムのうち$x$個が変異型だったとします。この時の割合$p=\frac{x}{2N}$を遺伝子頻度と呼びます。今回はこの割合$p$の時間変化をシミュレートします。

以下のシミュレーションでは祖先型であっても変異型であっても生存や生殖成功率に関して何の影響もないものと仮定します(中立モデル)。この場合、子世代の各個体の持つ遺伝子型は親世代集団からランダムに選ばれたタイプ同士が組みになったものと考えることができます。このように遺伝子の頻度に基づいて次世代の遺伝子の組成を決めるモデルをWright-Fisher modelといい、集団遺伝で標準的に用いられているモデルの一つです。

第0世代

いま仮に集団の大きさ(集団中のゲノム数)が100だとしましょう。二倍体の個体数だと50個体で100ゲノム、半数体の生物だと100個体で100ゲノムに相当します。環境は常に一定であり、集団の大きさは変わらないものとします。説明のため、第0世代目の生物集団ではすでに祖先型と変異型が混ざっている状態にあると考えます。第0世代目の変異型の数が40である場合について考えてみましょう。

この初期状態をプログラム中の変数で表します。集団の大きさを変数sizeであらわし、変異型の数をinitで表すものとすると、初期状態は以下のように設定できます

In [2]:
size=100
init=40

これが親世代です。この親世代から配偶子が作られ、遺伝子プールが形成されます。

では親世代の変異型の遺伝子頻度を求めましょう

In [3]:
p=init/size
print('Gene frequency (0th gen): p =',p)
Gene frequency (0th gen): p = 0.4

この親世代の遺伝子プールから集団サイズ分だけ遺伝子をランダムに抽出して、第一世代を作ります。

第1世代

この親世代の遺伝子プールからランダムに選ぶという操作をプログラム上では

  1. 一様乱数を発生する
  2. 一様乱数が変異型の遺伝子頻度よりも小さい時は変異型を選ぶ。大きい時は祖先型を選ぶ

という方法で実現します。

次世代の集団を作るプロセスをプログラムの流れで表現すると、

  1. 子世代の変異型の個数をnextで表す
  2. 最初は子世代の変異型の数を0にする、next=0
  3. 集団のサイズに達するまで(size回)、遺伝子のタイプを決定する
    1. 一様乱数を発生させる
    2. 親世代の変異型の頻度$p$と乱数の値を比較
    3. 乱数が$p$よりも小さければ変異型、そうでなければ祖先型とみなす
  4. 子世代の全てのタイプが決定し終わった時のnextの値が、子世代の変異型の数になる

このプロセスをプログラムで表現すると、

In [4]:
next=0          # initialize num of mutant alleles in the next generation

# pick up 'size' genomes 
for i in range(size):
    
    # determine allele type
    if( np.random.rand()<p ):
        
        # if derived
        next+=1

# print out results
print('1st generation')
print(size-next) # ancestral type
print(next)      # derived type
1st generation
59
41

第2世代

第二世代目の集団にとっては第一世代が親ですので、第0世代の頻度ではなく、第一世代目の遺伝子頻度に基づいて第二世代を作ります。第0世代目の遺伝子頻度とは同じこともあるし微妙に変化していることもあります。この場合は、

In [5]:
p=next/size
print('Gene frequency (1st gen): p =',p)
Gene frequency (1st gen): p = 0.41

親世代の頻度を反映した遺伝子プールからランダムにサンプルし、第二世代目を作ります

In [6]:
next=0
for i in range(size):
    if( np.random.rand()<p ):
        next+=1

# print out
print('2nd generation')
print(size-next) # ancestral type
print(next)      # derived type
2nd generation
58
42

第3世代

第三世代は第二世代から作られるので、第二世代目の遺伝子頻度が影響します

In [7]:
p=next/size
print('Gene frequency (2nd gen): p =',p)

next=0
for i in range(size):
    if( np.random.rand()<p ):
        next+=1

print('3rd generation')
print(size-next) # ancestral type
print(next)      # derived type
Gene frequency (2nd gen): p = 0.42
3rd generation
53
47

第4世代

In [8]:
p=next/size
print('Gene frequency (3rd gen): p =',p)

next=0
for i in range(size):
    if( np.random.rand()<p ):
        next+=1

print('4th generation')
print(size-next) # ancestral type
print(next)      # derived type
Gene frequency (3rd gen): p = 0.47
4th generation
54
46

第5世代

In [9]:
p=next/size
print('Gene frequency (4th gen): p =',p)

next=0
for i in range(size):
    if( np.random.rand()<p ):
        next+=1

print('5th generation')
print(size-next) # ancestral type
print(next)      # derived type
Gene frequency (4th gen): p = 0.46
5th generation
60
40

関数を定義する

ここまでの計算を振り返ると、各世代で同じ操作を繰り返していることがわかります。つまり、

  1. 前の世代の遺伝子頻度を計算する
  2. 子世代の各遺伝子のタイプを、親世代の遺伝子プールからランダムにサンプリングすることで決定する

プロセスそのものは毎世代同一です。違うのは親世代の遺伝子頻度とその結果得られる子世代の変異型の数です。

毎回同じプロセスを繰り返す部分を「関数」として定義して独立させます。プログラムを再利用してコーディングの手間を省き、エラーを防ぎます。

In [10]:
def make_next_gen(derived, size=100):
    '''
    Calculate the number of derived alleles in the next generation.
    
    Arguments:
        derived: number of derived alleles in the parent population
        size: size of the population
    Returns:
        the number of derived alleles in the child population
    '''
    p=derived/size
    next=0
    for i in range(size):
        if( np.random.rand()<p ):
            next+=1
    return next

シミュレーション

 この関数を使うと、先ほどのシミュレーションは少しすっきりします。第5世代目までの変異型の数をシミュレートしてみましょう。

In [11]:
init = 40
print(init)

next=make_next_gen(init, 100)
print(next)

next=make_next_gen(next, 100)
print(next)

next=make_next_gen(next, 100)
print(next)

next=make_next_gen(next, 100)
print(next)

next=make_next_gen(next, 100)
print(next)
40
42
35
35
31
35

こうやってみてみると、毎世代

  1. 親の変異型の数を与えて
  2. 子の変異型の数を受け取る

と、同じ操作を繰り返していることがわかります。

単純な繰り返し操作はコンピュータに任せましょう

In [12]:
# initialization
derived = 40
size = 100

# print out the initial number of derived alleles
print(derived)

# for 20 generations
for i in range(20):
    
    # make the next generation
    derived = make_next_gen(derived, size)
    
    # print out the number of derived alleles
    print(derived)
40
43
46
44
42
41
44
41
37
34
40
44
39
33
28
22
22
26
22
16
15

遺伝子数ではなく、遺伝子頻度でプリントアウトしましょう。

In [13]:
# initialization
derived = 40
size = 100

# print out the initial number of derived alleles
print(derived/size)

# for 20 generations
for i in range(20):
    
    # make the next generation
    derived = make_next_gen(derived, size)
    
    # print out the number of derived alleles
    print(derived/size)
0.4
0.39
0.4
0.41
0.35
0.33
0.26
0.23
0.22
0.22
0.19
0.13
0.12
0.17
0.2
0.23
0.11
0.12
0.14
0.13
0.08

シミュレーションの中心部分は完成です

さらに手を加えましょう。これまでは結果を画面に出力して、それで終わっていました。結果を解析したりグラフを作るためにはデータを保存する必要があります。

ここでは遺伝子頻度をリスト型の変数に保存します。リスト型とは一連の数字の列を保存するためのデータ型です。

trという名前のリスト型変数を作ります。入れ物だけ作り、中身は空っぽです。その後、初期頻度を追加(append)しています。毎世代得られた変異型の遺伝子頻度を追加していきます。

In [14]:
# initialization
derived = 40
size = 100

# initialize a list
tr = []

# append the initial frequency
tr.append( derived/size )


# 50 generations
for i in range(50):
    
    # make next generation
    derived = make_next_gen(derived, size)
    
    # append freq
    tr.append( derived/size )
    
In [15]:
print(tr)
[0.4, 0.39, 0.4, 0.34, 0.44, 0.37, 0.36, 0.37, 0.37, 0.28, 0.28, 0.3, 0.31, 0.32, 0.35, 0.28, 0.3, 0.38, 0.45, 0.48, 0.47, 0.47, 0.39, 0.43, 0.41, 0.45, 0.44, 0.38, 0.36, 0.37, 0.37, 0.41, 0.42, 0.43, 0.48, 0.51, 0.43, 0.39, 0.37, 0.45, 0.45, 0.44, 0.45, 0.4, 0.43, 0.4, 0.37, 0.4, 0.46, 0.5, 0.49]

可視化

横軸に時間を取り縦軸に頻度をとってプロットすると、頻度の軌跡を描くことができます

In [16]:
fig=plt.figure()
ax=fig.add_subplot(1,1,1)

ax.plot(tr)

ax.set_xlim(0,50)
ax.set_ylim(0,1)
plt.show()

軌跡一つだけでなく、何本も描いてみよう。やり方はシンプルです。ここまでで行ったtrajectoryの作製を、好みの回数だけ繰り返すだけです。

繰り返して得られた軌跡を、順次simという名前のリストにappendしましょう

sim.append(tr)

つまりリストのリストを作成します

In [17]:
# a list of lists
sim = []

# set initial conditions
derived = 40
size = 100

howmany = 1000
generations = 1000

# simulation starts
t_start = time.time()

# repeat simulation 'howmany' times
for j in range(howmany):
    
    # reset initial condition 
    derived = 40
    size = 100
    tr = [derived/size]
    
    # simulate frequency trajectory
    for i in range(generations):
                
        # make next generation
        derived = make_next_gen(derived, size)
        # append freq
        tr.append( derived/size )
    
    # append trajectory
    sim.append(tr)

# finished
t_end = time.time()
print( '{0:.3f} seconds'.format(t_end-t_start))
62.363 seconds

ここで得られたsim[0]が1回目のシミュレーション結果です。sim[1]が2回目のシミュレーション結果。sim[2]が3回目の...

各sim[ i ]は数値の並びからできているリストで、0世代目からの毎世代ごとの頻度が保存されています。

この結果をプロットしましょう。以下に先頭から10回目までのシミュレーション結果をプロットした例を示します

In [18]:
fig=plt.figure()
ax=fig.add_subplot(111)
for i in range(10):
    ax.plot(sim[i])
ax.set_xlim(0,50)
plt.show()

解説

上でやっていることは、最初にfigure領域を作製し

fig=plt.figure()

プロット領域を作製し

ax=fig.add_subplot(1,1,1)

順次プロットを載せて

for i in range(10):
    ax.plot(sim[i])

x軸の範囲を設定し、

ax.set_xlim(0,50)

画面にグラフ出力

plt.show()

頻度変化、固定確率、固定まで 消失までの待ち時間分布

以上のように頻度変化を求めることができました。

せっかくシミュレーションしたので、固定確率(割合)・固定あるいは消失までの待ち時間の分布とその平均を求めてみましょう。固定するまで、消失するまでの待ち時間を数学を使って解析的に求めることは難しいけれども、シミュレーションを使えば比較的容易に求めることができます

In [19]:
fixed = []
lost = []

for i in range(len(sim)):
    trajectory = np.array(sim[i])
    if(trajectory[-1]>0):
        fixed.append( trajectory[ trajectory<1 ].shape[0] )
    else:
        lost.append( trajectory[ trajectory>0 ].shape[0] )
        

p_fix = len(fixed)/len(sim)
ave_fixation_time = np.average(fixed)
ave_time_until_lost = np.average(lost)

print( "the proportion of fixed trajectories is {0:.2f}".format(p_fix))
print( "the average time until fixation of the allele is {0:.2f}".format(ave_fixation_time))
print( "the average time until loss of the allele is {0:.2f}".format(ave_time_until_lost))
the proportion of fixed trajectories is 0.40
the average time until fixation of the allele is 147.38
the average time until loss of the allele is 119.41

解説

まずはシミュレーションを固定した場合と消失した場合に分ける。最初に配列を初期化して、

fixed = []
lost = []

シミュレーション結果を一つづつ取り出す。そのままでは扱いにくいのでnumpyのnp.array形式に変換する。変換したデータ(i番目のシミュレーション結果)をtrajectoryという変数に保存している

for i in range(len(sim)):
    trajectory = np.array(sim[i])

次に各頻度変化の最後の部分を見て、その時の頻度を確認する。十分長い時間シミュレーションを行なっているので、全ての場合は固定あるいは消失しており、中途半端に多様性として残っている場合はないと判断する。もし固定していれば頻度が1になっているはずで、消失した場合は頻度が0になっているはず。

そこで各シミュレーションの1番最後を取り出し、その値が0よりも大きい時は固定したものと判断する。ここで1と等しいならば、という判断をしないことに注意。コンピュータ内部では二進法で計算しており、桁の末端まで厳密に数学的な実数の1.になるかどうかは保証できないから。

    if(trajectory[-1]>0):

固定していた場合、頻度変化データのうち固定する前まで(頻度が1よりも小さいもの)を取り出す。

trajectory[ trajectory<1 ]

は、trajectoryというlistデータのうち、1よりも小さいものという意味。

一般にnp.arrayデータの一部だけを取り出すときは

trajectory[ 0:100 ]

と、後ろの[]に範囲を指定する。範囲の代わりに条件を入れると、その条件にあったものだけを取り出してくれる。

trajectory[ trajectory<1 ].shape

で取り出したデータの縦横の次元を表示してくれる。実際のデータの個数はshape[0]で取り出すことができる。これが頻度1になる前までの頻度データの個数、ここでは世代ごとに頻度データを計算しているので世代数、に相当する。よって

    if(trajectory[-1]>0):
        fixed.append( trajectory[ trajectory<1 ].shape[0] )

は固定するまでの世代数をfixedに付け足す、という意味になる。

plot

さて、頻度変化と固定するまでの時間の分布、消失までの時間の分布を図にしてみよう。いろいろと付加的な情報も付け加えてみました。

In [20]:
xmax=800
ymax=70

fig = plt.figure(figsize=(6,10),dpi=80)

ax1 = fig.add_subplot(3,1,1)
for i in range(10):
    ax1.plot(sim[i])
ax1.set_xlim(0,50)
ax1.set_ylabel("Allele frequency")


ax2 = fig.add_subplot(3,1,2)
ax2.hist(fixed, bins=50, normed=True)
ax2.set_xlim(0,xmax)
ax2.set_ylim(0,.01)
ax2.set_ylabel("Count")
ax2.text(500, .0060, "Time until fixation")
ax2.annotate("average ({0:.1f} generations)".format(ave_fixation_time), 
             xy=(ave_fixation_time, .0030), 
             xytext=(ave_fixation_time, .0050),
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))


ax3 = fig.add_subplot(3,1,3)
ax3.hist(lost, bins=50, normed=True)
ax3.set_xlim(0,xmax)
ax3.set_ylim(0,0.01)

ax3.set_xlabel("Generation")
ax3.set_ylabel("Count")
ax3.text(500,.0060, "Time until lost")
ax3.annotate("average ({0:.1f} generations)".format(ave_time_until_lost), 
             xy=(ave_time_until_lost, .0030), 
             xytext=(ave_time_until_lost, .0050),
             arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))



fig.tight_layout()
#fig.savefig("Drift_simulation_1.eps")
plt.show()

おまけ

集団サイズが減少した時

In [21]:
# a list of lists
sim = []

# set initial conditions
size_ans = 1000
derived_ans = size_ans*.5


howmany = 10
generations = 100

# simulation starts
t_start = time.time()

# repeat simulation 'howmany' times
for j in range(howmany):
    
    # reset initial condition 
    derived = derived_ans
    size = size_ans

    tr = [derived/size]
    
    # simulate frequency trajectory
    for i in range(generations):
        
        # population size change
        if i==20:
            # reduction
            size = int(size/10)
            derived = int(derived/10)
            
            # growth
            #size *= 10
            #derived *= 10
        
        # make next generation
        derived = make_next_gen(derived, size)
        # append freq
        tr.append( derived/size )
    
    # append trajectory
    sim.append(tr)

# finished
t_end = time.time()
print( '{0:.3f} seconds'.format(t_end-t_start))


# plot
fig=plt.figure()
ax=fig.add_subplot(111)
for i in range(10):
    ax.plot(sim[i])
ax.set_xlim(0,50)
plt.show()
0.165 seconds

ボトルネック

In [22]:
# a list of lists
sim = []

# set initial conditions
size_ans = 1000
derived_ans = size_ans*.5

howmany = 10
generations = 100

# simulation starts
t_start = time.time()

# repeat simulation 'howmany' times
for j in range(howmany):
    
    # reset initial condition 
    derived = derived_ans
    size = size_ans
    tr = [derived/size]
    
    # simulate frequency trajectory
    for i in range(generations):
        
        # population size change
        if i==20:
            # reduction
            size = int(size/10)
            derived = int(derived/10)
            
        if i==40:
            # growth
            size *= 10
            derived *= 10
        
        # make next generation
        derived = make_next_gen(derived, size)
        # append freq
        tr.append( derived/size )
    
    # append trajectory
    sim.append(tr)

# finished
t_end = time.time()
print( '{0:.3f} seconds'.format(t_end-t_start))


# plot
fig=plt.figure()
ax=fig.add_subplot(111)
for i in range(10):
    ax.plot(sim[i])
ax.set_xlim(0,100)
plt.show()
0.521 seconds

© 2017 Kosuke Teshima