今だに位相補正で苦しむ

何でもない 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)) をかけて補正すればよいのでしょうか？ややこしいです。