共分散法 covariance

NMR スペクトルにおいて、共分散 covariance 法はこれまでのフーリエ変換 Fourier-transformation 法に代わるかもしれない プロセス法の一つとして注目を集めていました。最初は 2D 1H/1H NOESY, TOCSY などのように直接測定軸と間接測定軸が同じ種類のスペクトルにおいて、後者の分解能を前者にまで高める方法として提唱されました。その後、２つの異なるスペクトル（例えば (I, K) 相関スペクトルと (L, K) 相関スペクトル）の間で covariance をとることにより、(I, L) スペクトルを生成する方法として使われました。これは更に例えば主鎖帰属の HNCA と HNCOCA スペクトルなどにおいて 13Ca の化学シフトをもとに HN(i) と HN(i-1) を相関させるといった方法にも適用されました。

Harden, B.J., Frueh, D.P. (2018) Covariance NMR processing and analysis for protein assignment. Methods Mol. Biol. 1688, 353-373. doi: 10.1007/978-1-4939-7386-6_16.

Covariance の計算方法は実はたいへん簡単でして、上の例では (I, K) 行列と (K, L) 行列の内積をとるだけです。ここで後者の行列は転置しています。この共分散はフーリエ変換前の時間軸データでも、フーリエ変換後の周波数軸データでもどちらに適応しても構いません（細かな違いはありますが）。苦労する点といえば、４次元などのギガバイト大容量のデータに適用すると小さな PC が気の毒になる点ぐらいでしょうか？したがってプログラミングの FOR 文を使いまくるのではなく、できるだけ内積の高速演算ライブラリーなどを使ってコーディングした方が良いでしょう。

ただし、共分散法には false-positive な偽ピークが出てしまうという大きな欠点がありました。例えば、K 次元に沿ってちょっとだけずれた２つのピークがあったとします。目で見るとずれていることがすぐに分かるようなレベルでです。例えば、HNCA と HNCOCA とで 13Ca のピークが半値幅ぐらいずれていたとしましょう。この場合に、この２つのアミド基の帰属を相関させるような間違いは目視ではまず起こり得ません。ところが、共分散をとると、きっちりと相関を示す偽ピークとなってしまうのです。これでは使い物にならず、失望して止めてしまったものです。

ところが、最近はこれを克服する方法も出てきました。過去の論文にちらっと書いてあったことなのですが、共分散をとる前に微分をとっておくのです（もしかすると、分散波形でも良いのかもしれません）。すると、２つのピークが K 次元に沿ってぴったりと揃っている時に共分散が正の値を示し、逆に少しずれている場合には、共分散は小さな値、時には負の値を示します。これにより、false-positive な偽ピークを見分けることができます。

さらにエラーを減らす方法が提案されています。HNCA と HNCOCA とで K 次元に現れる 13Ca のピークをもとに HN(i) と HN(i-1) を相関させていったとします。しかし、この 13Ca だけだと間違いが多いでしょう。そこで普通は 13Cb, 13Co なども隣り合うアミド基を相関させるための共通の「糊」として使います。同じように (I', K') 相関スペクトルと (L', K') 相関スペクトルの間で covariance をとって (I', L') スペクトルを生成したとします。そして (I, L) も (I', L') も正の値のみを残しておいてから、両者を要素ごとに掛け合わせます。すると、いずれか一方に偽ピークがあった場合でも、掛け算によりそれは消えてしまいます。 (I, L) と (I', L') の両方に相関ピークがある場合にのみ、掛け算のスペクトルにピークが残るという仕組みです。ちょうど 13Ca, 13Cb, 13Co 全てを通して２つのアミド基どうしが相関を示した時に連鎖帰属を確定するのと同じです。二次元どうしですと、あまりメリットが感じられないかもしれませんが、３次元どうし（(I, J, K) と (L, M, K)）で K 次元に沿って共分散をとり４次元 (I, J, L, M) とすると、もしかすると有用かもしれません。もちろん今まで通りに処理すると多くの偽ピークが出てしまいます。しかし、covariance の前に K 次元に沿って微分をとり、covariance 後にスペクトルどうしを掛け合わせれば、上手く行くかもしれません。

Covariance や掛け算処理をする場合、各次元のデジタル分解能を合わしておくことは重要です。そのためには、同じマシンで同じスペクトル幅で一連のスペクトルを測っておきましょう。ポイント数は違っていてもスペクトル幅が同じであれば、0-fill 後のデータ数を同じ値にすることでデジタル分解能を調整することができます。異なるマシンで測った２つのスペクトルでは、いろいろなミスマッチが起こると考えられます。まず、マシンの絶対的な温度が異なります。さらにピークの位置も少しずれます。DSS のピークで 1H, 13C, 15N 次元を全て調整することでかなり揃えることができますが、それでも限界があります。本当は、HNCA と HNCOCA のペアなども interleaved-manner で測る方がよいのでしょう。それぞれを２日間ずつ連続して測定したとしても、それなりのずれは（ロックの不安定性などから）生じ得ます。

しかし、それでも false-positive な偽ピークが出てしまうそうです。これが起こるのはかなり線幅が異なるピークどうしで共分散をとった時です。例えば、HNCO では感度が高く大きな 13CO のピークが観えている一方、HN(CA)CO では感度が低く小さな 13CO のピークしか観えていないような時です。すると、共分散の前に微分をとったり、後に他の４次元と掛け算をしたとしても偽ピークが出勝ちです。そのような場合も想定して、一応は共分散をとる前のオリジナルスペクトルも少しチェックした方が良いとのこと（それでは covariance を活用する意味がないのですが）。また、４次元 (I, J, K, L) では (I, J) から (K, L) への相関と、その逆の (K, L) から (I, J) への相関の両方が共存しているかどうかを確かめることが重要とのことです。

ここで covariance の意味をあえて言うと、これまでの方法ではピークを拾い忘れていたり拾い方が悪ければ、それで終わりでした。拾った後の「化学シフトの値だけ」をもとに連鎖帰属をしていきますので、どのようにピークを拾うかが、その後の連鎖帰属の効率を決めているようなところがありました。一方 covariance 法では、false-positive が出る程に、そのような「拾い忘れ」が無いことが利点といえます。

