Tm値について
プライマーやプローブを設計するにあたり、Tm値は重要な指標の一つで、A、B二つの鎖からなる二重鎖（DNA:DNA、RNA:RNAあるいはDNA:RNAの複合体）の50%がAとBとに乖離する温度である。古い文献ではTm50と書かれていたりする（melting temperature 50)、というのはよく知られています。
ところが、こういわれると少しハッとします。
「WebにTm計算ツールがあって、プライマーの配列を入力したらTm値が出て便利なんだけど、Tmって二つの鎖があって初めて計算できるものですよね、なんで片方しか入力してないのに計算できるんだ？これってもしかしてselfに対するTm？？？」
すごくいい疑問ですね、これ。実際には入力した配列に対する相補配列が裏で生成されて、それとのTm値が計算されているというのが実情だと思います。入力配列自身との相互作用もプライマーのデザイン上では重要ですので、上のような疑問を持つのはとても良いことだと思います。また、相補配列が自動で生成されるのは便利なのですが、ミスマッチがあるケースでのTm値が知りたい場合には、その便利さがかえって不自由ということになってしまいます。
”片方しか入力してないのに計算でき”てしまっていることに疑問を持たなかった人も、たいてい最初に書いたようなTm値の定義については知っています。しかしそのままでは知識と実務上の理解が乖離してしまっています。それでは想像力豊かにPCRなどの最適化やアッセイ開発ができません。Biopythonの提供するTm計算関数を使って溝をうめていきましょう。
後半にも書きますが、二つの鎖を与えてTm値を計算させることができると、mismatchのある場合のTm shiftを調べることができ、応用（と心配）範囲が広がります。この稿ではWeb上でのサービスや簡単ツールだけではできないことをやってみることが目的です。
Tm計算その１：基本編
Bio.SeqUtils.MeltingTempを使います。いくつかのmethodが提供されていますが、事実上のde facto standardである、Nearest Neighbor 法に基づくTm_NN()を使用しましょう。Documentationをみると、このmethodは以下のようなargumentをとることがわかります。
Bio.SeqUtils.MeltingTemp.Tm_NN(
seq, check=True, strict=True,
c_seq=None, shift=0, nn_table=None, tmm_table=None,
imm_table=None, de_table=None, dnac1=25, dnac2=25,
selfcomp=False, Na=50, K=0, Tris=0, Mg=0, dNTPs=0, saltcorr=5
)
たくさんありますね。よく見ると、seq 以外のargumentにはdefault値が設定されています。動かしてみるだけならば、primerやprobeの配列を一本与えるだけで使えます。
import Bio.SeqUtils.MeltingTemp as mt
print(mt.Tm_NN(seq='CGATATTACGCGGCAAACTGCTA')) # テキトーな配列で試す
# 55.886758538428694
動きました。return valueをそのまま見たかったので、formatも何もしませんでした。
最低でも配列が一つ必要なのは当然ですね。一本だけで計算ができるということは、内部的に相補配列が生成されているとみることができ（なぜならNN法では、連続する2塩基分の枠に入る塩基対情報がパラメターテーブル参照に必須だからです）、実際にTm_NN()のコードをみてみると、そうなっています。Argumentにc_seq(complement sequence)があって、default値はNoneになっていますが、これがパートナー配列を与える場合のargument keywordです。もうひとつ、selfcomp=Falseがあるのも覚えておくといいかもしれません。これをTrueとすれば自身に対するTm値が計算できます。
seq=に与える値の型ですが、DocumentationにはstringあるいはBiopythonのSeq objectをとる、と書いてあります。Bio.Seq objectでもいいということは覚えておくといいかもしれません。
seq argumentはただ一つのpositional argumentなので、最初のargumentとして配列を与えればkeyword seq=を省略することができますが、この後みていくように実務ではほかのargumentも必須になると思いますので、配列をfirst argumentにする癖をつけるにせよ、常にseq=を書くようにしたほうがいいと思います。
さて、 ほかのargumentをどう設定するか、ですが、いま使っているツールでの設定に従うのがスジです。研究室内、開発部内、先輩後輩と私、過去の俺と今の俺、などの間で、Tm計算するときこれね、というルールがあるはず。Tm計算にはいろいろな方法があるので、整合性をとるためにも必ず計算方法とパラメターのセットを決めておくべきです。これをおろそかにすると、あれ？Tmは60度前後にしてね、って言ったのに、全部低すぎるんだけどなんで？というようなことがオリゴDNAの注文後に起こったりしてしまいます。計算方法（ツール）とパラメターは、必ず設定、記録しておきましょう。
例として、別ツール上での設定が、
[Na+] = 50 mM、[Mg++]= 1.5 mM、[dNTPs] = 0.2 mM、[oligos] = 200 nM
だとします。これに対応するコードは、以下のようになります。
import Bio.SeqUtils.MeltingTemp as mt
est_tm = mt.Tm_NN(seq='CGATATTACGCGGCAAACTGCTA',
Na=50, Mg=1.5, dNTPs=0.2,
dnac1=100, dnac2=100)
print(est_tm)
# 63.97649655581523
もし、対応するパラメターがない、あるいは、いままで一つの値として与えていた[oligos] = 200 nMをdnac1、dnac2にどう設定するかといったことは、実際に試して、別のツールとの間で得られる値との差が小さくなるように決めましょう。完全に一致させることにあまり意味はないと思いますので、えいやっ、で決めましょう。決めたら仲間、未来の自分と共有。Consistencyがすべてです。
そもそもTm値はあくまでも何かを説明、決定して前に進むための指標で、予測値です（計算のためのパラメターは、先人たちが実験的にもとめたものが使われています）。わからないfactorについては見てみないふりをして計算します。そうしないと計算ができない。少し考えてみましょう。たとえばPCRのプライマーについて考えるとき、サイクル0での標的の濃度はさっぱりわからない、あるいは毎回違う、時には0、しかも目的ではない余計な核酸の量も予測不能で、、、という具合で、実際計算のための条件がきまらない。仮にきまったとしても、サイクルが進むにつれ、標的濃度はあがり（もとの標的濃度は変わりませんが増幅産物としての標的は増える）、逆にプライマー濃度は産物に取り込まれて下がっていく。え？ということは厳密にいえばTmの計算はできないし、ものすごくダイナミック！そう。でも計算するんです。今までだってそうしてきたし、それに基づいた私のプライマーは動く！そうです！それでいいんです。
ですから、思い切って、プライマーと標的が同じ濃度だとするとき、というようにパラメターをとります。Tm_NN()のdnac1、dnac2を使えばいろいろに試すことができますので、様々な条件を設定してコードを走らせてみましょう。エラーがでることはあるにせよ、コンピューターが壊れてしまう、というようなことは起こらないと思います。大切なのは、手を動かして試し、どんな時にTm値がどう影響されるか、想像するためのヒントを得ることです。
Tm計算その２：ミスマッチを含む場合の計算（とlimitation）
さて、次はミスマッチを含む場合の計算です。この場合、調べたい二つの配列をseq=とc_seq=で与えます。ここで注意が必要なのは、c_seq=はcomplementary sequenceであるということです。うっかりreverse complementを与えないように注意しましょう。seq=に与える配列は5' -> 3'、c_seq=に与えるのは3' -> 5'の向きということになります。
試しにこれまたテキトーに打った以下の配列を使ってみます。これからミスマッチを導入する場所を大文字にしてあります。もとの配列ではA:Tのペアです。
AをT、C、Gなどと順に変えてTm値がどうシフトするのかをみてみましょう。
以下のようなコードを書いて実行してみます。
（ちなみに、塩基文字のcaseは内部で大文字に変換されるので、大文字と小文字が混ざっていても問題ありません。Caseでのちょっとしたマーキングは便利ですから、ある場所の塩基を入れ替えるような作業をする時には、cAsEでのマーキングを使いましょう）
import Bio.SeqUtils.MeltingTemp as mt
# 1. only seq= is given
est_tm = mt.Tm_NN(seq='tactgcgatttAgcatcgatgaa',
Na=50, Mg=1.5, dNTPs=0.2,
dnac1=100, dnac2=100)
print('1. only seq arg:', est_tm)
# 2. give c_seq, no mismatch
est_tm = mt.Tm_NN(seq='tactgcgatttAgcatcgatgaa',
c_seq='atgacgctaaaTcgtagctactt',
Na=50, Mg=1.5, dNTPs=0.2,
dnac1=100, dnac2=100)
print('2. A:T match:', est_tm)
# 3. T:T mismatch
est_tm = mt.Tm_NN(seq='tactgcgatttTgcatcgatgaa',
c_seq='atgacgctaaaTcgtagctactt',
Na=50, Mg=1.5, dNTPs=0.2,
dnac1=100, dnac2=100)
print('3. T:T mismatch:', est_tm)
# 4. C:T mismatch
est_tm = mt.Tm_NN(seq='tactgcgatttCgcatcgatgaa',
c_seq='atgacgctaaaTcgtagctactt',
Na=50, Mg=1.5, dNTPs=0.2,
dnac1=100, dnac2=100)
print('4. C:T mismatch:', est_tm)
# 5. G:T mismatch
est_tm = mt.Tm_NN(seq='tactgcgatttGgcatcgatgaa',
c_seq='atgacgctaaaTcgtagctactt',
Na=50, Mg=1.5, dNTPs=0.2,
dnac1=100, dnac2=100)
print('5. G:T mismatch:', est_tm)
実行結果はこうなりました。
1. only seq arg: 60.754317191625375
2. A:T match: 60.754317191625375
3. T:T mismatch: 57.1127923250296
4. C:T mismatch: 55.96567268904096
5. G:T mismatch: 58.33014311575482
１と２の結果が同じになっています。１では相補配列を与えませんでしたが、内部でperfect matchのものが生成されて計算が行われるので、明示的に相補配列を与えた２の場合と結果が同じになるわけです。
さて、ミスマッチの場合の結果を見てください。ミスマッチの種類によりずいぶん値が異なります。とくにC:Tのミスマッチでのシフトが大きい。一般にCがからむとシフトが大きめ、Gがらみは影響が小さめになる傾向があります。どの程度のシフトが出るかは、ミスマッチの両側の塩基コンテキストにも影響されます。興味のある方は、ミスマッチを固定し、両側の塩基をすこしづつ変えてTm値がどう影響されるかを調べてみてください。
こうしていろいろと計算してみますと、ミスマッチひとつでしょ、と侮れない場合があることがわかってきますし、逆にミスマッチで結合の差を出したいときには、どのミスマッチをとらせるのが効率的か、を考えることができるようになります。
ここでひとつ、計算上の制限を書いておきます。それは、連続する位置でのミスマッチについては計算ができないということです。
例えば以下のような場合です。
これで計算しようとしますと、
import Bio.SeqUtils.MeltingTemp as mt
# 6. G:T, A:C consecutive mismatch
est_tm = mt.Tm_NN(seq='tactgcgatttGAcatcgatgaa',
c_seq='atgacgctaaaTCgtagctactt',
Na=50, Mg=1.5, dNTPs=0.2,
dnac1=100, dnac2=100)
print('6. G:T, A:C consecutive mismatch:', est_tm)
こんなエラーが出て計算してくれません。
ValueError: no thermodynamic data for neighbors 'GA/TC' available
GA/TCに対するデータがないよ、と言っています。そうなんです。連続するミスマッチが２塩基のフレームに入ったときの値が提供されていないので、計算ができないのです。（いまのところ）Biopythonの提供するTm_NN()では、連続したミスマッチは扱えないと覚えておきましょう。
連続していなければ該当するパラメターがコード内のデータテーブル（dictの形になっています）にあるはずなので、複数のミスマッチを含む、以下のようなものでも計算できます。
応用編：バッチ処理で一気にTm計算
プライマーやプローブの配列が書かれたfastaファイルがあるとします。
>primer01 self annealing
CGCGATTGGCCAATCGCG
>primer02 mostly A and T
TTTTTAGAATTATATATAATATA
>primer03 G4 strech included
CCATATATCGGGGTTATATAA
このファイルをとってそれぞれの配列のTm値を計算し、id lineの最後に書き足してくれるコードを書いてみます。出力はSTDOUTにすることにします。
from Bio import SeqIO
import Bio.SeqUtils.MeltingTemp as mt
import sys
# fixed parameters for Tm_NN()
# change here for your use
NA = 50 # in mM
MG = 1.5 # in mM
DNTPS = 0.2 # in mM
DNAC1 = 100 # in nM
DNAC2 = 100 # in nM
# get fasta file path
in_mfa = sys.argv[1]
# main part
for record in SeqIO.parse(in_mfa, 'fasta'):
est_tm = mt.Tm_NN(seq=record.seq,
Na=NA, Mg=MG, dNTPs=DNTPS,
dnac1=DNAC1, dnac2=DNAC2)
print('>' + record.description + ' Tm=' + '{:.2f}'.format(est_tm))
print(record.seq)
走らせてみますと、こんな感じにTm値がid lineの最後にはいります。
一歩踏み込む：Tm値算出に使われた、dH, dSを引っ張り出す
さて、最後にTm計算で使われたdH、dSを引っ張り出すべく、Tm_NN()のコードをプチハッキングしましょう。まず、この関数が書かれているファイルのpathを以下のようにして調べます。mt.__file__でMeltingTemp.pyのありかがわかります。
Pathがわかったら、MeltingTemp.pyをeditorで開き、myMeltingTemp.pyと（でも）してオリジナルのファイルと同じdirectoryに保存します。これを編集します。
Tm_NN()の定義までずーっとスクロールしてください。850行のあたりから始まります。このコードをじっくり見ますと、どのようにしてTm値が計算されるのかがわかりますので、興味のある方はたどってみてください。また、ファイルの前半に、すべてのthermodynamics parametersがあります。それぞれもとになる論文も明記されています。
最後にsalt correctionがありますが、その手前にある、以下の式でTmが計算されています。delta_hとdelta_sが含まれていますね。
melting_temp = (1000 * delta_h) / (delta_s + (R * (math.log(k)))) - 273.15
最終的にreturn melting_tempでTm値だけが返されるので、ここを以下のように変更してdelta_hとdelta_sも返るようにします。
return melting_temp, delta_h, delta_s
こうしておくと、dH (kcal/mol)、dS (cal/mol k)を得ることができ、
dG = dH - ((273.15 + 37)*dS/1000)
でdG°37 (kcal/mol)を計算することができます。
Note
プチハッキングしたmyMeltingTemp.pyを使うときは、
# import Bio.SeqUtils.MeltingTemp as mt # ではなく、
import Bio.SeqUtils.myMeltingTemp as mt
とするのをお忘れなく。
おわりに
さて、いかがでしたでしょうか。コードでTm値の計算ができるとなると、プライマーとして使えそうなTm=60の候補ストレッチを、与えられたゲノム配列から網羅的に作り出すとか、ミスマッチありでの計算ができることを使ってSNP検出のプローブをデザインするとか、いろいろと応用が考えられます。
あと、ソースコードをじっくりたどると、Tmの計算について大変勉強になりますし、コードを書いてくれた方々への尊敬の念もあらためて湧いてきます。ぜひ、おおー、などと感心しながら使ってみてください。