NMRによる生体高分子の動的秩序形成解析に向けたベイズ推定に基づく構造最適化計算

Abstract

核磁気共鳴分光法 (Nuclear Magnetic Resonance; NMR) は，溶液状態や生きた細胞内のようなより生理的条件に近い状態の生体分子の3次元立体構造やダイナミクスを原子分解能で解析できる現在唯一のツールである．一方で，NMR測定により得られる直接的なデータは，複数のピーク信号からなるスペクトルであることから，ここから分子の構造情報を正確に抽出し，分子の3次元構造を NMR立体構造計算によって再構成する必要がある．特に不安定分子や不均一系の測定データは，感度が低く，多くのノイズを含むため，構造情報の抽出・解析が困難となる．従来の立体構造計算法では，こうしたスパースデータから正確な構造決定を行うことは難しかった．そこで我々は，Riepingらによって提案されたベイズ推定を用いた新しい立体構造計算手法を応用し，この手法に複数の改良を加えた新規手法を開発，NMR立体構造計算ソフトウェアCYANAに実装した．新しく開発した手法では，NOEクロスピークの自動解析とAmberの物理ポテンシャルを用いた構造サンプリングにより，構造アンサンブルからなる事前確率分布を得ることができる．サンプリング手法には，ギブスサンプラーによるマルコフ連鎖モンテカルロ法 (MCMC)と分子動力学計算を組み合わせたハイブリッドモンテカルロを採用した．さらに，タンパク質の構造座標変数は膨大で関数空間が極めて広いため，レプリカ交換MCと組み合わせることで計算の効率化を図った．本手法の性能は，既知構造から作成したシミュレーションデータと，実データをランダムに間引いて情報量を減らしたデータ，それぞれを作成し，従来法と比較することで検証した．検証の結果，本手法は従来法と比較して，情報量を大きく減らしたスパースデータに対しても十分に正確な構造決定ができることを示した．

1 はじめに

タンパク質，DNA, RNAなどの生体高分子は，特定の，もしくは複数の3次元立体構造とダイナミクスを持ち，様々な生命現象に必須の機能を生みだしている．したがって，生体分子の生命機能や疾病機構の分子レベルでの理解，およびこれを基にした創薬開発においては，それぞれの生命過程，或いは疾患に関与するタンパク質の立体構造研究が不可欠である. 現在,生体分子の立体構造を原子分解能で解析可能な主な手法は，X線結晶構造解析，クライオ電子顕微鏡，および核磁気共鳴 (NMR) 分光法の3つが挙げられる．この中でNMRは，より生理的条件に近い環境での分子の立体構造やダイナミクスの情報を得られることに際立った特徴を持つ． X線結晶構造解析のように分子を結晶化する必要がないことから，NMRは生体分子の構造解析法の1手法として確立し，現在世界で12,000 以上の立体構造が決定され，データベース (Protein Data Bank; PDB) にその座標データが登録されている．また系にほとんどダメージを与えることなく計測可能であることから，生きた細胞中のタンパク質構造やダイナミクスを直接観測するin-cell NMR法なども開発されている [1,2]．ナノ秒から数時間にわたる広範囲の時間スケールでの生体分子のダイナミクスを観察することもNMRの際立った特徴として挙げられる．一方で，溶液NMR分光法の問題の1つは，他の測定法に比べ感度が低く，長時間の測定が必要な点にある．したがって，解析に十分な分解能と感度をもつNMRスペクトルの取得には，高濃度で長い時間安定な試料を要求する．このため，不安定試料から得られた質の低いNMRデータにおいても，正確な立体構造決定が可能な計算手法が強く望まれてきた．

NMRデータには生体分子内の原子の様々な空間情報が含まれるが，観測データそのものはスペクトル上のピーク信号の位置と強度にすぎないため，ここから構造情報を抽出し，より解釈しやすい形式である分子構造に可視化する必要がある．この解析過程を一般にNMR立体構造計算と呼び，観測により得られた構造情報を満たす分子構造のモデルを構築する．一般的な従来法では，タンパク質の場合，アミノ酸鎖が直鎖状に伸びた構造を初期モデルとし，実験データから得られる原子間距離および二面角を満たす最適構造を分子動力学計算 (MD)を用いたシミュレーテッドアニーリング法 (SA) により探索する．

これまでに，ノイズが多いもしくはデータ量の少ない実験データにおいても，より正確な構造決定が可能な様々なNMR立体構造計算法が提案されてきた [3,4,5]．その中で特に，2005年にNilgesらによって提案されたInferential Structure Determination (ISD) [6,7]と呼ばれるベイズ推定を応用した立体構造計算は，従来法に比べて計算精度が大きく向上し注目を集めた．ベイズ推定は，事前知識を確率分布で表し，事後確率分布を求めることで，原因事象を確率的に推論できる手法として広く用いられている．ISDは，従来法と比較して計算時間を有するものの，スパースデータの欠損部の推定としても有効であり，潜在変数を事前確率分布として定義し，マルコフ連鎖モンテカルロ (MCMC) でサンプリングすることにより，実験データからより確率的に構造情報を引き出すことができる．