いったい covariance という処理によって、スペクトルから何が失なわれてしまうのでしょうか？個人的にはピークトップの情報ではないかと考えています。人の目でピークを判断する時は、等高線で描かれた楕円の画像全体を観ながら、その中心をピークトップとして認識します。必ずしも最高強度の地点とは限りません。楕円が歪んでいたりすると、実はもうひとつ別のピークとオーバーラップしているのではないかと脳は推測したりもします。すると、その重なり具合に応じてピークトップの位置を少しずらしたりもします。初心者の場合、ここが欠点となってきます。covariance をとる前に微分するという処理は、ちょうどこのピークトップをできるだけ目立たせようとしていることに相当するのではないでしょうか？微分値ではちょうどピークトップの箇所で大きく正負が入れ替わります。そのため、どれぐらいの幅で微分（差分）をとるかも重要な要素になってくるだろうと考えられます。このような点が議論されているかとも思ったのですが、下記のオリジナル論文も含め見つかりませんでした。

Bradley J. Harden, Scott R. Nichols, and Dominique P. Frueh (2014) Facilitated assignment of large protein NMR signals with covariance sequential spectra using spectral derivatives. J. Am. Chem. Soc. 136 (38), 13106–13109. DOI: 10.1021/ja5058407.

生きた博物館

過去に何かの記事で NMRBox を紹介しましたが、また改めて。

これは NIH からのサポートを受けている一種のクラウドシステムです。サーバにはいろいろな NMR に関するソフトがインストールされています。たいへん興味深いことは、もう古くなって使えなくなってしまったようなソフトもそこでは生（活）きていることです。もちろん今でも優秀なソフトはどんどん開発されてはいますが、どうしても新しいトピックスをテーマとしたソフトとなってしまっており、いわゆる古典的な機能をもったソフトが、気が付いたら消滅してしまっているということがよくあります。しかし、この NMRBox ではそのようなソフトが使えるのです。まるで博物館です。

また、non-uniform sampling で測定した４次元データなどは、プロセスするのも大変です。そもそも私のところも Windows-XP の古い PC などを Cent-OS linux などに使い回しているものですから、四次元 NUS をプロセスするとなると今でも大変です。さらに PC も古くなってくると、突如として壊れます。もうこの半年で３台が去ってしまいました。なんとかデータだけは Puppy を使って幸運にも救出できてはいるのですが。一方、NMRBox には NMRPipe の IST や Smile がありますし、メモリーも H/D も豊富に備えていますので、そのような資源について心配する必要はありません。つい先日まで四次元 NUS の Smile によるプロセスが自分の PC ではメモリー不足でうまく行かなかったのですが、やっと NMRBox でうまくいきました。

また、夜帰宅してからでも週末でもプロセスの続きを見たり、ピークを帰属することができます。実際にプロセスしているコンピュータはきっと NIH にあるわけで、そのディスプレーを Real VNC viewer を使って見ていることになります。家でも職場でも出張中でも講演中でも NMRBox にアクセスできるということは素晴らしいことです。

ちょっと変わった使い方ですが、環境のセットアップをうまく出来ない人がいた時に、こちらで NMRBox にログインして環境を整えてあげることができます。例えば、連鎖帰属に必要な全ての３次元４次元スペクトルをまず関連付けます。二次元 HSQC のピークをクリックすると、自動的に３次元４次元の対応するピークに飛んでいくように設定しておきます。その後でその人がログインすれば、もうすでに帰属のための環境は整っていますので、後は実際の帰属作業に集中してもらうことができます。これまでは、わざわざその人の PC のところに行って環境を整えていましたので、それが出来る時間が何かと限られてしまっていました。４次元を何個もセットアップしようとすると何日も通わないといけませんでした。しかし、今は自宅からセットアップできます。また、担当者が代わっても、次の人は環境をそのまま引き継ぐことができます。面白いことに２箇所から同時にログインすると、他者の作業の状況をリアルタイムで見ることができます。まあ、同時に同じピークをピックしないように気を付けないといけませんが。。。。ちょうど OneDrive や GoogleDrive を他人と共有しながら、一つの Word ファイルを同時に編集しているような感じになります。

惜しい点は、どうしても画面を日米間で転送するので、ちょうど国際電話をかけている程度の時差が生じてしまう点です。１秒以内なのですが。vi でテキストファイルを編集していると、どうもカーソルが１秒ほど遅れて着いて来るため、希望する行を行き過ぎることが多いです。まあこの辺りは仕方がないですね。

できれば、Cyana, MagRO, Topspin も入っていれば言うこと無しなのですが。Topspin はどうしても立ち上がらず、管理者に問い合わせたら、ちょっとしたミスマッチでもう少し待って欲しいとのことでした。期待したいと思います。Bruker さん、どうかよろしくお願いいたします。

１時間程前にスタートさせた四次元 NUS のプロセスが順調に進んでいるようです。それでは明朝の結果に期待して、そろそろお休みしたいと思います。もちろん私の PC の電源は切ってもプロセスは続きます。zzzZZZZ

今だに位相補正で苦しむ

何でもない NMR 実験のはずが、そのプロセスで意外にも時間を費やしてしまいましたので、その覚え書きとして残しておきます。その測定とは HCCH-COSY です。ちょっと高分子量の蛋白質になってくると、CCONH や HCCONH での側鎖のピークが急に観えなくなります。そこで、HCCH-COSY などアミド水素にまで磁化を移動させない方法で側鎖の帰属を試みました。

三次元の HCCH-COSY と言っても、FID としての直接測定軸（x）を除いた y, z 軸がどれに当たるのかが問題です。今回は HCcH-COSY と hCCH-COSY をとってみました。小文字の箇所が検出をスキップする核種を示しています。４次元に拡張すると間接測定次元をどの 1H, 13C に当てるかという問題はなくなるのですが、やはりそれなりに感度が必要になってきます（次元が一つ増えるごとに感度がルート２に反比例して落ちる）。

ところが、Br 社のパルスプログラムでは、どの HCCH-COSY が上記の各々に当たるのかをすぐには見つけられません。パルス図を見て初めて、そうか HCcH-COSY が hcchcogp3d で hCCH-COSY が hcchcogp3d2 だとやっと分かります。プログラムの後ろの「２」は、新旧のバージョン番号を示しているのだと思っていましたが、そうではありませんでした。

