はじめに
新型コロナウィルスに結合するDNAアプタマーを取得しようと思って実験しています。DNAアプタマーとは、特定の何かに強く結合するDNAのことで、基本的には一本鎖DNAです。結合する配列が得られれば、ウィルスの検出や治療など、様々な方面で利用できるかもしれません。
一本鎖DNAは、化学的に合成することができ、合成受託会社から買うことができます。購入時、DNA配列の中央部分をランダムな配列にしておきます。A、C、G、またはTをNで表すことにすると、例えばランダムな20塩基は、N20と表せます。ランダム配列は、20塩基程度の2つの別の配列で挟み込むようにしておきます。これら2つの配列は後にPCR増幅する時に必要です。ランダムな配列をN20とすると、配列としては全部で4の20乗つまりおおよそ10の12乗(2の10乗は1024、つまりおよそ10の3乗)の種類(=1兆種類)があることになります。色々な配列が混ざり合ったものをDNAプールと呼びます。
現在行なっているのは、新型コロナウィルスのスパイクタンパク質に結合するDNAアプタマーの取得です。スパイクタンパク質(こちらも買うことができます)を支持体に結合させ、これと中央部分がランダムである一本鎖DNAを混ぜ合わせます。結合していないDNAを洗い流したら、結合していたDNAを回収します(色々なやり方があります)。回収したDNAはPCRで増幅します。PCR増幅後は二本鎖DNAになっているので、なんらかのやり方で一本鎖DNAのみを精製します(こちらも色々なやり方が開発されています)。このサイクルを繰り返すことで、よく結合するDNA配列が濃縮されてくる、というのがDNAアプタマーの取得の方法です。
N20のDNAを出発材料として実際にやってみましたが、なかなか濃縮が進みません。そこでアプタマー取得についての理解を深めるため、pythonを使ってシュミレーションをしてみることにしました。
解離定数 Kdの復習
タンパク質とDNAの結合がどれくらい強いかを表す指標として解離定数というものがあります。
解離定数は小さいほど、結合が強いことを表します。
結合するタンパク質とDNAは、以下の式で表される平衡状態になると考えられます。
[DNA] + [Protein] ⇆ [DNA・Protein]
カッコ([])は濃度であることを示し、[DNA]はfreeなDNAの濃度、[Protein]はfreeなProteinの濃度、[DNA・Protein]はDNAとProtein複合体の濃度を表します。
ここで、左辺から右辺へ右向きの矢印があります。freeのDNAとfreeのProteinが結合して、複合体を作る反応は、[DNA]と[Protein]に比例し、反応速度定数K1を置けば、K1[DNA][Protein]となります。DNAとProteinの複合体が解離する速度は、[DNA・Protein]に比例し、反応速度定数K2を置けば、K2[DNA・Protein]となります。この速度が釣り合っている時、平衡状態にあるといいます。平衡状態では、
K1[DNA][Protein] = K2[DNA・Protein]
であり、この時の
\frac{[DNA][Protein]}{[DNA・Protein]} = \frac{K2}{K1} = Kd
を解離定数としてKdと表します。
この場合、Kd値の単位は、上の式が表す通り「濃度」です。
解離定数はしっかりとDNAに結合するタンパク質で1 nM程度、そうでもないもので100 nM程度かと思います。
全タンパク質濃度、全DNA濃度、解離定数から複合体濃度を計算する
全DNA濃度を$a$、全タンパク質濃度を$b$、解離定数を$c$とします。複合体濃度$[DNA・Protein]$を$a$、$b$、$c$で表します。
[DNA] + [DNA・Protein] = a
[Protein] + [DNA・Protein] = b
[Protein][DNA] = c[DNA・Protein]
これを、$[DNA・Protein]$について解くと
[DNA・Protein] = \frac{a+b+c\pm\sqrt{(a+b+c)^2-4ab}}{2}
となります。ルートの前に$\pm$がありますが、DNA濃度(あるいはタンパク質濃度)が0の時に、複合体濃度は0になりますので、この$\pm$はマイナスです。
pythonでこれを返す関数を作ります。詳細は別記事にありますが、普通にコードを書くとDNA濃度がそれほど低くないのに、計算値が丸められてしまって返り値(複合体濃度)が0になってしまうと言う問題がありました。Decimalクラスを使えば大丈夫なようです。なお全タンパク質濃度は既知であり、全DNA濃度は実験時に容易に測定可能な値です。
from decimal import *
getcontext().prec = 100
def complex_conc(a,b,c):
#aは全DNA濃度、bは全タンパク質濃度、cはKd値
sum = Decimal(a + b + c)
v = (sum - (sum**2 - Decimal(4*a*b))**Decimal('0.5'))/2
return float(v)
ランダムなDNAプールのKd値分布
10の12乗の種類のDNA配列プールがあるとして、仮にその全ての種類についてKd値を求めることができたとしましょう。この時、Kd値はどのように分布しているでしょうか？シミュレーションをするにあたって、この分布について仮定する必要がありそうです。
大多数のDNA配列は、ほとんど標的に結合しないはずなので、Kd値が極めて大きいものが大多数を占めているとおもいます。kd値が10 nM程度のものが取れれば万々歳ですが、ここでは、含まれているDNAのうち、一番結合力が強いものの解離定数が1 nMであるとし、解離定数が大きいものは解離定数が大きいほど指数関数的に多い、と仮定します。
実際には1つ1つのDNA種ごとにKd値は異なるはずですが、計算コストを減らすため、Kd値の階級分布を考えることにします。階級値を1, 10, 20, 30, ...100, 200, 300,....1000, 2000, 3000,....とすることにし、各階級に初期濃度を割り振ります。計算時には、各階級に含まれるDNA種すべてが、所属する階級の階級値をKd値として持つとして計算することにします。
以下のコードでは、(特に理由はないですが)Kdの最大の階級値は90,000,000となっています。
#Kd値の階級値を含むndarrayを作成。1,2,3,..10,20,30,...100,200,300....,70000000,80000000,90000000とする。
list = []
for j in range(8):
for i in range(1,10):
list.append(i*10**j)
kd = np.array(list)
続いて、各階級の初期濃度に相当する値を含むndarrayを返す関数initialConcentrations()を準備します。この関数では、Kd値が増えるに従って指数関数的に濃度が増えるように、また、DNA濃度の総和が指定した値になるようにします。また、後のシミュレーションのために、一定のKd値以下の濃度を0にできるよう、cut off値を設定できるようにします。
まず、指数関数的になるように各階級に仮の濃度を与え、これら仮濃度の総和をとって、想定する実際のDNAの全濃度との比を、仮濃度に掛け合わせて、濃度にします。仮の濃度は定数をkとして、以下の式で計算します。
tentativeConcentration = k^{Kd}
#初期値をndarrayで返す関数。
#kd_listは、kd値のリスト、dはDNA濃度の総和、kは定数、cutoffは、この値以下のKd値のDNAの濃度を0にするために指定
def initialConcentrations(kd_list,d,k,cutoff):
tc = np.array([k**n for n in kd])#仮の濃度
for n in range(len(kd_list)):#ある程度低いKd値のものが存在しないように(濃度を0に)することができるようにする。
if (kd_list[n] <= cutoff):
tc[n] = 0
return tc*d/sum(tc)#総和が全DNAの濃度となるように変換
定数kはどの範囲が現実的か考察します。以下にk = 1.0000001からk = 1.0000005で計算した様子をグラフにしてみます。DNAの全濃度は10 nMとします(1 pmol / 100 µl)。
import matplotlib.pyplot as plt
d_total = 10
#各階級値の濃度を決める。
concentrations1 = initialConcentrations(kd,d_total,1.0000001,0)
concentrations2 = initialConcentrations(kd,d_total,1.0000002,0)
concentrations3 = initialConcentrations(kd,d_total,1.0000003,0)
concentrations4 = initialConcentrations(kd,d_total,1.0000004,0)
concentrations5 = initialConcentrations(kd,d_total,1.0000005,0)
fig, ax = plt.subplots()
ax.plot(kd, concentrations1,'.',label="1.0000001")
ax.plot(kd, concentrations2,'.',label="1.0000002")
ax.plot(kd, concentrations3,'.',label="1.0000003")
ax.plot(kd, concentrations4,'.',label="1.0000004")
ax.plot(kd, concentrations5,'.',label="1.0000005")
ax.ticklabel_format(style='plain',axis='x',scilimits=(0,0))
ax.set_yscale('log')
plt.legend()
plt.xlabel("Kd")
plt.ylabel("conc (nM)")
plt.show()
こんなグラフになりました。縦軸は対数目盛にしてあります。
kが大きいほど厳しい見積り(=初期DNAプールに良いものがすくない)です。
ここで5通りのkの値がおおよそ良いか考えてみます。仮に100 µlに1分子あるとすると(=1リットルに10000分子)、 アボガドロ数$6\times10^{23}$を使って濃度を計算すると$1.67\times10^{-10} nM$となりますので、k=1.0000001から1.0000005で大体良い範囲を押さえられているように思います。以降k=1.0000003で実行します。
結合したDNAの回収に相当するメソッド
アプタマー取得実験の各サイクルでは、標的に結合したDNA分子を回収してきます。このとき、各階級ごとにKd値に応じて一定の割合が回収されてくることになります。Kd値が低い階級では、高い効率で回収されてくる一方で、Kd値が高い階級では、低い効率で回収されてきます。これを表すメソッドbinding()を作成します。DNA濃度のリスト、Kd値のリスト、全タンパク質濃度を渡すと、支持体に結合し回収されたDNAの濃度が返されます。内部で使用するメソッドcomplex_conc()については冒頭で作成した通りです。
def binding(d,p,k):#dはDNA濃度のリスト、kはkd値のリスト、pは全タンパク濃度
r = []
for i in range(k.size):
r.append(complex_conc(d[i],p,k[i] ))
return np.array(r)
DNAのPCR増幅に相当する処理
回収されてきたDNAを鋳型にPCR増幅を行います。PCRでは偏りなく増幅されるとし、また反応スケールなどを調整すれば、必要量を十分合成できます。これを10 nMとかいった決まった濃度に調製しますので、次のサイクルに用いる各kd階級のDNA濃度は以下のように更新できます。
result = binding(concentrations3,p_total,kd)
concentrations3 = result * d_total / result.sum()#PCR増幅で元の濃度に戻すことに相当。
シミュレーション
購入したタンパク質を、メーカーが推奨する濃度(なぜかはわからないが600 ng/µl)にすると分子量から計算しておおよそ25 pmol/µlとなり、これを100 µlの反応系で1 µl使用していましたのでタンパク質濃度は250 nMとしました。5サイクル回すことにして、また、どれくらいのDNAが各サイクルで回収されてきたかをプロットすることにすると、コードは以下のようになります。
import numpy as np
import matplotlib.pyplot as plt
from decimal import *
getcontext().prec = 100
#タンパク質の濃度 25 pmol / 100 µlなら 250 nM
p_total = 250
#全DNAの濃度。nM
d_total = 10 #100 µlで1 pmolこれは、6.0E+11分子。
#Kd値の階級値を決める。1,2,3,..10,20,30,...100,200,300....とする。
list = []
for j in range(8):
for i in range(1,10):
list.append(i*10**j)
kd = np.array(list)
#初期濃度をndarrayで返す関数。
#kd_listは、kd値のリスト、dはDNA濃度の総和、kは定数、cutoffは、この値以下のKd値のDNAの濃度を0にするために指定
def initialConcentrations(kd_list,d,k,cutoff):
tc = np.array([k**n for n in kd])#仮の濃度
for n in range(len(kd_list)):#ある程度低いKd値のものが存在しないときを表現できるようにする。
if (kd_list[n] <= cutoff):
tc[n] = 0
return tc*d/sum(tc)#総和が全DNAの濃度となるように変換
#DNA濃度とタンパク質濃度とkd値を指定すると、形成された複合体の濃度を返す
def complex_conc(a,b,c):
#aは全DNA濃度、bは全タンパク質濃度、cはKd値
sum = Decimal(a + b + c)
v = (sum - (sum**2 - Decimal(4*a*b))**Decimal('0.5'))/2
return float(v)
def binding(d,p,k):#dはdna濃度のリスト、1kはkd値のリスト、pは全タンパク濃度
r = []
for i in range(k.size):
r.append( complex_conc(d[i],p,k[i] ))
return np.array(r)
#各階級値の濃度を決める。
k=1.0000003
concentrations = initialConcentrations(kd,d_total,k,0)
fig, ax = plt.subplots()
ax.plot(kd, concentrations,'.',label=str(k))
ax.ticklabel_format(style='plain',axis='x',scilimits=(0,0))
ax.set_yscale('log')
plt.legend()
plt.xlabel("Kd")
plt.ylabel("conc (nM)")
plt.show()
#初期値をプロット
fig, ax = plt.subplots()
ax.set_xscale('log')
ax.plot(kd, concentrations,label="start")
#シミュレーション開始
recoveryRate=[]
for i in range(20):
result = binding(concentrations,p_total,kd)#結合したものの濃度に変換
recoveryRate.append(result.sum()/d_total)
concentrations = result * d_total / result.sum()#PCR増幅で元の濃度に戻すことに相当。
plt.plot(kd,concentrations,label="round_"+str(i+1))
plt.legend()
plt.xlabel("Kd")
plt.ylabel("conc")
plt.show()
plt.plot(recoveryRate)
plt.xlabel("cycles")
plt.ylabel("recovery rate")
plt.show()
実行結果
各サイクルで、inputのDNA濃度に対してどれくらいのDNAが回収されたかの割合です。
わずか5サイクルで良い感じに濃縮されています。それでもkd値が1 nMから数百nM程度までが幅広く存在していることがわかります。サイクル数を20までにして、kd値が低い方20階級分をプロットしてみます。
3サイクル目以降、Kd値が低いものの増加が緩やかになってきています。タンパク質濃度を下げるとどうでしょうか？250 nMから25 nMに下げてみます。
DNAプールに良いDNA配列は含まれているのだろうか？
上で示したように濃縮が進むのが理想ですが、実際はそうではないのでシミュレーションを始めたのでした。DNAプール中に、強い結合能を持つような(=Kd値が十分に低い)ものが存在しない場合はどうなるのか？仮に一番良いKd値が1000だったら、どのように濃縮が進むでしょうか？
サイクルを経るにつれてkd値が低いものの割合が増加するように集団構成が変化していくようです。
この時の回収率をプロットすると以下のようになります。
回収率に着目すると、DNAプールの中に果たして良いものがあるのかどうか、目安になるようです。
Kd値が10000以上のものしか含まれない場合:
シミュレーションをやってみた結果
タンパク質を結合させた担体への結合など、ここでのシミュレーションで考慮されていないことがたくさんありますが、
1. DNAの集団構成は、サイクルを経るにつれて、Kd値が低いものが増えるように変化していく。
2. タンパク質濃度を下げると、濃縮がより速くなる。
3. Kd値が低いものがない場合でも、濃縮は進む。ただし回収率は低いところで頭打ちになる。
4. DNAプールの中に良いものがあるかどうかは、回収率で判断できる。
ということについて、より理解が深まったと思います。python便利だな。
終わりに
N20のプールでは回収率が低いところで頭打ちになりました。このプールには良いものは無いと判断して、N50のプールを調べてみることにしました。こちらも進めていますが難航中です。washのやり方と、支持体へDNAが結合してしまうのが問題のようですが、機会があればまた別記事をあげたいと思います。