本研究で用いるCYANAプログラム [8,9] は，現在世界の数百の企業，大学で使用され，生体系NMR領域の業界標準としてこの分野を代表するソフトウェアである．また，生体分子の分子モデリング等にも広く使用されている [10]．2面角系分子動力学計算 (torsion angle MD; TAMD) [11]を採用しているため，デカルト座標空間におけるMDよりも長いステップサイズを用いることができ，収束した最終構造をより高速に得ることができる特徴を持つ．ここでは，CYANAに実装したベイズ推論に基づくNMR構造最適化計算法 (CYANA BAYesian inference; CYBAY) の実装と結果について報告する．NMR信号のピーク位置と強度には，多数の分子構造に関する情報を含んでいるが，ここでは，NMR立体構造計算において最も広く用いられている核オーバーハウザー効果 (nuclear Overhauser effect; NOE) [12] データにのみベイズ推定を適用した．本手法は，NOE自動解析と粗い全体構造を迅速に得る計算の初期段階と，この構造を元にベイズ推定により構造最適化を行う2つの段階とに分けられる．

2 理論

NMR法において最も重要な計測データの1つであるNOEクロスピーク強度には，近距離 (おおよそ5Å以下) の2原子間距離情報を含むため，ここから分子内の原子の相対位置を推定することが可能である．しかしながら，得られる構造情報は，近距離のみであることと，様々な要因によってデータの欠損も起こるため，スパースな情報となっている． NMR立体構造計算では，こうしたスパースな局所構造情報を多数組み合わせることで，最終的に分子全体の構造を構築する．問題は不良設定であることから，分子のモデル化と事前知識 (正則化項) の活用が必要である．最も単純な形では，実験値とモデルからの推定値との二乗誤差 χ² と，正則化項として分子の物理ポテンシャルエネルギーEからなる目的関数Tを最小化する．   

$T\left(\theta \right)={\chi }^2\left(\theta \right)+wE\left(\theta \right)$(1)

ここでθは2面角空間内の原子座標を示し，wは各項のバランスを調整する重み因子である．

ベイズ推定における構造アンサンブル評価のための目的関数は，事後確率分布の形で置き換えられる．   

$P\left(\theta |D\right)\propto P\left(\theta \right)P\left(D|\theta \right)$(2)

ここでのDは観測値を表す．本稿では，最初の試みとして，NMR立体構造解析の中核となる計測データであるNOEについてのみベイズ推定の枠組みで計算を行い，ピークの位置情報から得られるタンパク質主鎖2面角情報や水素結合情報等については，計算に使用しなかった．したがって，ここでは，(2) 式のDはNOEデータを指している．NOEクロスピークのピーク強度I_klは，2つの原子核スピンkとlの間の双極子相互作用によって引き起こされる．NOEピーク積分強度から求まる二原子間距離は，以下のように定式化される．   

$I_{kl}=\frac{\gamma }{r_{kl}{}^6}$(3)

ここで， r_klは原子核スピンkとlの間の距離，γは補正定数である．NOEピーク強度比のダイナミックレンジは通常かなり大きいため，尤度関数Lは対数正規分布を用いて以下のように記述する．   

$P\left(I_{kl}|\theta \right)=\prod _{i=1}^n{L}_i\left(I_{kl}|\theta \right)$(4)

  

$L_i\left(I_{kl}|\theta \right)=\frac{1}{\sqrt{2\pi {\sigma }^2}{I}_{kl}}\exp \left(-\frac{1}{2{\sigma }^2}{\ln }^2\left(\frac{I_{kl}}{\mu }\right)\right)$

ここで，nはデータ数，σは標準偏差，μはピーク強度の平均値を表す．また，潜在変数および分子座標の事前確率分布は，以下のようにモデル化した．   

$P\left(\sigma ,\theta ,\gamma \right)=P\left(\sigma \right)P\left(\theta \right)P\left(\gamma \right)$(5)

  

$P\left(\theta \right)=\frac{1}{Z\left(\beta \right)}\exp \left(-\beta E\left(\theta \right)\right)$

  

$\begin{array}{c}\sigma ~G\left[a,b\right]\\ \gamma ~LN\left[{\mu }_{\gamma },{\sigma }_{\gamma }\right]\end{array}$

ここで，P(θ)は，現在の拘束条件を満たす分子立体構造のカノニカルアンサンブルを示し，分解関数Zと逆温度βにより表現される．aとbは，ガンマ分布関数Gの形状と尺度母数，LNは対数正規分布を示す．