日本では地震が多いため、最近は皆？ NUS で測定します（たとえ 100% sampling であったとしても）。地震で変になった箇所のデータを後から除けるためです。それに私は窒素やヘリウムを入れるスケジュールをすぐに忘れて測定をスタートしてしまうので、三次元測定の最中にやむなく測定を止めないといけない事もしばしばです。そのような時に NUS で測定していると、一応は途中までのデータを捨てずに有効活用することができます（もちろん感度は落ちますが）。また、重水素デカップリングをかけた実験の場合には、NUS にしないと autoshim との干渉が起きてしまい、翌朝にはロック信号が底辺をごそごそと蠢いているという事態を見ることになります。

さて測定が終わり、まずは hCCH-COSY (hcchcogp3d2) をプロセスすることにしました。ここで困ったことはパルスプログラムには aqseq 312 と書かれていることです。これは F3, F1, F2 の順にサンプリングすることを意味します。In0=Inf1/2 という表現から F1 には D0 インクリメントが対応しており、in10=inf2/2 という記述から F2 には D10 インクリメントが対応していることが分かります。FID での検出（F3）を Hi とすると、D0 が 13Cj に D10 が 13Ci に対応することもパルスプログラムで調べておきます。さて aqseq 312 ですので、NMRPipe の fid.com では y と z のカラムがひっくり返るはずなのです（と思い込んでいました）。しかし、実は NUS の場合にはきっとこれが起こらないのです。というのは、NUS では y と z のどちらを先にインクリメントするのかという概念がそもそも無くなるためです。それに nuslist でも y と z のインデックス番号のカラムをひっくり返してはいませんでした。しかし、それに気付かずに fid.com での y と z カラムをわざわざ交換してしまいました。普通でしたら間違いにすぐ気付くはずなのですが、この時は y, z ともに 13C 軸なので違いがすぐには分からないのです。

こうして、Hi を横軸に Cj を縦軸に二次元として表示すると、スライス軸は Ci になります。するとピークが Cj 次元軸に平行に点々と縦に並ぶはずです。実際には感度が悪くて並ばなかったので、ますます気付けなかったのですが。中には何故 Hi と Ci を二次元表示しないのかと不思議に思われる方もおられるかもしれません。感覚的にはその方が分かり易いかもしれません。４次元ではどの軸を後ろに持っていくかはさらに興味深い問題です。理屈よりも実際に試してみてその時の良し悪しを実感してみた方が面白いでしょう。将来は目の前の空間に球が浮かんだような感じの表示（ホログラム）になるので、今のような議論は消えるでしょう。ホログラムで四次元はどう表示するのかは問題ですが。

次に HCcH-COSY (hcchcogp3d) のプロセスです。これは更に悲惨でした。F1 が 1Hj に F2 が 13Cj と異なる核種に対応しているので、両軸をひっくり返してしまうという間違いはないのですが、F1 軸は TPPI-States に加えてちょっと特殊な位相回しが入っていました。パルスプログラムの下の方に calph(ph3, -90) という記述があります。これは D0 をインクリメントした時に ph3 を 90 度逆に回すことを意味します。せっかく TPPI-State で x から y に +90 度進めたのに何故また逆転させるのか？これは一種の TPPI 法に似た概念を利用していまして、TPPI-States を基本としながらも、さらに t1 インクリメントごとに観測座標を 1/4 回転だけ逆に回すのです。すると、磁化ベクトルがスペクトル幅の 1/4 だけ遅く回っているように見えるのです。これを FT する際にはスペクトルの中心をスペクトル幅の 1/4 だけあえて高磁場側に移します。HCCH-COSY ではアミド水素は観えませんので、4.7ppm から低磁場側はほとんど必要なくなるのです。その代わり環電流シフトを受けたメチル基なども検出するために高磁場側までスペクトル幅を広げる必要があります（ピークを折り返してやってもよいのですが、ややこしくなりますので）。

普通はこのような場合、周波数シフト fq を施して 4.7 から 3.0 ppm ぐらいに中心周波数を移せばよいのですが、何故 Br 社はあえてこのような奇妙な方法をとっているのでしょう？しかし、まあここまではよくあるパターンでした。ところが hcchcogp3d をよく見ると、d0=inf1/4 と書かれています。Br 社にしてはめずらしい。これは t1 の初期値を t1/2 に設定することによって、折り返しのピークを負に逆転させる方法です。昔はこれしか使わなかったのですが、もう最近はスペクトル幅を広くとって折り返しをあまり利用しないようになってきました。それで NMRPipe の位相を ph0=-90, ph1=180 に設定しました。Topspin では 90, -180 なので、これも記憶がごちゃごちゃになってしまう要因です。さらに echo-antiecho のプログラムによっては ph0=0 もあり得ます（足し算と引き算の結果のどちらを虚数側に選ぶかで変わってくる）。

ところがスペクトルを見てみると、かなり位相がずれているのです。そもそも NUS データに普通の FT を施して位相だけを観ようとするのですから悲惨です。めちゃめちゃな位相のピークが羅列しているのです。しかも軽水溶媒で測定したので、水のベースラインのうねりもかぶってきて、もうまるでゴミが一杯浮かんだ昔の荒れた大阪湾のような情景です。HCCH-COSY も普通の COSY のようにピークが桜模様になるんだっけ？と思いながら（そんなことはありません）ph0 = 0, 90, ph1=180, -180 さらに F2 の方が狂っているのかもと F2 も同様に変えると、その組み合わせだけで 16 通りです。この間違いに気付くには一晩を要したのですが、実はスペクトル中心を 1/4 だけずらしているので、ph0 も 1/4 だけずらさないといけなかったのです。何の 1/4 か？もちろん ph1=180 の 1/4 です。したがって、ph0= 90+180/4=135, ph1=180 が正しかったのです。NUS でなければ、マニュアルで位相補正した際に何気なく気付くのかもしれませんが、NUS データでは FT 後のデータは錯乱状態ですので、その中でかろうじてピークの姿を醸し出しているのを選んで位相調整することになります。もう二度目の間違いはごめんですので、ここに覚え書きしておいて次は参照することにします。

