Hatena::ブログ(Diary)

biochem_fan's note

リハビリを兼ねて、興味深く思ったことや考えたことを書いていこうと思います。
何もかもつまらなく感じて辛くなることがあっても、今、面白いと思ったことは本当なので。
持病により、時々「鬱発言」がありますのでご注意を。詳細は 精神不調をご覧ください。
カテゴリー

2099-12-31

はじめに

興味深く思ったことを書いていこうと思う。何もかもつまらなく感じて辛くなることがあっても、今、面白いと思ったことは本当なので。その一方で、持病が悪い相にある時は、愚痴っぽいことを書いていることも多い。そんな恥をあえて晒している理由は、精神不調に書いた。

全角のシャープ(#)は、もっと書くべき・要検討・あるいは直したいが、とりあえず公開してしまう場合につける目印である。本来#がなくなってから公開すべきであるが、完成度100%を待っていると下書きが増えるばかりで、吐き出せるものも吐き出せなくなるので。特に、しばらくは、下書きに溜まった大量の原稿を処分することから始めたい。

学問・技術的な話については、twilog を発掘していただいたほうが良いかと。

2018-08-26

ベクトルを転置したものの積について

自明なことだけれど、こういう表記でつまづく人がいるのはもったいないので、ここに明記しておく。

 ¥mathbf{v} = ¥begin{pmatrix}v_1 ¥¥ v_2 ¥¥ ¥vdots ¥¥ v_n ¥end{pmatrix},  ¥mathbf{w} = ¥begin{pmatrix}w_1 ¥¥ w_2 ¥¥ ¥vdots ¥¥ w_n ¥end{pmatrix} とするとき、
 ¥mathbf{v}^T¥mathbf{w} = ¥begin{pmatrix}v_1 && v_2 && ¥cdots && v_n ¥end{pmatrix} ¥begin{pmatrix}w_1 ¥¥ w_2 ¥¥ ¥vdots ¥¥ w_n ¥end{pmatrix}= v_1 w_2 + v_2 w_2 + ¥cdots + v_n w_n
内積である。一方、
 ¥mathbf{v}¥mathbf{w}^T = ¥begin{pmatrix}v_1 ¥¥ v_2 ¥¥ ¥vdots ¥¥ v_n ¥end{pmatrix}  ¥begin{pmatrix}w_1 && w_2 && ¥cdots && w_n ¥end{pmatrix} = ¥begin{pmatrix}v_1 w_1 && v_1 w_2 && ¥cdots && v_1 w_n ¥¥ v_2 w_1 && v_2 w_2 && ¥cdots && v_2 w_n  ¥¥ ¥vdots ¥¥ v_n w_1 && v_n w_2 && ¥cdots && v_n w_n  ¥end{pmatrix}
は要素ずつの積を並べた行列である。
この行列の各列は  ¥mathbf{v} = ¥begin{pmatrix}v_1 ¥¥ v_2 ¥¥ ¥vdots ¥¥ v_n ¥end{pmatrix} w_i をかけたものに過ぎないから、お互いに定数倍である。したがって、その rank は 1 である。各行についても同じことがいえる。

2017-09-02

実装して学ぶデータ処理システム その3 - Kabsch transform の意味

今回は XDS における最大の特徴、Kabsch transform の意味について説明する。論文の 2.3 節に対応する。

Kabsch transform は、goniostat 座標系から reflection 座標系への変換である。これまでに出てきた検出器系や実験室系との変換と違い、この変換は局所的である。局所的というのは、1つ1つの逆格子点ごとに固有の reflection 座標系があって、その点の近傍でしか変換が意味をなさないということだ。数学が得意な人は、多様体における局所座標系の貼りあわせをイメージするかもしれない。ただ、検出器上の1つのピクセルは、最近傍の反射がもつ局所座標系だけに属することになっており(pixel labeling)、局所座標系同士が重なりあうことはないから、多様体よりも単純である。

ある反射 hkl の固有座標系は、逆空間(goniostat 座標系で考えるのがよい)でその反射に対応する格子点を原点とし、基底ベクトル

 ¥mathbf{e_1} = ¥frac{¥mathbf{S} ¥times ¥mathbf{S_0}}{|¥mathbf{S} ¥times ¥mathbf{S_0}|}
 ¥mathbf{e_2} = ¥frac{¥mathbf{S} ¥times ¥mathbf{e_1}}{|¥mathbf{S} ¥times ¥mathbf{e_1}|}
 ¥mathbf{e_3} = ¥frac{¥mathbf{S} + ¥mathbf{S_0}}{|¥mathbf{S} + ¥mathbf{S_0}|}

によって定義される。1つ目の式は、¥mathbf{e_1}¥mathbf{S}, ¥mathbf{S_0}直交する長さ1のベクトルだと言っている。二つ目の式は 、こうしてできた ¥mathbf{e_1}¥mathbf{S}直交する長さ1のベクトル ¥mathbf{e_2} だと言っている。どちらも散乱ベクトル ¥mathbf{S}直交するから、この2つのベクトルは Ewald 球の接平面上にある。¥mathbf{e_3} の向きを考えるために、¥mathbf{e_1} との内積を取ると、

 ¥mathbf{e_1} ¥cdot ¥mathbf{e_3} = ¥frac{(¥mathbf{S} ¥times ¥mathbf{S_0}) ¥cdot (¥mathbf{S} + ¥mathbf{S_0})}{|¥mathbf{S} ¥times ¥mathbf{S_0}||¥mathbf{S} + ¥mathbf{S_0}|} = ¥frac{(¥mathbf{S} ¥times ¥mathbf{S_0}) ¥cdot ¥mathbf{S} + (¥mathbf{S} ¥times ¥mathbf{S_0}) ¥cdot ¥mathbf{S_0}}{|¥mathbf{S} ¥times ¥mathbf{S_0}||¥mathbf{S} + ¥mathbf{S_0}|} = 0

となるから、¥mathbf{e_1}¥mathbf{e_3}直交している。同様に ¥mathbf{p^*} = ¥mathbf{S} - ¥mathbf{S_0} との内積をとっても、

 ¥mathbf{p^*} ¥cdot ¥mathbf{e_3} = ¥frac{(¥mathbf{S} - ¥mathbf{S_0}) ¥cdot (¥mathbf{S} + ¥mathbf{S_0})}{|¥mathbf{S} + ¥mathbf{S_0}|} = ¥frac{¥mathbf{S} ¥cdot ¥mathbf{S} + ¥mathbf{S} ¥cdot ¥mathbf{S_0} - ¥mathbf{S_0} ¥cdot ¥mathbf{S} - ¥mathbf{S_0} ¥cdot ¥mathbf{S_0}}{|¥mathbf{S} + ¥mathbf{S_0}|} = 0

となって、こちらも直交している。ここで、内積が可換であることと、|¥mathbf{S}| = |¥mathbf{S_0}| (Laue 条件) を用いた。

反射座標系の成分のスケールとしては、角度(度)を用いることになっている。長さの単位として角度というのは不思議に感じられるかもしれない。これが可能となるのは、¥mathbf{e_1}¥mathbf{e_2} は Ewald 球の接平面上にあり、ごく近傍では Ewald 球の表面(曲面)とも同一視できるからである。Ewald 球面上の弧の長さは、半径 |¥mathbf{S_0}| = ¥frac{1}{¥lambda}と中心角(ラジアン)をかけたものだから、長さの基準として角度を用いても問題ないのだ。

数式はこのくらいにしておいて、意味について考えよう。検出器上で観察される spot profile は、still shot の場合、次のような影響をうける。

  • Ewald 球の曲率に由来する広がり (高角のほうが大きい)
  • スペクトル幅による広がり (高角のほうが大きい)
  • sensor 厚みによる位置のズレ

振動写真の場合は、さらに次のような影響も加わる。

  • 回転軸から遠い( ¥rho が大きい)反射は、大きい線速度を持つ

これらを理解するために、まず、逆空間における1つ1つの「反射球」は、一定の大きさを持っていることを思い出そう。1つの単位胞フーリエ変換したもの、つまり 連続的な molecular transform に無限サイズの格子に由来する逆格子を掛け算したものが理想結晶の逆空間を与える。実際の結晶は有限のサイズを持つから、その shape transform を逆格子点の各点に畳み込んだものとなる。詳しくは◯◯を参照。

一方、検出器上で観測される現象は、#未稿

実装して学ぶデータ処理システム その4 - 座標変換の実際

前回(本記事公開時には未公開)、Kabsch 変換の意味について説明した。今回は、実装上の問題点について述べる。DIALS のレポジトリにある James の文章 が参考になる。

逆空間の一点とKabsch 変換後の反射座標系での一点は、数学的には一対一に対応している。しかし、プログラムとして実装する上では、それぞれを有限の大きさを持ったボクセルの集合として表現するので、多対多の関係になってしまう。つまり、反射座標系の 1 つのボクセルに対して、複数のフレーム上の複数のピクセルが対応しうる。逆に、あるフレームのあるピクセルが、反射座標系では複数のボクセルにまたがることもある。

このような複雑な対応関係を処理するために、XDS は複数の戦略を取ってきた。Kabsch, Wolfgang. "Evaluation of single-crystal X-ray diffraction data from a position-sensitive detector" Journal of Applied Crystallography 21.6 (1988) で述べられている初期のバージョンでは、1 つの検出器ピクセルの値は、反射座標系で近傍の 2 * 2 * 2 = 8 ボクセルに線形補間で配られる。たとえば、検出器ピクセルの中心座標が、プロファイル座標系で (3.2, 2.4, 1.5) に対応するならば、(3, 2, 1) から (4, 3, 2) までの 8 つのボクセルに配るのである。この方法は、反射座標系の grid の刻みが細かくなって、1つの検出器ピクセルがもっと多くのボクセルにまたがるようになると破綻してしまう。

2010 年の論文 Kabsch, Wolfgang. "Integration, scaling, space-group assignment and post-refinement" Acta Crystallographica Section D (2010): 133-144 で説明されているように、新しいバージョンの XDS では、1 つの検出器ピクセルを 5 * 5 = 25 個のサブピクセルに分割して、それぞれを対応するボクセルへ配るようになった。ただし、各ボクセルの間で補間することはしない。この方法でも、1 つの検出器ピクセルが 25 個以上のボクセルにまたがるとうまくいかないはずだが、実用上は問題ないようである(BEAM_DIVERGENCE などをうまく取っている?)。

理論的には、1 つの検出器ピクセルは、反射座標系で凸多面体に対応する。この多面体と交差するボクセルを列挙し、交差する体積を解析的に求めることは、複雑だが可能である。DIALS における XDS アルゴリズムの実装は、この方向で進められていると聞いている。#対応するコードと資料を要確認

XDS では、プロファイル座標系でのボクセル数は NUMBER_OF_PROFILE_GRID_POINTS_ALONG_ALPHA/BETA, NUMBER_OF_PROFILE_GRID_POINTS_ALONG_GAMMA パラメータで設定される。デフォルトは 9 で、最大は 21 だ。これがプロファイル空間でどれだけの幅にマップされるかを決めるのが BEAM_DIVERGENCE と REFLECTING_RANGE である。 BEAM_DIVERGENCE_E.S.D と REFLECTING_RANGE_E.S.D は、三次元正規分布で表される解析的なプロファイルの幅(標準偏差)であり、強い反射のスポット形状から決定される。両側に 3 sigma 取れば正規分布の 99 % をカバーするから、CUT が 2 の場合、 BEAM_DIVERGENCE と REFLECTING_RANGE は 各 E.S.D. の 6 倍あれば十分である。実際には、正規分布よりも広がりがある(over dispersion) かもしれないから、少し多めにとっておいたほうが無難である。実際の XDS がどのような基準で決めているかは論文では明記されていない(邪悪!)が、実は 7 倍らしい。

2017-04-28

RELION の classification で、クラスの割合を抽出する

RELION の 2D / 3D classification で、iteration ごとに各クラスの割合がどう変化したかをプロットしたいと言われたので、スクリプトを作った。

for i in `seq 0 20`; do 
 j=`printf %03d $i`;
 echo -n $j" ";
 relion_star_printtable run_it${j}_model.star data_model_classes _rlnClassDistribution | tr '\n' ' ';
 echo;
done > dist.dat

relion_star_printtable がクラスごとに改行してしまうので、tr を使って置換するなど、ちょっとダサいことになっている。

こうすると、run_it_0NN_model.star (NN = 00 .. 20) から、割合を切り出して行に並べたデータファイル dist.dat

000 0.200000 0.200000 0.200000 0.200000 0.200000
001 0.200000 0.200000 0.200000 0.200000 0.200000
002 0.199695 0.201143 0.198778 0.200882 0.199502
003 0.214579 0.171917 0.207016 0.180553 0.225935
004 0.228429 0.149543 0.176681 0.134442 0.310905
005 0.105184 0.321839 0.107021 0.308757 0.157199 
(後略)

ができるので、gnuplot

plot for [i=2:6] 'dist.dat' using 1:i with line title "Class ".(i-1)

などとしてプロットする。for を使うには、比較的新しいバージョンが必要。
4.6 ではできたが、4.2 ではダメだった。

なお、ファイル数 (20) やクラス数 (5) を指定しなくて済むようにするのを、読者の練習課題とする。

2017-03-11

Jupyter notebook で cctbx などを使う

DIALS を bootstrap.py からビルドしたあと、base/bin/pip を使って Jupyter notebook をインストールしても cctbx を import できない。環境変数を設定する wrapper となっている dials.python などと違い、PYTHONPATH や LD_LIBRARY_PATH が設定されていないようだ。

次のように設定して解決した。

jupyter notebook --generate-config

として、設定ファイル ~/.jupyter/jupyter_notebook_config.py を生成し、そこに以下のような記述を付け足す。何を足したらいいかは、dials.python の中で当該環境変数を調べ、Jupyter notebook 上での値と比べれば分かる。

import os
os.environ['PYTHONPATH'] = "/opt/dials/modules/cctbx_project:/opt/dials/modules:/opt/dials/modules/cctbx_project/boost_adaptbx:/opt/dials/modules/cctbx_project/libtbx/pythonpath:/opt/dials/build/lib:" + os.environ['PYTHONPATH']
os.environ['LD_LIBRARY_PATH'] = "/opt/dials/build/lib:/opt/dials/base/lib64:/opt/dials/base/lib:" + os.environ['LD_LIBRARY_PATH']