観測データ表現のための潜在変数 σと γ，および二面角空間上の原子座標θのサンプリングは，それぞれMCMCとMDにより交互に計算され，全体をハイブリッドモンテカルロ (もしくは，ハミルトニアンモンテカルロ; HMC)として実行する．タンパク質は各原子が極めて密に詰まった内部コア構造を形成しており，コア領域では原子間にほとんど移動空間がない．こうした密空間の座標探索はMCよりMDの方が効率の良いことが知られている [13,14]．従って，原子座標θはMDでサンプリングを行い，潜在変数σ，μおよびγは，ギブスサンプラーにより探索する．一般的に，ベイズ統計におけるHMCは，説明変数が高次元の際のメトロポリス判定の採択率を向上させる目的で，仮想的な質量を定義して加速度を求めることにより効果的なステップ方向を決定する手法を指す [15]．本手法は，原子の質量が分子モデルから明確に定まっているため一般的なHMCとは一部相違する．MD計算の目的関数は以下のように定義される．   

$T\left(\theta \right)=\beta E\left(\theta \right)+L\left(\theta ,\gamma ,\sigma \right)$(6)

  

$L\left(\theta ,\gamma ,\sigma \right)=\frac{1}{2{\sigma }^2}\sum _{i=1}^n\ln \left(\frac{I_i{r}_i{}^6}{\gamma }\right)$

  

$\begin{array}{c}E\left(\theta \right)={\sum E}_{dihedral}+\sum E_{vdw}+{\sum E}_{vdw14}\\ +{\sum E}_{electro}+{\sum E}_{electro14}+{\sum E}_{GB}\end{array}$

ここでE_dihderal, E_vdw, E_vdw14, E_electro,E_electro14, および E_GB は，それぞれ，2面角項，Lennard-Jones項 (LJ)，1–4 LJ項，静電相互作用項，1–4静電相互作用項，generalized Born (GB) モデル [16] による溶媒和項を表す．TAMD計算では，原子間の結合長，結合角は固定されているため，これらのエネルギー項は上記式には含まれない．ポテンシャル関数および関数パラメーターは，Amber ff03 [17] に定義されたものを用いている．一方で，E_vdw項はTAMDにおける原子間結合長と結合角の固定により，エネルギー項が極端に大きくなることから，TAMDに最適化させるため，”1–5”,”1-6”,”1-7” ファンデルワールス相互作用のみ斥力成分を弱めたソフトポテンシャル項 [18] を用いた．   

$E_{\text{soft}}=\epsilon \left[C\left({\left(\frac{r-r^{*}}{r_0-r^{*}}\right)}^2-1\right)+\left(1-C\right)\left({\left(\frac{r-r^{*}}{r_0-r^{*}}\right)}^3-1\right)\right]$(7)

ここで，ε は各原子種に対応するLJポテンシャルの最小エネルギー値, r* はLJ関数が最小となる原子間距離，r₀ はLJ関数が0となる原子間距離，Cはスケーリング因子 (0 ≤ C ≤ 1) である．本計算では，スケーリング因子Cを”1-5”,”1-6”,”1-7” ファンデルワールス相互作用についてそれぞれ0.8，0.6，1.0に設定した．”1-5”,”1-6”,”1-7” 共有結合を介した2原子間距離rがr*より小さい際に，(7) 式がAmber ff03のLJ関数に置き換わって計算される．

また，ピークの重なり等により一意に原子ペアを確定できなかったピークについては，m個の候補原子ペアに対して，以下のように原子間距離が計算される [19]．   

$\overline{r}={\left(\sum _{j=1}^m{r}_{k_j{l}_j}{}^{-6}\right)}^{-1/6}$(8)

　生体分子のサイズが大きくなる，およびデータに曖昧さが増えることで自由度が増大する問題については，複数の異なる温度を持つ独立したMC計算 (レプリカ) を実行し，メトロポリス判定によりレプリカ間の温度を任意の間隔で交換するレプリカ交換MC [20] により空間探索の効率化を図った．

3 手法とデータ

Figure 1に本手法の主な流れを示す．本手法の第一段階では，CYANAプログラムによるNOE SpectroscopY (NOESY) クロスピークの自動解析と立体構造計算により，タンパク質のおおよその全体構造を得る．この段階では， NOEピーク積分強度の補正定数γは，スペクトル中のすべてのピーク強度の中央値を任意の距離 (通常4.0Å) に一意に決める．あるスペクトルから得られる全ての距離制限拘束は，この補正定数を用いてピーク強度から計算する． この条件をもとに，NOESYクロスピーク解析 (NOE帰属)と立体構造計算を再帰的に7回繰り返し，NOE帰属と立体構造を段階的に収束させる [21]．この計算過程は，初期値を変えた計算を独立に100回実行し，得られた100構造の中から最も目的関数値の低い上位20構造を選択する．次に，この 20個の構造に対して，Amber ff03ポテンシャルパラメーターを元に原子間の結合長と結合角を最適化させる．こうして結合長，結合角を Amber ポテンシャルに補正した2面角座標を元に，レプリカ交換HMCを実行しベイズ推定に基づく構造最適化を行う．

Figure 1.

 Structure refinement by Bayesian inference.

3.1 従来法によるタンパク質立体構造計算