ところで４次元の NUS のプロセスについてですが、自分の PC で行うと２日経っても終わらず、さらに Word すらも重い状況になってしまいました。しかし、NMRBox を活用してから、その問題が一挙に解決しました。NMRBox は一種のクラウドで、無料で過去のさまざまな NMR ソフトを活用できます。先日もちょっと質問を送ったところ、翌日には回答が返ってきました。もう本当に「ありがとう！」です。おかげで４次元に躊躇しなくなりました。

２匹でやめておきたい時

何気なく今まであやふやにして誤魔化していたことを少し確かめてみました。どのような実験であれ誤差を出さないといけません。英語では uncertainty, standard deviation (sigma), error bar, root mean square deviation (RMSD) などと表します（細かな違いがあるとは思いますが）。ここでは（平均値）+- （１標準偏差）を表すことにします。「プラス・マイナス」にご注意ください。しばしばエラーバーの長さは？と尋ねられ、プラスもマイナスも含めた全体のバーの長さのことなのか？それとも正か負か片方だけの長さのことなのか？で迷います。前者は後者の２倍の長さになります。しかし、以下では後者の方（片方だけの長さの絶対値）を指すことにします。

もし実験を 100 回も 10,000 万回も繰り返すことができれば、それらの値を統計処理して、（平均値）+- （１標準偏差）を出すことができます。これが実験の上ではもっとも正確な誤差でしょう。しかし、一回の測定に云百万円もかかるとすると、せいぜい２回程度（１回？）しかできないことになります。この２回の実験で得られた値を使ってエラーバーを書いてもよいのかどうか？という点が長らくの疑問でした。

ここから下に書く内容は、次のことと多いに関係しています。標準偏差を出す時の手順です。まず、データ a が 10 個あるとすると、その 10 個の平均値 v を出します。そして、それぞれのデータ a から v を引きます。さらにそれぞれの差 (a-v) を２乗します。それらの値 (a-v)^2 を足して n=10 で割ります。これを分散と呼びますが、最後にこの分散のルート（平方根）をとれば標準偏差になります。（標準偏差）= Sqrt（分散）sqrt とは root のことで、しばしば数学のプログラムで使われています。

この処理はちょうど RMSD と同じです。RMSD は後ろから読んでいきます。D は differnece ですので差です。つまり、平均値 v からの差をとります (a-v)。Square は２乗です (a-v)^2。Mean は平均です。R は root です。よって「RMSD = 標準偏差」になりそうな気がします。

ところが、標準偏差の式をいろいろと調べてみると、(a-v)^2 を全て足すところまでは良いのですが、その後 n で割るのではなく (n-1) で割っている式も時々出てきます。昔から不思議だったのですが、この (n-1) はいったい何？

実は Wikipedia に答がちゃんと書かれていました。無限回数の実験をこなした時（いわば理想的な状況）での標準偏差と、事情により数回の実験しかできなかった時（現実的な状況）での標準偏差のことでした。そして、

無限回数での分散 : 限定回数での分散 = n : n-1

となるのだそうです。標準偏差を２乗すると分散になりますので、上記の公式をルート（平方根）すると標準偏差の式にそのまま変身します（しかし、そうとも言い切れないので下記を参照してください）。

ここで現実的な状況を考えてみましょう。分散を計算する時に n で割り算したとします。つまり (R)MSD の計算と同じです。すると、それはあくまでその限定された回数における分散であって、仮に無限回数実験をしたと仮定した時の分散値よりも小さくなってしまうということです。それでは、少ない回数であたかも無限回数実験したかのように計算するにはどうすればよいのか？n で割るのでなく (n-1) で割るということになります。すると、ちょっと大きめの値が出てくると思います。これで良いのだそうです。なぜ (n-1) で割ってちょっと大きめにするかは、調べてみるといろいろ載っていましたが、ここでは複雑ですので記さないことにします。

よって数匹の鼠 n=3 でしか実験しなかった状況で (R)MSD をまじめに計算すると（３で割ると）、もし１億匹の鼠 n=100,000,000 で実験したと仮定した場合の (R)MSD よりも小さめの値が出てしまうということです（n-1=2 で割るとよい）論文にそのようなエラーバーを出すとこれはデータを都合の良い方に解釈したことになり、データ捏造になるのでしょうか？

ややこしい問題は使っているソフトにもあります。もし、そのソフトの標準偏差のツールが、すでにあえて n ではなくて (n-1) を使っていたとすると、その標準偏差は無限回数への拡張を想定していたことになります。一方、RMSD のように n で割っていたとすると、それはその制限回数内のみでの標準偏差を表していることになり、ちょっと小さい値なので得をします。実験を数多く繰り返して n を大きくすると、n と (n-1) の違いが小さくなってきますので、両者の RMSD はどんどん似てきます。

n と n-1 が問題になってくるのは n の値が小さい時です。可哀想な話ですが、もし鼠 100 匹で実験したとすれば、無限回数と制限回数での標準偏差はたったの 0.5% しか違いません。もし、n=10 ですと 5% ぐらいです。幸いうちではバクテリアしか使いませんので、制限されているといえどもすでに数え切れない位の大きな n と言えます。

さて、問題は２回しか実験しなかった場合です。ソフトによってはエラーバーの大きさが 1.4 倍（sqrt(2)）違ってきますので、ちょっと問題です。

２つのデータ a, b があり、それの RMSD を計算しました。v = (a+b)/2 ですので、(a-v)^2 と (b-v)^2 を計算し、これを足して２で割ります。最後に平方根をとります。Mathematica で

Simplify[Sqrt[Mean[{(a - Mean[{a, b}])^2 , (b - Mean[{a, b}])^2}]]]

と打つと、これが abs(a-b)/2 と同じ値であることが分かります。abs とは絶対値のことです。

一方、２個のデータから無限回数を想定して統計をとるとすると、標準偏差は abs(a-b)/sqrt(2) となり、先ほどの２個だけで閉じた標準偏差よりも 1.4 倍だけ大きな値になります。この値をしぶしぶ文献に載せるべきでしょうか？

Simplify[Sqrt[Plus[(a - Mean[{a, b}])^2, (b - Mean[{a, b}])^2]/(2 - 1)]]

では Excel ではどうなっているのだろう？と気になり、調べてみるとちゃんと２個の関数が別々に定義されていました。具体的には 5.3 と 9.2 という２個の数値を使いました。STDEVP 関数は n を使っているので 1.95 となりました。一方 STDEV は n-1 を使っているので 2.76 となりました。両者は 1.4 倍違います。ここでは 1.95 を使いたいところですが、ちゃんとプラスマイナス 2.76 の長さのエラーバー（合計長さは 2.76*2）を記載しました。P は parent の P です。限定された少ない標本を選び出してきて、その少ない数があたかも母集団そのものである（他にデータは何もない）と考えるとよいのだそうです。

本当にこれで良いのかどうかを確かめるべく、シミュレーション実験してみました。乱数が正規分布になるようにしました。例えば（平均値３, 標準偏差２）になるような乱数を 10,000 個発生させました。つまり、10,000 万回実験したことになります。次に任意の２個を抜き取ってきます。これが実際には２回しか実験できなかったことを意味します。そして、この２個の数値を処理した時に果たして（平均値３, 標準偏差２）になるかどうかです。

まず平均値３については、まったく問題なく再現できました。もちろんたまたま選び出した２点の平均値は３からずれています。そかし、それを多数回おこない、それらの平均をとると３に極めて近づきました。しかし、標準偏差２についてはちょっと上記で理解したこととは異なってしまいました。

実際に任意の２点の差の絶対値をたくさん求めてみた結果、その平均値は 2.27 になりました。つまり、ルート２で割る必要がないぐらいです。それではと中央値 median をとってみると 1.93 でした。これもルート２で割るまでもありません。それでは２点の標準偏差はどうなるのかと思い計算してみました。Mathematica の StandardDeviation 関数のですので無限回数を想定しての標準偏差（n-1=1 を使っている）です。その平均値は 1.60 でした。ますます変な値になることは分かっていながら RMSD を計算してみました。n=2 を使う方式です。これは 1.13 と、やはり予想どおり小さな値（1/1.4）になってしまいました。

これはおかしいなと思い、また Wikipedia をしっかりと見てみると答が載っていました。

分散の場合
　無限回数での分散 : 限定回数での分散 = n : n-1
標準偏差の場合
　無限回数での標準偏差 : 限定回数での標準偏差 = sqrt(n) : sqrt(n-1.5)

これはややこしい。単純に平方根をとればよいと思っていたのに、無限回数を想定した分散と標準偏差とでは計算の仕方が違うのだそうです。前者は n-1 で割り、後者は n-1.5 で割ると理想値に近づくのだそうです。n=2 のとき (n-1.5) = 1/2 になります。そのようにして２個の標準偏差を計算すると、ちょうど差の絶対値と同じ値になります。

Simplify[Sqrt[Plus[(a - Mean[{a, b}])^2, (b - Mean[{a, b}])^2]/(1/2)]] は Abs(a-b) と同じです。

以上から言えることは、２点間の差の絶対値をそのまま１標準偏差 s にするとまずまず良いのではないかということです。平均値を v とすると、無限回数測定したと仮定した場合の統計は v +- s となります。n=2 以外の場合、STDEVP 関数で計算した値には sqrt(n/(n-1.5)) をかけ、STDEV 関数で計算した値には sqrt((n-1)/(n-1.5)) をかけて補正すればよいのでしょうか？ややこしいです。

ファンデルワールス相互作用を観る

今までも水素結合 (1999 年), CH/π 相互作用, 静電的相互作用に J-coupling が見つかり NMR で、その相互作用による交差ピークが検出されてきましたが、とうとう蛋白質 GB3 でも van-der-Waals 相互作用が検出されたようです（低分子ではすでに見つかっていた）。

Li J, Wang Y, An L, Chen J, and Yao L. (2018) Direct observation of CH/CH van der Waals interactions in proteins by NMR. J. Am. Chem. Soc. 140, 3194-3197. doi: 10.1021/jacs.7b13345.

NMR 測定では 2D 1H-13C HMQC を利用しています。メチル基の 1H から 13C1 に磁化を移動させて 2C1yHx を作った後 125ms (2T) かけて J(CC)-coupling を展開させます。すると、2C2zC1x が少しできます（Hx がずっと付いてくるのですが、ややこしいので以下では省略します）。そこに位相 y で 90 度パルスを 13C にかけると 2C2xC1z となります。最初の 13C1 から J-coupling した 13C2 へと磁化が移動したわけです。いわゆる 13C-13C の COSY です。その２個目の 13C2 の化学シフトを t1 で検出した後、もと来た道を戻ります。そして１個目の 13C1 に付いた 1H の FID を検出します。試料はメチル基のみが 1H 化、それ以外は 2H 化されています。

水素結合での J-coupling の検出（3D HNCO）においてもそうですが、ピークを少しでも観ようとすると、とにかく J による展開時間 2T を長くとならないといけません。そのため大きめの蛋白質になると、今回の場合は 13C の横緩和のために急に信号が観えなくなってしまいます。今回、横磁化になるのはメチル基の 13C ですので、交差相関緩和（methyl-TROSY 効果）によって横緩和時間がかなり長くなっています。

行きの磁化移動を三角関数を使って表現すると以下のようになります。

交差ピーク：C1y　→　2C2zC1x * sin(pi*J*2T)

対角ピーク：C1y　→　C1y * cos(pi*J*2T)

交差ピークと対角ピークの比をとると tan^2(pi*J*2T) となります。tan が２乗になっているのは行きと帰りの磁化移動が二重にかかってくるためです。そして、実際のピークの強度比をこの式に当てはめると J-coupling の値は 0.1-0.5 Hz ぐらいの大きさになったそうです。2T = 150 ms ほど取ったとしても tan^2(pi*J*2T) =0.06 ぐらいの小ささです。

どうもこのような測定を見ていると、もしかして 13C と 13C の間の NOE がアーティファクトとして入って来てはいないだろうかと勘ぐってしまいます。もちろん NOE が起こるためには 13C 磁化が縦 z 方向にないといけないのですが、いずれのパルスも完全とは言えませんので、どうしても縦磁化が少し混じってきてしまいます。それらは位相回しやグラジエントできっちりと消してやらないと、13C-13C NOE が観えてしまいます。また、1H-1H NOE は重水素化している蛋白質では起こらないはずですが、これもパルスの不完全性と重水素化の不完全性から生じてしまうこともあります。上記のように 2T を変えた時にピークの強度比が tan^2 の曲線に載ったので、おそらく J-coupling によるものと思われます。