新規手法との比較のため，従来法による立体構造計算も同時に行った．従来法は，CYANAに実装されている標準手順を実行した [8]．ランダムに生成した初期構造100個についてNMR実験データより得られた距離・角度拘束条件とともに，独立にTAMDを実行する． TAMDは 10000ステップ回行い，SA法により徐々に構造を収束させた．構造の最適化は，上記の基準と同様に100構造の中から 20構造を選択し，OPALpプログラム [22] により行った．OPALpでは，デカルト空間での共役勾配法によるエネルギー最小化を行う．ここでのポテンシャル関数は， Amber ff94 [23] を，水モデルにはTIP3P [24] を用いた．計算に用いた実験データは，比較のため，従来法と新規手法で同一のものを入力として用いた．

3.2 データセット

本手法の性能評価のために，2種類のデータを準備した．(1) 既知の構造から再構築された完全にシミュレートされたデータ， (2) 実測のスペクトルからピークをランダム削除することによって構造情報をスパースにしたデータ．(1) のデータセットで本手法の理論性能を評価し，(2) のデータセットでは，より現実的なデータに対する性能評価をする．完全にシミュレートしたNOEピークリストは，ヒト由来のユビキチンタンパク質 (76 aa, PDBID 1D3Z) の立体構造から作成した．ここでは，それぞれのタンパク質に対して最も一般的なデータセットである3D ¹H-¹³C NOESY-HSQC (¹³C-separated NOESY), 芳香環シグナル選択的3D ¹H-¹³C NOESY-HSQC (¹³C-separated aromatic NOESY)，および 3D ¹H-¹⁵N NOESY-HSQC (¹⁵N-separated NOESY) クロスピークを生成する． NOEは，通常2.4–5.0Åの範囲内に存在する水素原子間で観測されることから，既知の原子座標をもとに (3) 式の関係に従って，この範囲にある原子間相関ピークを作成した．個々のピークの積分強度は，実験値に一定幅のノイズが含まれることを仮定し，ランダムに，任意のγ値を平均としたガウス分布より作成した．また，従来法では解析の難しいシグナル・ノイズ (S/N) 比の低いデータを想定したデータセットも別途作成した．ここでは，ピーク強度の弱い信号はノイズに埋もれやすいと仮定して，ランダムかつ積分強度の弱いピークが優先して指数関数的に間引かれたデータを作成した．一般的に，強度の弱いNOEピークは，空間的に遠距離の原子間距離情報を含むため，強度の強いピーク (近距離情報) に比べ，情報量が大きく全体構造の決定における比重も大きい．したがって，強度の弱いNOEピークが消失しやすい傾向は，構造決定により価値の高いデータをより多く失うことを意味する． Table 1は，こうして得られたユビキチンのNOEクロスピークの数を示す．本研究では，NOEデータのみに新規手法を採用したため，計算には，通常のNMR立体構造計算で用いられる水素結合情報や2面角角度拘束は使用しなかった．またこの性能評価試験では，従来法と新規手法との純粋な構造計算の精度比較を行うため，NOEピークはすでに解析 (帰属) された段階から開始し，CYANAによるNOE自動解析機能は用いていない．

Table 1. Parameter sets for creating simulated NOESY peak lists of ubiquitin

		γ			σ
	peaksa	setb	conv.c	Bayesian	setb	conv.c	Bayesian
original data:
13C-NOESY	3217	22.87	22.58	22.47	2.29	-	2.34
13C-aro NOESYd	110	18.14	17.75	17.84	4.67	-	4.57
15N NOESY	823	16.46	16.23	16.22	3.75	-	3.70
filtered data:
13C-NOESY	124	20.49	19.87	20.47	2.23	-	1.79
13C-aro NOESYd	17	12.65	12.67	12.76	3.97	-	2.28
15N NOESY	89	13.31	14.04	12.98	3.41	-	2.24

^a number of used peaks

^b γ and σ values used for creating the fully simulated NOESY peak lists

^d conventional structure determination by CYANA-OPALp

^{c 13}C-separated aromatic NOESY