J-coupling は隣り合う原子間でそれらの電子雲が重なっていると生じます。量子力学計算によると、CH3-CH3 の方が CH-CH に比べて炭素間の距離は短くなるようです。CH-CH の場合は C-H-H-C と一直線になるのに対して（直線配置）、CH3-CH3 では両者からの水素原子が組み合わさったような構造（交差配置）をとるためのようです。DFT 計算によると、この J-coupling の大きさは C-C 間距離が狭くなるほど大きくなるようです。蛋白質はできるだけパッキングした構造をとろうとするので、CH3-CH3 の交差配置の構造の方を採りやすく、そのため J(CC)-coupling も大きく出そうな気がします。しかし、同じ距離ならば、C-H---H-C と直線に並ぶほど J は大きくなるのだそうです。その方が電子雲の重なりが大きくなるためですが、重なりが大きいほどお互いに反発し合って蛋白質の構造の点では不安定になります（交換反発エネルギーが大きくなる）。これらの特徴のために、今回の感度では全体の 40% ほどしか J が観えなかったのかもしれないと著者らは分析しています（蛋白質は電子軌道の重なりを避けて反発力が小さい最密充填構造をできるだけ採ろうとするが、その結果、小さな J-coupling となってしまう）。このように距離や角度に依存する点は水素結合での J-coupling の特徴と似ています。

ほどけたままくっ付くらしい

正電荷をたくさん持ったある蛋白質と負電荷をたくさん持ったある蛋白質とが、いずれも決まった構造をもたずに、お互いに相互作用していたというお話です。

Borgia, A., Borgia, M.B., Bugge, K., Kissling, V.M., Heidarsson, P.O., Fernandes, C.B., Sottini, A., Soranno, A., Buholzer, K.J., Nettels, D., Kragelund, B.B., Best, R.B., and Schuler, B. (2018) Extreme disorder in an ultrahigh-affinity protein complex. Nature 555(7694), 61-66. doi: 10.1038/nature25762.

細胞（核）内の雑多な環境内にあるにもかかわらず、その二つは相手を見つけて？遠くからどうしでも相互作用します。両者ともに intrinsically disordered proteins（IDP）です。普通は相互作用した途端に induced fit 等が起こって決まった特定の立体構造に固まり、いわゆる鍵と鍵穴のような相補的な関係によって噛み合うものと期待します。しかし、今回の例では相互作用してもお互いに特定のアミノ酸どうしで組み合うわけでもなく、両者がフレキシブルに動いたダイナミックな状態で相互作用していたということです。特定の鍵と鍵穴の間の立体構造の関係によらない、このような相互作用でしたので、Nature に採り上げられたようです。

正電荷の蛋白質：linker histone H1.0 (H1) +53

負電荷の蛋白質：prothymosin-α (Pro-Tα) -44

しかし、それでもあの雑多な中で相手を間違えないというのは驚きです。普通は鍵と鍵穴の相補的な関係を通して「特異的に」相互作用します。もし、両者の立体構造が噛み合わなければすぐに離れてしまい、複合体は相互作用と呼べるほどの滞在時間を保てません（koff が大きい）。もちろん、この二つの蛋白質の電荷の絶対値は両者ともにかなり大きいですので、強い静電的相互作用によりお互いに離れないのでしょう（pM の解離定数をもつ）。また、静電的相互作用は疎水的相互作用よりも遠くまで及ぶので、この２つの蛋白質が少しぐらい離れていてもお互いに引き合うのでしょう。

ただし、本文には次のように書かれています。このように決まった構造をとらずに相互作用する目的は、非常に強い親和性を持ちながら、それでもなお速く付いたり離れたりする必要があるためと。たしかに Pro-Tα はリンカーヒストン H1 のシャペロンですので、H1 がクロマチンと相互作用する時の親和性に競合するほどの親和性を H1 に対して持っていなければなりません。しかし、Pro-Tα がずっと H1 に付いたままですと、H1 をクロマチンから離して次の場所に再配置するといった「交換」を行えません。遠くからでも効く高い親和性をもちつつ、高速に付いたり離れたりするには、このような静電的相互作用と同時に disorder が必要なのだと書かれています。この考察は個人的にはちょっとよく分かりません。それよりかは、クロマチンも含めた３者複合体を作った時に初めて induced-fit が起こるのではないだろうか？第三者とはクロマチン、修飾酵素などが考えられます。例えば Pro-Tα:H1 には複数の複合体の形があり、次にやって来る第三者の（鍵）構造に応じて、そのうちのどれかの（鍵穴）構造が選ばれる（population selection）。しかし、NMR も含め現在の物理的計測法では、この複数の、しかも速く入れ替わっている構造を区別してとらえることは難しく、全ての構造の相加平均として計測されてしまっているのではないだろうか？と想像してしまったりもします（CD の値や 13Ca の chemical shift deviation が複合体の形成の前後であまり変わっていないので、やはり disorder したままなのか？あるいは、大半は disorder していても、実は第三者が来た時にピタッと当てはまる鍵穴構造がいくつか含まれているのかもしれません。そのモル比が小さいので、相加平均をとると disorder に見えてしまう？）。

細胞内には他にも電荷をたくさんもった蛋白質や核酸成分があり、それらが無差別に（非特異的に）くっ付いて来て離れなくなってしまっても良さそうです。もちろん、教科書にも載っているように、核内で局所的に存在したり、ある時期にだけ特別に発現したりしているのかもしれません。これについては、同じ号の

Tight complexes from disordered proteins (NEWS AND VIEWS) Berlow, R.B. and Wright, P.E.

に次のようにコメントされています。

"Pro-Tα and H1 form an archetypal fuzzy complex that involves a large ensemble of possible bound protein conformations, many of which are adopted by only a small number of individual complexes and occur with approximately equal probability."

相補的に噛み合った立体構造の箇所（とまでは明記されてはいませんが）が存在するのだが、ある瞬間をみると、そのフィットし合った領域はたいへん狭く、その領域は時間が経つと共にあちこちに移ります。H1 のある一つの狭い領域が Pro-Tα のある一つの狭い領域といつも１：１で相互作用しているのではなく、 H1 の複数の狭い領域と Pro-Tα の複数の狭い領域とが、いろいろな組み合わせでお互い交代しあいながら相互作用しているのかもしれません。しかもその小さな領域どうしが結び付いている時間は全てが同じぐらい短いのでしょう。すると、それらの複合体の構造がお互いに素早く交換していくので、観察者にはフレキシブルに動きながら非特異的に相互作用しているように観えるのかもしれません。そして、常に余った正電荷が H1 にあるために（Pro-Tα の負電荷によって H1 全ての正電荷が同時に相殺されるわけではない）、H1 がクロマチン（負電荷）とも相互作用できるのかもしれないと書かれています。

核内にある仁（核小体）は、何か膜のようなもので囲まれているわけではなく、蛋白質と核酸が相互作用し合いながら寄り集まっています。これは液-液相分離と呼ばれていますが、このような相互作用も上記のような仕組みで起こっているのかもしれません。ただし、今回の H1:Pro-Tα の系では濃縮しても液-液相分離は起こっていません。著者らは、ちょうど長さの点でも電荷の点でも２者間の相互作用でお互いに打ち消し合うような関係にあるためではないかと推察しています。また、もし疎水的な残基や芳香環がもっとあれば、これにさらに疎水的相互作用やカチオンパイ相互作用などが加わり、液-液相分離に進んでいくのかもしれないと書いています。

そもそも２つの蛋白質はお互い遊離状態であっても複合体状態であっても CD スペクトルがほとんど変わりませんでした。ということは、二次構造（特に α ヘリックス）の量は、複合体になったからといって増えているわけではありません。これ以上の構造情報を得ようとすると（結晶は当然のように出来ませんので）NMR しかありません。しかも二次元 1H-15N HSQC だけでかなりの情報が得られます。大腸菌発現系であれば、15N 標識蛋白質の発現が比較的容易にできます。IDP の二次元 1H-15N HSQC スペクトルは特徴的でして、横軸（1H）の 8.6 ppm より左にピークが現れません。これは「ハムの壁」と呼ばれています（この語呂は日本語だけで成り立ちます）。H1 には一部 fold した領域がありますので、そこのアミノ酸のアミド 1HN からのピークのみ 6-12 ppm 範囲に散らばります。一般的に水素結合が強いほど大きな 1H 化学シフト値（スペクトルでは左側, 低磁場側）をとる傾向があります（NMR では、軸の向きが常識とは逆になっていることに注意）。

ここで、２つの蛋白質を混ぜ合わせると、両者ともにピークが移動しました。一般的に電荷や芳香環をもった相手方リガンドが近づいてくると、化学シフト値が大きく変わります。しかし、この混合で移動したピークが一面に散らばる方向に動けば、複合体を形成した際にはっきりとした立体構造が誘導された（induced-fit）と分かりますが、移動したピークも依然せまい領域に固まっていました。つまり、二次構造の量が増えたわけでもないのです。また、13C の化学シフト値は二面角の大きさに強く影響を受けますので、二次構造の判定にしばしば使われます。しかし、その 13Ca の化学シフト値を見ても、α ヘリックスや β シートが誘起されたようではありませんでした。

このように帰属をしなくてもかなり詳しい構造情報が NMR から得られます。もし、15N, 13C で二重標識できれば、それぞれのピークがどのアミノ酸由来であるかを帰属できますので、もっと細かく調べることができます。実は、IDP の帰属は技を使えば可能な場合が多いのです。なぜならば、主鎖、側鎖が揺れ動いているために見かけの分子量が小さくなったような効果が生じ、感度が著しく上がるためです。複合体のスペクトルで、H1 の globular-domain からのピークが消えたのは、まさに溶液内での回転運動が遅くなったためです。しかし、それでもフレキシブルな、構造をとっていない領域は観測が可能でした。もちろんピークは全体的に強度が下がったそうです。さすがに複合体の状態では、単量体でいる時よりも回転運動が遅くなったのでしょう。また、反対の電荷どうしの相互作用が速く入れ替わったので、Rex も大きくなったのでしょう。

「技」というのは高次元化や高分解能化です。つまり、多くのピークがスペクトルの狭い領域にひしめき合うので、分解能を上げた測定をしてあげる必要があります（普通の蛋白質と同じパラメータで測定すると失敗します。それこそ non-uniform sampling, NUS が有効でしょう。さらに帰属の作業には少しばかりの根気が必要）。しかし、感度さえあれば、この高分解能化を達成する方法はいろいろと紹介されており、むしろ NMR の超複雑な最新パルス技術とプロセス技術を存分に発揮できますので、オタクにとっては実は嬉しい系であったりもします。

渡り鳥

渡り鳥がどうして道を間違えずに地球を半周ほども回れるのか？という話についてです。当然のように方位磁石（コンパス）*1 の代わりになる何かを持っているのでしょうが、それが何でどのような物理的原理によるのか、よく分かっていません。個人的に興味を持っていて、このような記事をたびたび拾い読みします。もっともよく登場するのは、網膜にある cryptochrome（クリプトクローム）と呼ばれる蛋白質です。

*1) ある種の細菌は磁鉄鉱をもっており、実際にそれが方位磁石のような働きをするそうです。鳥類にも磁鉄鉱が見つかるのですが、コンパスの働きはないとされています。

今回は驚いた事に、DNA 修復酵素の一種である光回復酵素が、このコンパスの働きをしているかもしれないという記事が出ました。

Zwang, T. J., Tse, E. C. N., Zhong, D. P., and Barton, J. K. (2018) A compass at weak magnetic fields using thymine dimer repair. ACS Cent. Sci. 2018, DOI: 10.1021/acscentsci.8b00008

この記事は下記にも紹介されています。

P. J. Hore (2018) Sensitivity of DNA repair enzymes to weak magnetic fields may have relevance to the mechanism by which birds sense the Earth’s magnetic field. ACS Cent. Sci., DOI: 10.1021/acscentsci.8b00091

Hore さんといえば、NMR でも有名な先生ですが、ついでに下記も紹介しておきます。

NMR入門: 必須ツール 基礎の基礎 (Chemistry Primer Series) 