2番目の評価試験として用いた実データには，高度好熱菌由来のTTHA1718タンパク質 (TTHA1718; 66 aa)，シロイズナズナ由来のENTHA-VHS domain At3g16270 (ENTH; 140 aa), およびヒトfeline sarcoma oncogene Fes 中のSH2ドメイン (SH2; 114 aa) の3種類のNMRクロスピークリストを用いた．ピーク除去は，S/N比の低いスペクトルを想定し，上記と同様に強度の弱いピークから優先的に除かれるように行った．削除するピーク数は，従来法で得られる構造が，正解構造からのタンパク質主鎖 (H^N, N, C^α, C') の root mean square deviation (RMSD) 値 1.0–4.0Å程度ずれる数に設定した．本研究では， RMSD値を指標としたときに，従来法でこの程度構造が不正確となるデータセットの計算の改善を対象としている．タンパク質分子の運動性や，データの不完全性を考慮すると，正解構造 (参照構造) に対するタンパク質主鎖RMSD値が1.0 Å未満の構造を目指すことは意味が薄く，正解構造からのRMSD値が4.0 Åよりずれた構造を修正するのはタンパク質の原子座標の広大な関数空間を考慮すると容易ではないためである．ピークの位置情報 (化学シフト値) は，NMR実験データのデータベース (Biological Magnetic Resonance Data Bank; BMRB) に登録されている値を用いた．BMRB登録IDは，TTHA1718，ENTH，SH2それぞれ11035，5928，6331である．

4 結果と考察

4.1 シミュレーションデータによる検証

初めに，ユビキチンタンパク質の立体構造から，逆計算により作成したNOEクロスピークリストを用いて，本手法の性能を検証した．従来のCYANA-OPALp法と新規ベイズ法は，比較のために同一のデータセットを用いて独立して実行した．Table 1は，既知構造からの逆計算に使用した任意の補正定数γおよび標準偏差σ値，および従来法 (CYANA-OPALp) とベイズ推定を用いた新規手法 (CYBAY) ，それぞれで推定した値を示す．ピーク強度には，距離情報だけでなく，タンパク質の運動性や実験ノイズ等複数の要因も含まれていると仮定し，σのガウスノイズを加えて作成した． γとσは，MC計算の中で事後確率値の変動がほぼ一定となったステップ回数からサンプリングを開始した (データは省略) ．また，レプリカ交換MCでは，最高温度および最低温度を有するレプリカが1回以上交換されており，構造座標を含めた関数空間のサンプリングが十分であると判断するステップ回数の計算を実行した．γとσの得られた分布から推定される値および得られた分布をTable 1 と Figure 2cに示す．

Figure 2.

 Bayesian structure calculations with fully simulated data of ubiquitin. (a) Superpositions of the 20 reference structures from the full data sets (red) conformations yielded by the conventional (light green) and Bayesian method (blue), showing the backbone (N, C^α, C'; left) and all heavy (right) atoms. (b) Distributions of γ and σ for the ¹³C-separated NOESY of ubiquitin. Insets show trajectories of γ and σ in the HMC sampling. (c) Histograms and distribution curves of γ for the ¹³C-, ¹⁵N-, and ¹³C-separated aromatic NOESY spectra. The distributions of γ in the input peaks (black) and predicted by the Bayesian method (red), as well as the values determined by CYANA (green arrows) are shown.

既知構造から逆計算により得られた元のピーク数は，4000以上あり一般的なNMR立体構造計算には十分な数である．しかしながら，ここでのピーク強度にはσ値分のガウスノイズを加えてあり，ピーク強度から補正定数γと原子間距離情報を正確に推定する必要がある．

Table 1に示すように，従来法とベイズ法いずれもかなり正確にγ値推定している．従来法であるCYANA標準の推定法 ("手法とデータ”を参照) は，単純で高速な計算にも関わらず十分な正確さでγ値を推定している．しかしながら，この方法ではγの値を一意に決定するのみであるので疑似ノイズの大きさを示すσは推定できない．一方で，ベイズ法では分布の推定を行うため，データの分布を示すσも同時に推定されている．本提案法は，説明変数γの元の分布も極めて正確に推定できている (Table 1 & Figure 2c) ．これによりピーク強度からの距離の見積もりがかなり正確に推定できることになる．

得られた構造は，従来法では，ピーク強度から推定される距離拘束は十分であるはずにも関わらず，ピーク強度の含まれるノイズを加味できていないため，疑似ピークを作った際に用いた参照構造 (赤) からわずかに構造 (緑)のずれが観察される (Figure 2a). 参照構造 (正解構造) に対する主鎖原子 (H^N, N, C^α, C') のみのRMSD値 (1–70残基領域)は1.94 Åであり，側鎖原子を含むすべての重原子では2.12 Åであった (Table 2)．一方，新規手法では，ノイズの大きさσ値も計算に加味されているため，タンパク質側鎖まで参照構造とほぼ同一の構造 (青) が得られている (Figure 2a)．ユビキチンのC末端領域 (71–74残基) のみ参照構造からのわずかなずれが見られるが，これは，この領域が構造の中心から突き出たような配置をとっており遠距離の拘束条件がないためである．本稿で対象としているNOEデータでは，5Å程度の近距離の原子間距離情報のみしか得られないため，他の領域から孤立してるC末端領域には周辺原子との距離情報が得られない．通常のNMR実験では，このような領域に対しては別の種のNMRデータで情報を補完するが，本計算ではモデルの単純化のためにこれらは用いていない．よって，この領域(71–74残基) の構造の違いは本研究では，評価の対象から外している．また，タンパク質には溶液中で一定のダイナミクスを持っており，このサイズのタンパク質分子では1Å程度の揺らぎがあることを加味すると，新規手法により得られた参照構造に対するRMSD値が主鎖原子のみで1.01Å，側鎖原子を含むすべての重原子では1.31 Åの値となる構造は，十分に正確な構造を推定できている．

Table 2. RMSD values (1-70) of ubiquitin structure calculation with simulated data

	original data		filtered data
	conv.a	Bayes	conv.	Bayes
RMSD to mean (Å)a	0.33/0.56	0.24/0.42	1.52/1.99	1.08/1.40
RMSD to reference (Å)b	1.94/2.12	1.01/1.31	2.92/3.28	2.15/2.53

^a RMSD to the mean structure of the structure ensemble (backbone/all heavy atoms)

^b RMSD to the reference structure (backbone/all heavy atoms)

^c RMSD to the mean structure of the structure ensemble (backbone/all heavy atoms)

次に，S/N比の悪いスペクトルを仮定し，強度の弱いピークを優先してランダムに間引いたデータセットを用いた計算を行った．このデータは，ピーク数の合計が230個となっており，得られる距離拘束も大幅に減っている (Table 1) ．これにより従来法で得られた構造アンサンブルは，主鎖RMSD値が1.52 Åと構造の収束がかなり悪くなっている．一方で，提案法では主鎖RMSDが1.08Åであり，構造の収束性が改善している．また，参照構造に対する主鎖RMSD値も，2.92 Åから2.15 Åへと0.8 Å程度の構造の改善が見られた．

これらの結果から，ベイズ推定を用いた立体構造計算は，NOEクロスピーク強度に一定のノイズを含む場合やピーク数が減少することによる距離拘束条件の減少したスパースデータに対しても，従来法より高精度の構造決定が可能であることを示してる．

4.2 不完全な実測データを用いた検証

既知構造からシミュレーションした理想データに対しては，本提案手法の従来法に対する優位性を示すことができたことから，次に実測データを入力として性能評価を行った．タンパク質ENTH，TTHA1718およびSH2の実験データをもとに， NOEクロスピークの約55∼90％をランダムに欠失させた3つのデータセットをそれぞれのタンパク質について用意した．Table 3は，計算に用いたデータセットと構造計算により得られた最終構造のRMSD値を示す．いずれの計算でも，提案手法は従来法に比べ最終構造に対するRMSD値が小さくなっており，ベイズ推定を用いた本手法は，実測データに対してもより正確な構造を得ることができることを示している．

Table 3. RMSD values of three proteins

run	1		2		3
	conv.a	Bayes	conv.	Bayes	conv.	Bayes
ENTH: (1–130)
Peaksb	1229/6823		1010/6823		903/6823
RMSD to mean (Å)c	1.48/2.01	0.69/1.06	1.92/2.46	0.74/1.18	2.13/2.80	1.24/1.70
RMSD to reference (Å)d	1.92/2.21	1.25/1.86	2.41/2.74	1.04/1.65	3.51/3.77	1.86/2.47
SH2: (8–110)
Peaks	2832/5422		2175/5422		2033/5422
RMSD to mean (Å)	1.71/2.28	0.65/0.96	2.19/2.76	1.13/1.57	2.92/3.45	1.38/1.80
RMSD to reference (Å)	1.39/1.90	1.27/1.82	2.24/2.76	2.28/2.83	3.10/3.55	2.25/2.65
TTHA1718:
Peaks	472/3205		348/3205		299/3205
RMSD to mean (Å)	1.47/2.03	0.84/1.32	2.13/2.83	0.82/1.33	2.43/3.18	1.52/2.06
RMSD to reference (Å)	1.57/1.89	1.55/1.97	2.40/2.77	1.32/2.00	3.21/3.59	2.36/2.85

^a conventional structure determination by CYANA-OPALp

^b number of used/original peaks

^c RMSD to the mean structure of the structure ensemble (backbone/heavy atoms)

^d RMSD to the reference structure (backbone/all heavy atoms)

Figure 3では，それぞれのタンパク質の3回の構造計算のうちの1つの最終構造を，従来法と提案法について示している．ここで表示しているリボンモデル構造は，それぞれの計算の目的関数値もしくは事後確率の値が最も良い代表構造である．ENTHタンパク質の計算では，6823個の実験ピークのうち85%程度を削減した入力データを用いた結果について示している．データ数がこれだけ削られた場合，従来法 (緑)では参照構造 (赤) に対する主鎖原子のRMSD値 (8–130残基領域) が3.51Åとなっており，4本のαヘリックスの相対位置やループ領域の構造が大きくずれている．一方で，ベイズ推定を用いた本手法 (青) は明確にこれらの構造のずれが修正されている．また，Figure 3の右図には，サンプリングで得られたタンパク質主鎖C^α原子の座標分布から標準偏差1σ内の位置をチューブモデルで表している．従来法は，100回の独立したSAを実行し目的関数値の低い上位20構造を選択している．本提案法はMCでサンプリングした数千構造の分布を表示しているため，2つの手法を単純に比較できないものの，新規手法はC^α原子の座標分布の広がりが明確に小さくなっており，スパースデータを用いた際でも最終構造の曖昧さが小さくなったと解釈できる．

Figure 3.

 Structure refinement by Bayesian inference. Structures bundles of three proteins calculated from reduced data sets by the conventional CYANA method (green) and Bayesian refinement (blue) superimposed on conventionally determined reference structures from the full data set (red), ENTH (a), SH2 (b), TTHA1718 (c). Deviations of C^α atoms in the conformers yielded by the conventional and Bayesian method are shown as tube models.

SH2タンパク質の計算では，特にサンプリングされたアンサンブル構造内のRMSD値 (7–110残基領域)が3回の計算すべてで小さくなっており，構造がよく収束していることを示す． Run 3の計算では，3205個の元のクロスピークのうち60%程度データを削除したデータセットで計算を行っている．ピーク数をかなり削減しているため，従来法の計算では，βストランドやαヘリックスの相対位置が大きくずれており (Figure 3b)，参照構造に対する主鎖RMSDが3.10Åであったが，提案法により得られた構造では構造の明確な改善が見られる．

TTHA1718タンパク質の計算では，3回の独立した計算のうちの2回は明確な構造の修正が見られ，参照構造に対する主鎖原子のRMSD値は，2.40Åから1.32 Åと3.21 Åから2.36 Åとより正確な構造が得られた．3回の計算のうちの1つは，従来法で参照構造に対する主鎖RMSDが，1.57 Åから1.55 Åとわずかに値が小さくなっているが，これは計算誤差の範囲内であり，この計算についてはほぼ従来法と同一の結果となった．TTHA1718のRun 2の計算は，3205個の元のクロスピークのうち90%程度データを削除し348個の入力データで計算を行ったものであり，この構造についてFigure 3cに示した．参照構造に対する主鎖RMSDが1 Å程度小さくなっているため，リボンモデル構造で表示した際も，構造の明確な改善が見られる．さらに，右図のチューブモデル構造においても最終構造のアンサンブルが，従来法で得られた構造と比較してかなり収束していることが分かる．

5 結論

本稿では，CYANAプログラムにベイズ推定を用いた新規立体構造計算手法を実装し，NOEクロスピーク強度にノイズを含む，もしくはピーク数の少ないスパースデータにおいても従来のCYANA/OPALp法と比べ，得られる最終構造が大きく改善可能であることを示した．従来法では，NOEピーク強度から原子間距離への変換補正項や目的関数内の重み係数をユーザーが任意に定義する，あるいは単純な計算法で一意に決定していた．一方で，これら補正項や重み係数の決定は，タンパク質の取り得る構造アンサンブルに強く関連しており，逆に計算で得られる構造サンブルはこれらのパラメーターによって導出される距離・角度拘束条件によって制限される．したがって，実験データを元に，構造および複数の説明変数を同時に探索することは理にかなっている．NOESYスペクトルの自動解析は，通常，構造計算とスペクトル解析の繰り返し計算により達成され，計算コストが大きい．よって，本手法では，計算の初期段階では大まかな全体構造を得るために，単純化した目的関数を用いて従来の構造計算を行い，その後ここで得られた構造をベイズ推定を用いて最適化する計算手順を採用した．

本研究では，新規手法の有効性を証明した点に加え，生体系NMR分野の業界標準であるCYANAプログラムに本法の実装が成功した点も，使用者の観点からの意義が大きい．多くのNMR分光学者にとって，パラメーターやデータフォーマットが異なる複数のソフトウェアを用いて構造解析を行うことは，個々のパラメーターの違いの補正や，ファイルフォーマットの変換等，手間のかかる作業である． CYANAプログラムはNMRスペクトルの自動解析から構造最適化計算まで1つのソフトウェアで全自動処理可能であり，ユーザーに簡便なデータ解析を提供する．現バージョンのCYANAは， NMRスペクトルを得た後の，ピーク認識，原子核共鳴シグナルの解析，NOESYクロスピーク解析，立体構造計算，および構造最適化計算を含む一連のステップをすべて処理可能である．一連の処理の自動化は，全体処理を簡便化するのみならず，これまでヒトによって独立に解析されていた各ステップを，前後のステップの結果をもとに，プログラムが再評価可能であり，解析の信頼度を向上させることができる．例えば，構造最適化計算により得られた立体構造の情報をもとに，再度ピーク認識段階の確度の評価が可能となり，これを再帰的に繰り返すことでそれぞれの段階の質を向上できる．

本稿で開発したベイズ推定を用いたNMR立体構造最適化計算の利点の1つは，レプリカ交換およびギブスサンプラー法によって，他の説明パラメーターを用いてより広い関数空間を探索し，最終構造を確率分布の形でアンサンブル構造として得られることである．さらにこの手法を，固体NMRのデータ，生きた細胞内のタンパク質解析 (in-cell NMR) などのS/N比の低い様々なNMRデータに適用することでさらなる応用も期待できる．今後は，こうした計測データに適用し，性能評価することで，本手法の一層の改良を進めていく．

Acknowledgment

本研究は，科学技術振興機構 (JST) の戦略的推進事業 (CREST)，日本学術振興会 (JSPS) の科研費基盤C (nos. JP25440032 and JP15K06979)，新学術領域 (nos. JP24113720, JP26102538, JP25120003, JP16H00779, JP15H01645 and JP16H00847) の研究費支援を受けて行われた．

参考文献

• [1]   D. Sakakibara, A. Sasaki, T. Ikeya, J. Hamatsu, T. Hanashima, M. Mishima, M. Yoshimasu, N. Hayashi, T. Mikawa, M. Wälchli, B. O. Smith, M. Shirakawa, P. Güntert, Y. Ito, Nature, 458, 102 (2009). doi:10.1038/nature0781419262674
• [2]   T. Ikeya, A. Sasaki, D. Sakakibara, Y. Shigemitsu, J. Hamatsu, T. Hanashima, M. Mishima, M. Yoshimasu, N. Hayashi, T. Mikawa, D. Nietlispach, M. Wälchli, B. O. Smith, M. Shirakawa, P. Güntert, Y. Ito, Nat. Protoc., 5, 1051 (2010). doi:10.1038/nprot.2010.69
• [3]   C. D. Schwieters, J. J. Kuszewski, N. Tjandra, G. Marius Clore, J. Magn. Reson., 160, 65 (2003). doi:10.1016/S1090-7807(02)00014-912565051
• [4]   W. Rieping, M. Habeck, B. Bardiaux, A. Bernard, T. E. Malliavin, M. Nilges, Bioinformatics, 23, 381 (2007). doi:10.1093/bioinformatics/btl58917121777
• [5]   T. Ikeya, D. Ban, D. Lee, Y. Ito, K. Kato, C. Griesinger, Biochim. Biophys. Acta, Gen. Subj., 1862, 287 (2018). doi:10.1016/j.bbagen.2017.08.02028847507
• [6]   W. Rieping, M. Habeck, M. Nilges, Science, 309, 303 (2005). doi:10.1126/science.111042816002620
• [7]   M. Habeck, M. Nilges, W. Rieping, Phys. Rev. E, 72, 031912 (2005). doi:10.1103/PhysRevE.72.031912
• [8]   P. Güntert, Prog. Nucl. Magn. Reson. Spectrosc., 43, 105 (2003). doi:10.1016/S0079-6565(03)00021-9
• [9]   P. Güntert, C. Mumenthaler, K. Wüthrich, J. Mol. Biol., 273, 283 (1997). doi:10.1006/jmbi.1997.12849367762
• [10]   K. Patel, K. R. Srivastava, S. Durani, Bioorg. Med. Chem., 18, 8270 (2010). doi:10.1016/j.bmc.2010.10.00321035349
• [11]   A. Jain, N. Vaidehi, G. Rodriguez, J. Comput. Phys., 106, 258 (1993). doi:10.1016/S0021-9991(83)71106-X
• [12]   K. Wüthrich, NMR of Proteins and Nucleic Acids, Wiley (1991).
• [13]   H. Yamashita, S. Endo, H. Wako, A. Kidera, Chem. Phys. Lett., 342, 382 (2001). doi:10.1016/S0009-2614(01)00613-3
• [14]   S. H. Northrup, J. A. McCammon, Biopolymers, 19, 1001 (1980). doi:10.1002/bip.1980.3601905067378544
• [15]   M. Girolami, B. Calderhead, J. R. Stat. Soc. Series B Stat. Methodol., 73, 123 (2011), (Statistical Methodology). doi:10.1111/j.1467-9868.2010.00765.x
• [16]   N. A. Baker, Curr. Opin. Struct. Biol., 15, 137 (2005). doi:10.1016/j.sbi.2005.02.00115837170
• [17]   Y. Duan, C. Wu, S. Chowdhury, M. C. Lee, G. Xiong, W. Zhang, R. Yang, P. Cieplak, R. Luo, T. Lee, J. Caldwell, J. Wang, P. Kollman, J. Comput. Chem., 24, 1999 (2003). doi:10.1002/jcc.1034914531054
• [18]   V. Katritch, M. Totrov, R. Abagyan, J. Comput. Chem., 24, 254 (2003). doi:10.1002/jcc.1009112497604
• [19]   M. Nilges, J. Mol. Biol., 245, 645 (1995). doi:10.1006/jmbi.1994.00537844833
• [20]   K. Hukushima, K. Nemoto, J. Phys. Soc. Jpn., 65, 1604 (1996). doi:10.1143/JPSJ.65.1604
• [21]   T. Herrmann, P. Güntert, K. Wüthrich, J. Mol. Biol., 319, 209 (2002). doi:10.1016/S0022-2836(02)00241-312051947
• [22]   R. Koradi, M. Billeter, P. Güntert, Comput. Phys. Commun., 124, 139 (2000). doi:10.1016/S0010-4655(99)00436-1
• [23]   W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell, P. A. Kollman, J. Am. Chem. Soc., 117, 5179 (1995). doi:10.1021/ja00124a002
• [24]   W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, M. L. Klein, J. Chem. Phys., 79, 926 (1983). doi:10.1063/1.445869