P.J. Hore (著),‎ 岩下 孝 (翻訳),‎ 大井 高 (翻訳),‎ 楠見 武徳 (翻訳)

DNA は紫外線を受けるとしばしば損傷します。その中でもよく知られている損傷がチミンダイマーです。隣り合うチミンどうしが結合して二量体になってしまうのです。チミンはピリミジン環をもつので、ピリミジンダイマーとも呼びます。

そこで、これを修復する酵素 photolyase（光回復酵素）が登場します。この photolyase は内部に FAD を持っています。これが完全に還元された形が FADH-（FADH-**）です。これに青色の光が当たると励起されて、チミンダイマーに電子を１個与え、その二量体を壊します。その際に （FADH*）側と、チミンダイマー（TT-*）側のそれぞれにラジカルができます。このラジカルペアですが、これが一重項状態（αβ-βα）になると緩和時間が伸び、さらに周りの磁場の向きによって DNA 修復の速度が変わるらしいのです。しかも、その磁場は地磁気ほど小さくても良いのだそうです。

詳しいことはよく分かりませんが、このラジカルペアは、内部では核スピンと hyperfine 相互作用を持ち、外部とは磁場との Zeeman 相互作用を持つため、singlet（αβ-βα）と triplet（αβ+βα, αα, ββ）の間で交換します。その最終的な割合が外部磁場との相対角度に依存するらしいのです。これが効率よく起こるためには、緩和時間が長くないといけませんし、また磁場も非常に強くないといけません。そこがまだあまり解明されていない問題点のようです。

実は cryptochrome は photolyase の先祖とされていますので、あながち急に photolyase が飛び出してきたわけではありません。Cryptochrome も FAD を持っていて磁場の向きに応じてラジカルペアを生成します。そこで、これが磁気コンパスではないかと言われています。

しかし、普通はこのような酵素は細胞内（核内）でブラウン運動により回転するため、磁場に対する向きもランダムに動いてしまいます。したがって、まだまだ解明されたとは言えない状態のようです。

渡り鳥だけでなく、鯨, Kujira、鮭, sake、鰻, MagRO、海老蟹 ロブスター、蝶々など、いろいろな生き物が長い距離を旅します。これらがどのようなコンパスを持っているのか、まだよく分かっていないようです。磁気、海水の香りなどさまざまな説があります。そして、もしかして人間も？と思ってしまいます。

ちょっと前までそれはまずあり得ないと思っていました。ところが、先日、とある遮蔽されていない超高磁場 NMR の下に潜って頭を動かすと、まるで車酔いしたような気分になりました。磁石の下から這い出すと、同じ姿勢でも全く問題ありません。他の数人もいっしょに何度試しても皆そのようになるので、もしかして何かある！と感じました。また、NMR 室に見学者を招待すると 100 人に１人ぐらいの割合で気分が悪くなる人が出てきます。毎年そのような事態になるので、いつも「閉所恐怖症ですか？」と尋ねるのですが、そうでもなく、その人達もたいへん不思議だと答えます。そういえば、目隠ししていても方角が分かる人がテレビで紹介されていました。どうなっているのでしょう？一応、血液の中にはヘモグロビンがあり、そのヘム鉄により強い磁場の中ではヘモグロビンは磁場に対して配向します。しかし、血流が乱す力の方が配向よりも圧倒的に大きいので、このまるで residual dipolar coupling, RDC を測る時の現象が、方向感知に効いているとは考えにくいでしょう。

ところで、上記のラジカルペアについてですが、これができる前は、この２つの電子はもともとは同じ核に所属しています。そのため、一方の電子スピンが上向き（α 状態）にあれば、他方の電子スピンは下向き（β 状態）にあります。ちょっとややこしいことに、一方が必ず α 状態にあるという意味ではなく、本当は α 状態と β 状態の両方の状態に同時にあります（α or β ではなく α and β）。これを量子力学の重ね合わせ状態と呼びます。そして、どちらの状態にあるのかを知るために観測すると、そのとたんにどちらかに収縮します。そして、もし β 状態に収縮して観測されたならば、相棒の電子スピンは即座に α 状態に収縮します。

興味深いことに、一つの電子が核を離れて別の核に移動し、ラジカルペアに変身しても、まだこの関係が保たれる場合があります。奇妙な量子もつれと呼ばれるそうです。もちろんそのようなコヒーレンスが保たれる時間が問題ですが、仮にそのようなもつれ状態がずっと続くとすると、なお不思議なことに、原理上は２つの電子がどれだけ離れていても（宇宙の両端に離して置いても）この関係が保たれるのだそうです。その場合、一方の電子スピンを観測して α 状態に収縮すると、遠く離れた電子スピンにも遠隔作用が及んで β 状態に収縮します。

しかし、ラジカルペアの一重項状態に三重項状態がある比率で混ざってくると、αβ 状態だけでなく αα や ββ 状態も混ざってきます。それらの比率がラジカルペアの向きと磁場の向きとの間の相対角度で決まるのだそうです。そして、渡り鳥がそれぞれの電子スピンの状態を観測（認識）することにより、一重項と三重項状態の間の比率が分かり、しいては磁場の向きが分かるという仕組みなのだそうです（ちょっと理解に自信がありませんが）。

どうも物理の話が濃すぎて、何が何だかよく分かりません。さらに最近は多世界解釈という思想？も入ってきました。これによると、αβ 状態と βα 状態という２つの世界があり、どちらかを観測したに過ぎないのだそうです。つまり、一方を観測した際にたまたま α 状態に収縮したので、その瞬間に他方は β 状態に収縮したという解釈ではなく、たまたま αβ 状態の世界の方を観測したに過ぎないという解釈です（あるいは観測した途端に世界が αβ と βα に分離した？）。どちらも嘘みたいな話ですが、個人的にはどちらかというと、この多世界解釈の方がしっくり来ます。

しかし、渡り鳥の目の中で本当にそのような事が起こっているのでしょうか？実は cryptochrome は植物も含めさまざまな生物が共通して持っており、実際に体内時計に関与することはよく知られています。よって、渡り鳥のような高等生物に進化する過程で生み出されてきたのではなく、はるか何億年も昔、生物の初期段階にすでにあったのではないかと言われています。