「RELION3チュートリアル」の版間の差分

提供: Eospedia
移動: 案内検索
(7.2 Analysing the results)
(7.2 Analysing the results)
行2,627: 行2,627:
  
 
[[File:rln-ctfrefine-03.png|x500px]]
 
[[File:rln-ctfrefine-03.png|x500px]]
 +
 +
 +
ビーム傾斜量の推定値は CtfRefine/job024/beamtilt_0.txt に記録されています。(※beamtilt_0.txtと書いてありますが beamtilt_iter-fit_class_0.txtおよびbeamtilt_lin-fit_class_0.txtのことでしょうか?あるいは本家チュートリアルはRELION 3.0のものなので、古いバージョンの挙動でしょうか。)
 +
 +
※ iter-fitとlin-fitはどちらもビーム傾斜量の推定値が記録されていますが、iter-fitの方が最終的な推定結果だと思われます。ソースコードを見るとまずlin-fitが推定され、それを初期値とした最適化計算によりiter-fitが推定されているようにみえ、かつCTF refinementジョブで生成されるparticle starファイル particles_ctf_refine.starの rlnBeamTiltX, rlnBeamTiltYはiter-fitに記載の推定値とおなじになっています。
 +
 +
※ beamtilt_iter-fit_class_0.txt の中身は、以下のようになっています。ビーム傾斜の推定値はX方向に 0.02 mrad, Y方向に -0.2 mrad です。
 +
<pre>
 +
$ cat CtfRefine/job049/beamtilt_iter-fit_class_0.txt
 +
 +
beamtilt_x = 0.0236249
 +
beamtilt_y = -0.214688
 +
</pre>
 +
 +
 +
また、このビーム傾斜推定に使用された位相ずれを表した画像は以下のコマンドで確認できます。
 +
 +
(※CTF補正済みの参照像の位相と比較したときのズレと思います。ビーム傾斜はコマ収差による位相ズレをモデルでフィッティングすることにより推定されています。位相ずれが最も大きくなる方向がビーム傾斜の方向です。)
 +
 +
<pre>
 +
$ relion_display --i CtfRefine/job049/beamtilt_delta-phase_per-pixel_class_0.mrc --scale 3 &
 +
</pre>
 +
 +
[[File:rln-ctfrefine-04.png|x400px]]
 +
 +
 +
※若干ですが、Y軸方向の下側ほど白く(位相ずれが正に大きい?)、Y軸方向の上側ほど黒く(位相ずれが負に大きい?)なっています。X方向にはそのような傾向は見られません。ビーム傾斜の推定値と対応しています。
 +
 +
 +
フィッティングされたビーム傾斜(位相ずれ)モデルは以下のコマンドで確認できます。(※本家チュートリアルでは lin-fit を表示していますが、最終的なフィッティング結果は iter-fitの方な気がするので、ここではiter-fitを示しておきます。)
 +
 +
<pre>
 +
$ relion_display --i CtfRefine/job049/beamtilt_delta-phase_iter-fit_class_0.mrc --scale 3 &
 +
</pre>
 +
 +
[[File:rln-ctfrefine-05.png|x400px]]
 +
 +
 +
もしビーム傾斜が大きい場合、これらの画像は円の片側が他よりも暗く、反対側が他よりも明るくなります。今回の場合は推定されたビーム傾斜量は X方向に 0 mrad, Y方向に -0.2 mrad でした。今回のデータの分解能では、これはあまり問題にならないくらい小さな値です。ビーム傾斜の推定無しで Ctf refinementをやり直してもいいでしょう。ここでは特にそうしたことはせず、このまま先へ進むことにします。
  
 
==8 Bayesian polishing(ベイジアンポリッシング)==
 
==8 Bayesian polishing(ベイジアンポリッシング)==

2019年6月4日 (火) 03:16時点における版

前書き

  • 本ドキュメントは、RELION3のチュートリアルを日本語訳し、図や補足説明などを追加して充実化したものです(Sjors Scheres博士の許可を得ています)。
  • 誤訳等がある場合には、ご連絡いただけますと幸いです。少しずつ、良いものにしたいと思います。
  • 本ドミュメントは随時更新中です。webブラウザにキャッシュが残っていると古い版の図が表示される場合があるようです。キャッシュを消してから見ていただくと安全です。


実行環境

CTF estimationジョブを除き、1回目のParticle extractionまで

Ubuntu 16.04.6 LTS

  • 3.6 GHz Intel Xeon E5-1650 (6core/12thread)
  • 64 GB 2133 MHz DDR4
  • NVIDIA GeForce GTX 1080 8 GB x 4台
    • Driverバージョン 418.74
    • CUDAバージョン 10.1

RELION version 3.0.5

  • Build from source


※GctfバイナリがCUDA10.1に対応しなかったのと、MPI関係で問題が生じたので、以下の環境に切り替えました。

CTF estimationジョブと、1回目の2D classification以降

CentOS 7 on Docker 18.09.6 with NVIDIA Docker2

  • ハードウェアとRELIONバージョンは上記と同じ
  • CUDAバージョンは 9.2

(参考) DockerでRELION


目次

1 Preprocessing (前処理)

1.1 Getting organized (解析の準備)

まず構造解析プロジェクトのディレクトリを作成します。1つの構造解析プロジェクトにつき1つのディレクトリを作成することを推奨します。これをプロジェクトディレクトリと呼びます。RELION GUIはプロジェクトディレクトリ直下から起動するようにしてください。次に、プロジェクトディレクトリの中に、未処理の電顕画像または動画を格納するためのディレクトリを作成し、そこへ画像または動画を格納します。サポートされているファイル形式はMRCまたはTIFFです。ディレクトリ名は、全ての動画を1つのディレクトリに納めるならばMovies/、例えば撮影日別でいくつかのディレクトリに分けて納めたいならばMovies/15jan16、Movies/23jan16のようにすると良いです。何らかの理由によりRELIONのプロジェクトディレクトリの中に動画を置きたくない場合は、動画を格納したディレクトリのシンボリックリンクをプロジェクトディレクトリ内に作成すれば良いです。

画像ファイル名の拡張子は .mrc、動画ファイル名の拡張子は .mrcs または .tiff とする必要があります。チュートリアル用のテストデータを解凍すると、中にMovies/ディレクトリが存在し、その中にTIFF形式の動画24本、ゲインリファレンスファイル(gain.mrc)、データ収集条件に関する情報が記載されたNOTESファイルが含まれています。

RELION GUIを起動してみましょう。前述のように、GUIは常にプロジェクトディレクトリからの起動が必要です。誤ってプロジェクトディレクトリ外から起動することを防ぐため、GUIは新しいディレクトリで初めて起動された場合、ユーザーに確認を要求します。その為、新しいディレクトリで初めてGUIを起動するときは、"&"を付けてバックグラウンド起動することは避けてください。プロジェクトディレクトリの中にいることを確認し、次のコマンドでGUIを起動します($はコマンドプロンプトを表します)。

$ relion

確認ダイアログに対し"y"と入力すると、RELIONプロジェクトがセットアップされ、GUIが立ち上がります。


RELION3のGUIウィンドウ.png


まずはじめに、動画ファイル達を解析パイプラインにインポートします。ジョブタイプブラウザから”Import”を選択すると、以下の様な入力画面が出てきます。


Importジョブ.png


以下の様にパラメータを指定します。


Importジョブのパラメータ
説明
Input files Movies/*.tiff Movies/以下のすべての.tiffファイル(動画ファイル)をインポートします。
Node type 2D micrograph movies 各ファイルを2次元の画像がスタックされた動画ファイルとしてインポートします。


Current job: Give_alias_here と書かれたフィールドに、このインポートジョブに対するエイリアス名を入力できます(例: movies)。【Run!】ボタンを押すとジョブが実行され、Import/job001 という名前のディレクトリが作成されます。エイリアス名を指定した場合はImport/moviesというシンボリックリンクも作成されます。作成されたディレクトリ内にはSTARファイルと、インポートされた全ての動画ファイルが格納されています。次のコマンドでSTARファイルの中身を確認してみましょう。

$ less Import/job001/movies.star

data_
loop_
_rlnMicrographMovieName
Movies/20170629_00021_frameImage.tiff
Movies/20170629_00022_frameImage.tiff
Movies/20170629_00023_frameImage.tiff
Movies/20170629_00024_frameImage.tiff
Movies/20170629_00025_frameImage.tiff
Movies/20170629_00026_frameImage.tiff
Movies/20170629_00027_frameImage.tiff
Movies/20170629_00028_frameImage.tiff
Movies/20170629_00029_frameImage.tiff
Movies/20170629_00030_frameImage.tiff
Movies/20170629_00031_frameImage.tiff
Movies/20170629_00035_frameImage.tiff
Movies/20170629_00036_frameImage.tiff
Movies/20170629_00037_frameImage.tiff
Movies/20170629_00039_frameImage.tiff
Movies/20170629_00040_frameImage.tiff
Movies/20170629_00042_frameImage.tiff
Movies/20170629_00043_frameImage.tiff
Movies/20170629_00044_frameImage.tiff
Movies/20170629_00045_frameImage.tiff
Movies/20170629_00046_frameImage.tiff
Movies/20170629_00047_frameImage.tiff
Movies/20170629_00048_frameImage.tiff
Movies/20170629_00049_frameImage.tiff

このSTARファイルは動画ファイル名(パス)を記録しただけの単純なものです。無名のデータブロック"data_"があり、"loop_"により _rlnMicrographMovieName をラベルとし動画ファイル名を値とするラベル-値ペアを24個記述したSTARファイルとなっています。RELIONで定義されているラベル名の一覧はこちらをご覧ください→RELIONの定義済みラベル


もしRELION以外のソフトウェアで単粒子切り出しなどを行った場合、以下の前処理ステップを飛ばし、"Import"ジョブタイプを使って単粒子のSTARファイル、3D参照像、3Dマスク等をインポートすることができます。ただし、これはRELIONを使う上で推奨される方法ではなく、適切なSTARファイルの作成はユーザー自身で行う必要があることに注意してください。

1.2 Beam-induced motion correction (電子線誘起性ドリフトの補正)

Motion correctionジョブタイプは、動画のフレームアラインメントプログラムであるUCSF MOTIONCOR2を用いて、電子線誘起性ドリフトを補正します。RELIONバージョン3.0以降はMOTIONCOR2と同じアルゴリズムがRELION自体に実装されており、そちらを使うこともできます。ただしRELION実装版はCPUでのみ動作します。

Motion correctionジョブタイプの設定はI/O、Motion、Runningの3つのタブに分かれています。それぞれパラメータを入力した状態の画面を以下に示します。

MotionCorrection IOtab.png

MotionCorrection Motiontab.png

MotionCorrection Runningtab.png


パラメータの一覧とその説明を次の表に示します。


Motion correctionジョブのパラメータ
タブ パラメータ名 説明
I/O Input movies STAR file Import/job001/movies.star 動画ファイルの一覧が記録されているSTARファイルをインポートする。BrowseボタンでImport/movies/movies.starを選択すると、エイリアス名moviesがjob001 (真のジョブの名前)に勝手に置き換わった。また、Browseボタンでファイル検索する際、movieのSTARファイルしか表示されないことにも注目。各ジョブタイプに対応したSTARファイルのみが表示される仕様の様。
I/O First frame for corrected sum 1 アラインメント後の積算画像作成に含める最初のフレーム。
I/O Last frame for corrected sum -1 アラインメント後の積算画像作成に含める最後のフレーム。0以下であれば、動画ファイルの一番最後のフレームを指す。Firstに1、Lastに-1を指定したので、つまり動画内のすべてのフレームを積算画像作成に使用することを意味する。
I/O Pixel size (A) 0.885 検出器の1画素が試料面上で何Åに相当するか[Å/pixel]。検出器の画素の大きさ(μmオーダー)を撮影倍率で割れば求まる。
I/O Voltage (kV) 200 加速電圧 [kV]
I/O Dose per frame (e/A2) 1.277 各フレームでの、試料面上での電子線のフルエンス。1フレーム分の露光中に、試料1平方Å当たりに何個の電子が当たったか[electron/A^2]で表す。ドース(dose)とも呼ばれる。
I/O Pre-exposure (e/A2) 0 動画の記録を始める前に試料に当てられた電子線のフルエンス。動画の最初の数フレームを捨てた場合などに使用。(歴史的には、まだ電子直接検出器がなかった時代、液体ヘリウム冷却式のクライオ電顕での撮影時に、電子線誘起性ドリフトを軽減する目的で、撮影前に一定量の電子線を当てることがあり、それをpre-exposureと呼んだ。)
I/O Do dose-weighting? Yes アラインメント後の積算画像作成の際に、累積フルエンスに基づくフレーム重み付け(dose-weighting)を行うか否か。累積フルエンスが多くなるほど試料に損傷が蓄積し、高周波数のシグナルが失われることがわかっているため、累積フルエンスが大きくなる後ろの方のフレームは高周波数成分の寄与を落とすような処理を行う。
I/O Save non-dose weighted as well? No 上でYesを選択した場合、dose-weightingを行なっていないアラインメント後積算画像も保存するか否か。dose-weightingを行わない方がCTF推定がうまくいく場合もある。今回はディスク容量節約のためNoとする。(こうしたことをせずとも十分良質なデータなので、そもそも不要ではある)
Motion Bfactor 150 各フレームに適用するガウス関数型ローパスフィルタのパラメータ[pixel^2]。Motion Correction v1の論文(Xuemin Li et al., Nature Methods, 2013)によれば、Bfactor(B)の式はexp(-Bk^2)で(k:空間周波数)、これがフレームのフーリエ変換に掛け算されて(Bが正の時は)ローパスフィルタとなる。Bfactorが大きいほどローパス効果が強い。高周波数のノイズがアラインメントの邪魔となるため、このような処理をする。Gatan K2-summitで超解像モードを使用した場合はより大きな値を指定した方が良い。
Motion Number of patches X, Y 5, 5 異方的な(場所ごとに異なる)電子線誘起性ドリフトを補正するために、フレームを5x5のパッチに分割し、各パッチ別でドリフト量を推定し、それを多項式函数でフィティングして画素単位でドリフト量を推定、アラインメントする。X, Y共に1とすればフレーム全体のズレのみ推定、アラインメントされる。
Motion Group frames 1 この値が例えば2だと、ドリフトの推定処理前に隣接2フレームが平均化(Grouping)され、その平均化フレームを用いてドリフト推定される。S/Nが低い動画のドリフト推定では有効なオプションだろう。
Motion Binning factor 1 画像のビニング(縮小)率。フーリエ空間の切り出しにより行う。この段階でビニングをするとそれを後々復調するのは難しくなるが、超解像モードで撮影された動画をダウンスケール(ファイルサイズ圧縮?)するのには役に立つ。浮動小数点数を指定できる。ただし、ビニング後の画像のサイズ(画素数)が偶数となる必要がある。超解像モードの場合、Binning factor 2 がよく使われる。
Motion Gain-reference image Movies/gain.mrc カメラのゲイン補正をしない状態で動画を撮影した場合、ゲインリファレンスをここで与えることでゲイン補正できる。
Motion Gain rotation No rotation(0) ゲインリファレンスを回転する必要がある場合はこのオプションを使う。
Motion Gain flip No flipping(0) ゲインリファレンス回転後にその鏡像をとる必要がある場合はこのオプションを使う。
Motion Use RELION's own implementation? Yes MotionCor2のRELION実装版を使用するか否か。RELION実装版はCPUでマルチスレッド動作する。動画フレーム数がスレッド数の整数倍だと効率的に計算される。4k x 4kの動画で6〜12スレッド使えば、GPU実装のMotionCor2と同等の処理速度となる。RELION実装版を使わない場合(Noを選択)は、MotionCor2の実行可能ファイルパス、使用するGPUデバイス番号を指定する必要がある。
Running Number of MPI procs 1 MPIプロセス数。
Running Number of threads 12 今回の動画は24フレームのため、それを割り切れる2, 3, 4, 6, 8, 12, 24スレッドのいずれかを指定するのが効率が良い。ここでは12を指定した。
Running Submit to queue? No ジョブ・キューイングシステムで実行する場合に使用するオプション。以降ではこれに関連する項目は省く。


パラメータ入力後、【Run!】ボタンを押すとジョブが実行されます。ソースからビルドしたRELION使用時、GUI下部のウィンドウに"ERROR: TIFF support was not enabled during compilation"と表示されてジョブが停止した場合はこちらを参照ください → RELIONトラブルシューティング


ジョブの実行にかかる時間は、最近のそれなりのマシンであれば12スレッド実行で約5分程度になります(※本記事執筆者の環境(ページ上部の実行環境をご覧ください)では4.6分でした)。


【Display:】の右側のドロップダウンメニューを押すと、"out: logfile.pdf"または"out: corrected_micrographs.star"を選択できます。"out: logfile.pdf"を選択すると、電子線誘起性ドリフトの推定結果とその統計値を記録したファイルを見ることができます。(※環境によってはpdfファイル生成関係や表示関係でエラーが出て見れないかもしれません。ghost scriptとevinceをインストールすれば大体大丈夫と思います。)


motioncorr logfile button.png


logfile.pdfの内容はこちらで説明します → Motion correctionジョブのlogfile.pdf


また、【Display:】で"out: corrected_micrographs.star"を選択すると、Relion display GUIが立ち上がります。

motioncorr display1.png

パラメータを適当に入力して【Display!】ボタンを押すと、アラインメント後(ドリフト補正後)の積算画像を見ることができます。今回はScale: 0.3(画像縮小)、Sigma contrast: 3 (高コントラスト化)、Nr.columns: 3 (1行あたり3枚の画像を並べてタイル表示)程度が適当でしょう。

motioncorr display2.png


ちなみに各画像上で右クリックして表示されるポップアップメニューで、"Show metadata"を選択すると、corrected_micrographs.starに記述されているその画像のメタデータが表示されます。

motioncorr display3.png


注意点として、このディスプレイを使って画像の選別を行うことはできません。データセットから消したい画像がある場合は、"Subset selection"ジョブタイプを使いましょう。


Motion correctionジョブタイプのstarファイルについてはこちらで説明します → Motion correctionジョブのstarファイル

1.3 CTF estimation (CTF推定)

CTF estimationジョブタイプでは、ドリフト補正後の画像のCTF推定を行います。CTF estimationジョブタイプはKai Zhangが開発したGctfのラッパープログラムになっています。GctfはGPUで動作します。適当なGPUがない場合は、Alexis RohouとNiko Grigorieffが開発したCTFFIND4.1を使うことができます。CTFFINDはCPUで動作します。

以下ではGctfを使った場合について説明します。まずジョブタイプブラウザから"CTF estimation"を選択します。タブはI/O, Searches, CTFFIND-4.1, Gctf, Runningの5種類あり、そのうちGctfで使うのはI/O, Gctf, Runningの3種類です。SearchesタブとCTFFIND-4.1タブは無視されます。以下にパラメータを入力した状態の画面を示します。(Gctfがインストールされていない場合はこちらへ→ Gctf)


gctf-param-1.png

gctf-param-2.png

gctf-param-3.png


パラメータとその説明を以下の表に示します。

CTF estimationジョブのパラメータ
タブ パラメータ名 説明
I/O Input micrographs STAR file MotionCorr/job002/corrected_micrographs.star Motion correctionジョブによるドリフト補正後の画像をリストアップしたSTARファイル。
I/O Use micrograph without dose-weighting? No CTF推定にdose weightingしていない画像を用いるか否か。今回はMotion correctionジョブでdose weightingしていない画像の作成、保存をしていないため、No一択。
I/O Spherical aberration (mm) 1.4 対物レンズの球面収差係数[mm]
I/O Voltage (kV) 200 加速電圧[kV]
I/O Amplitude contrast 0.1 コントラスト全体(位相コントラスト+振幅コントラスト)に対する振幅コントラスト成分の寄与の程度を表す数値(Kai Zhang, J. Struct. Biol (2016)の式(1)を参照)。実際には0.1より小さいことが知られているが、単粒子解析的には0.1程度とした方が良い結果を得られる傾向がある。これは現状のCTFが非弾性散乱由来の低周波数成分をモデル化していないことによる誤差と考えられている。
I/O Magnified pixel size (Angstrom) 0.885 検出器の1画素が試料面上で何Åに相当するか[Å/pixel]。検出器の画素の大きさ(μmオーダー)を撮影倍率で割れば求まる。
I/O Amount of astigmatism (A) 100 想定される対物非点収差の大きさ[Å]。適切に軸調整された電顕ならばこの値で概ね上手くいく。
Gctf Use Gctf instead? Yes CTF推定にGctfを使用するか否か。Gctfを用いないならば代わりにCTFFIND-4.1を使う必要がある。
Gctf Gctf executable Gctfバイナリのファイルパス Gctfバイナリのファイルパスを指定。環境変数$RELION_CTFFIND_EXECUTABLEでファイルパスを指定済みであれば、それがデフォルト値として使われる。参考までに、本記事執筆者環境では/home/kttn/softwares/Gctf/Gctf_v1.18_b2_sm61_cu9.2
Gctf Ignore 'Searches' parameters? Yes Yesの場合、SearchesタブのCTF推定用パラメータは無視され、Gctfのデフォルトパラメータが使用される。Searchesタブでパラメータを指定してGctfに渡したい場合はこれをNoとする。
Gctf Perform equi-phase averaging? Yes 推定したCTFモデルに基づく等位相平均化(equi-phase averaging, EPA)を施したパワースペクトルを生成するか否か。EPAによりトーン(Thon)リングのS/Nが向上し、目視でのCTF推定クオリティの判定に役立つ。あくまでも人間が目で判定するために使うもので、CTF推定計算自体には影響を与えない。Kai Zhang, J. Struct. Biol(2016)の2.9節を参照。
Gctf Other Gctf options Gctfに渡したい追加のパラメータがあればここに記述。
Gctf Which GPUs to use 0 使用するGPUのデバイス番号(複数台のGPUを使用する場合は複数指定可能)。GPUデバイス番号は普通は0から始まり、GPU1台につき1つの番号が割り当てられる。今回は小規模データセットのためGPU1台で十分。
Running Number of MPI procs 1 MPIプロセスを増やしてジョブを並列化可能だが、今回はGPU1台、シングルプロセスでも全ての画像を処理するのに30秒程度しかかからないため、1で十分。


CTF estimationジョブタイプの結果生成されるSTARファイルについてはこちらで説明します → CTF estimationジョブのSTARファイル

ジョブが完了すると、CtfFind/job003/Movie/ディレクトリ以下に、各画像ごとに複数のファイルが生成されます。例えば20170629_00021_frameImage.mrcを例にとると、以下の様な画像が生成されます。(EMAN2のe2display.pyで表示しています。)


20170629_00021_frameImage.ctf : 左側がモデルのパワースペクトル 、右側が画像(実データ)のパワースペクトル 。白いリングは多分_rlnCtfMaxResolutionです。

ctf-ctf.png


20170629_00021_frameImage.epa : 左側がモデルのパワースペクトル 、右側が画像(実データ)のパワースペクトルを等位相平均化(EPA)したもの。EPAによりThonリングがよりよく見えます。

ctf-epa.png


20170629_00021_frameImage.pow : 画像(実データ)のパワースペクトル 。20170629_00021_frameImage.ctfの右半分はこれの右半分と同じもののはず...ですが、.ctfの方は何らかの都合により階調が潰れていて見えづらく、困ります。

ctf-pow.png


また、各画像ごとに.logという拡張子のテキストファイルも生成されます。これはGctfのログファイルになります。


RELION GUIの【Display】ボタンを押し、out: micrographs_ctf.starを選択すると、*.ctfを一覧表示できます(左半分がCTFモデルで右半分が画像のパワースペクトル のやつ)。Thonリング(白いリング)の間の暗環(ゼロ点と呼びます)とCTFモデルの暗環は一致しているべきです。RELION display GUIのSort imagesにチェックを入れて、on: で好きな指標を選択すると、その指標に基づいて並び替えられた状態で表示されます。

以下は例です。_rlnCtfMaxResolutionの大きいものから順に並べています(混乱しそうになりますが、大きい方が要はフィッティングの質が悪いです)。フィッティングの質に基づいてスクリーニングしたいときなどに便利です。(このチュートリアルのデータではなぜかパワースペクトルが潰れて見えないので難しいですが。)

ctf-disp1.png

ctf-disp2.png


また、【Display】ボタンからlogfile.pdfを選択すると、micrographs_ctf.starに記録されている各種推定値のプロットとヒストグラムを確認できます。推定値の説明はこちらを参照ください→ CTF estimationジョブのSTARファイル。撮影時に狙った通りの分布になっているかなどを確認しましょう。


CTFモデルの推定(フィッティング)が悪いと思われる画像がある場合は、それらの画像のみパラメータを変えて再推定することができます。まずCtfFind/job003/Movies/の中から、再推定したい画像の.logファイルを削除します。次に、GUIの【Finished jobs】リストの中からCtfFind/job003を選択し、パラメータを変更し、【Continue now】ボタンを押下すれば新規パラメータで再推定が実行されます(【Run!】ボタンが【Continue now】ボタンに変わっています)。CTF推定結果に満足するまでこのプロセスを繰り返すことができます。もしパラメータを変更して再推定してもうまくいかないものがあれば、【Subset selection】ジョブタイプでそれらの画像を以降の解析から除外しましょう。

1.4 Manual particle picking (手作業による単粒子ピッキング)

次のManual pickingジョブタイプは、(ドリフト補正により平均化された)画像から手作業で単粒子の座標を指定する(ピッキングする)ために使います。自分のデータの素性をよく知るためにも、少なくとも数枚の画像で手作業単粒子ピッキングを行うのが良いでしょう。このジョブタイプの典型的な使い方は、ここで手作業で指定した単粒子を用いて参照像無しで2Dクラス平均画像を計算し、それを参照像としてデータセット全体に対し自動単粒子ピッキングを行うというものです。しかしRELION 3.0から、ラプラシアン-オブ-ガウシアン(LoG)フィルターによる参照像無し自動単粒子ピッキングが実装されています。これまでのところ、ほとんどの場合において、LoGによる自動ピッキングでまずまずの初期座標を得られており、それゆえManual pickingジョブタイプはスキップできます。RELION付属の全自動単粒子解析スクリプトであるrelion_it.pyもLoGによる自動ピッキングを使用しています。今回このチュートリアルでは、あくまでも説明目的でManual pickingジョブタイプを触ってみることにし、その後LoGベースの自動ピッキングで単粒子の初期座標を得ることにします。


単粒子の手動ピッキングは経験あるのみです!もしRELIONのピッキングインターフェイスが使いにくければ、次のプログラムを代わりに使っても良いです(RELIONはそれらが出力する座標ファイルのフォーマットにも対応しています):Jude Shortが開発したXIMDISP (座標ファイルの拡張子は任意)、XMIPP-2.4(座標ファイルの拡張子は任意)、Steven Ludtkeが開発したe2boxer.py(座標ファイルの拡張子は.box)。それらを使用する場合には、以下のことを守ってください。

  • 座標ファイルはテキストファイルとして、画像または動画のインポート元ディレクトリ(最初のImportジョブタイプで動画を読み込む時に指定したディレクトリ)に保存すること。
    • 今回の場合はMovies/ディレクトリです。
  • 座標ファイルのrootname(ファイル名から拡張子を除いた文字列)は対応する画像ファイルのrootnameと同じにし、拡張子のみ変更すること。
    • 例として、20170629_00021_frameImage.tiff の座標ファイルの場合、Movies/20170629_00021_frameImage.box, 20170629_00021_frameImage.star など。

ピッキングができたら、Importジョブタイプより、Node type: "2D/3D particle coordinates"を選択し、Input filesに該当する座標ファイルをワイルドカードで指定(Movies/*.box、Movies/*.starなど)して座標をインポートしてください。


ここではRELIONのピッキングインターフェイスについて説明します。(e2boxer.pyによる方法はこちら)


まずManual pickingジョブタイプのI/Oタブにて、【Browse】ボタンを押し、CtfFind/job003/以下に生成されたmicrographs_ctf.starファイルを選択します。Colorsタブはここでは無視します。Displayタブには以下の通り入力します。

manualpick-display.png

Manual pickingジョブのパラメータ
タブ パラメータ名 説明
Display Particle diameter (A) 200 単粒子(蛋白質分子)の大きさ[Å]。ピッキング時に表示される円の大きさを指定するパラメータで、解析には影響を与えない。みやすい値を選べば良い。
Display Scale for micrographs 0.25 画像の表示スケール。使用しているモニタの大きさ、解像度に合わせて適宜調整する。
Display Sigma contrast 3 典型的な電顕画像は"シグマ(標準偏差)コントラスト"により見やすくなる。3とすると、画像の輝度値の平均 - 3*標準偏差が黒、平均 + 3*標準偏差が白となる様、かつ線形にヒストグラム伸長して画像を表示する。詳細はRELION wikiのDisplayImagesエントリを参照。
Display White value 0 ヒストグラム伸長の際に白となる値を自分で指定する場合はここに入力する。その場合、Sigma contrastは0とすること。
Display Black value 0 ヒストグラム伸長の際に黒となる値を自分で指定する場合はここに入力する。その場合、Sigma contrastは0とすること。
Display Lowpass filter (A) 20 画像にローパスフィルタを施してから表示する。カットオフをÅで指定する。ノイズの大きい画像で単粒子が見やすくなる場合がある。
Display Highpass filter (A) -1 画像にハイパスフィルタを施してから表示する。カットオフをÅで指定する。画像全体に乗っている明暗の勾配を除去するのに有用。-1はハイパスフィルタ無効化を意味する。
Display Pixel size (A) 0.885 ピクセルサイズ[Å/pixel]。周波数フィルタのカットオフと単粒子の大きさをpixel単位に換算するために使用。(画像のピクセルサイズと検出器のピクセルサイズのどちら??)
Display Pick start-end coordinates helices? No 線維状蛋白質をピッキングする場合に指定するオプション。
Display Scale for CTF image 1 画像の【CTF】ボタンを押下した際に表示されるパワースペクトル画像の表示スケールを指定。


【Run!】ボタンを押すと、RELION manual-picking GUIウィンドウが立ち上がります。画像名、ピッキング済みの単粒子数、デフォーカス量がリストアップされています。これを利用して単粒子ピッキングするわけですが、今回のチュートリアルではLoGによる自動ピッキングを使って解析を進めるため、ピッキングはしなくても良いです(ただしこのジョブを実行(Run)したことにより生成されるstarファイル自体は使用するので、ここまではやってください)。

pick1.png


以降の解析では使いませんが、ピッキングの方法を説明しておきます。【pick】を押せば画像が表示されます。左クリックで単粒子ピッキング、中クリックでピッキングの取り消しができます。また、右クリックでメニューが表示されます。ピッキングを終えたら右クリックメニューの"Save STAR with coordinates"を押し、座標ファイルを保存する必要があります。座標ファイルの保存は画像ごとに個別に行う必要があります。

pick2.png


座標ファイルを保存してさえいれば、ピッキングは一時中断して後からまたやり直すことができます。その場合、RELION GUIのFinished jobsペインからManualPick/job004を選択し、【Continue now】を押せばいいです。

手作業ピッキングした単粒子を使って2D classificationジョブタイプでまずまずのクラス平均像を得るためには、およそ500-1000個の単粒子が必要です。

1.5 LoG-based auto-picking (LoGに基づく自動ピッキング)

ここではラプラシアン-オブ-ガウシアン(LoG)フィルタを用いて、参照像(reference)無しでの自動ピッキングを行います。ここでピッキングした単粒子を次の2D classificationジョブタイプに渡し、Auto-pickingジョブタイプで用いる参照像を作成するというのが後の流れです。この段階ではそれほど多くの単粒子は必要無いので、今回は最初の5枚の画像についてのみLoGベースのピッキングを行うことにします。ただし、これはチュートリアルの時間を短縮するためです。実際のプロジェクトでは、より良い参照像を得るために、使える画像全部でLoGピッキングを行なっても良いでしょう。(LoGピッキングについてはZivanov et al., eLife(2018)を参照。)

まずはLoGピッキングを行う画像を選択するため、"Subset selection"ジョブタイプを使用します。I/Oタブの"OR select from picked coords"で、先ほどの"ManualPick/job004/coords_suffix_manualpick.star"を指定し、Current job:にエイリアス名 "5mics" と入力し、【Run!】を押します。 (※ Manual pickingジョブタイプでRun!するだけで生成されるSTARファイルです。単粒子を手動ピッキングしていなくても生成されています。)

subset1.png


Manual pickingジョブタイプと同様のRELION manual-picking GUIが表示されます。デフォルトでは全ての画像のチェックボックスにチェックが入っており、すなわち全ての画像が選択された状態です。manual-picking GUIのFileメニューでInvert selectionを押すと、全てのチェックボックスが外れて画像が選択されていない状態になります。最初の5枚の画像だけチェックを入れて、FileメニューでSave selectionを押します。

subset2.png


すると、Select/5mics/micrographs_selected.starに、選択した画像のファイル名やCTF情報が保存されます。これを次のLoGピッキングに用います。


LoGピッキングを行います。Auto-pickingジョブタイプを選択し、以下の様にパラメータを入力します。また、エイリアス名としてLoG_basedを指定しておきましょう。Referencesタブ、helixタブはここでは無視します。

log1.png

log2.png

log3.png

LoG pickingのパラメータ
タブ パラメータ名 説明
I/O Input micrographs for autopick Select/job005/micrographs selected.star 選択した5枚の画像情報(CTF情報含む)を記録したSTARファイル。
I/O Pixel size in micrographs (A) -1 上記でCTF情報を含むSTARファイルを指定した場合、ここで負の値を与えれば、STARファイルに記録されているピクセルサイズ[Å/pixel]が使われる。
I/O 2D reference LoGピッキングは参照像(reference)無しなので空欄にする。
I/O OR: provide a 3D reference? 同上
I/O OR: use Laplacian-of-Gaussian Yes LoGピッキングを実行する。
Laplacian Min. diameter for LoG filter (A) 150 単粒子の最小サイズ [Å]を指定すれば良い。e2display.pyなどを使って画像からおおよその値を見積もる。
Laplacian Max. diameter for LoG filter (A) 180 単粒子の最大サイズ [Å]を指定すれば良い。e2display.pyなどを使って画像からおおよその値を見積もる。
Laplacian Are the particles white? No 単粒子が背景よりも白い場合(コントラストを反転したり、CTFの位相フリッピングを行なった場合など)はYesとする。今は単粒子が背景よりも黒いのでNoとする。
Laplacian Maximum resolution to consider (A) 20 ここで指定した分解能が画像内の最大の分解能となる様に画像を縮小した上でLoGピッキングが実行される。今回はデフォルトの20Åのままで良い。
Laplacian Adjust default threshold 0 正の値にすると、すなわち閾値を高くすると、ピッキング数が減る。負の値にすると、すなわち閾値を低くするとピッキング数が増える。良い値はおそらく[-1, 1]の範囲にあるが、多くの場合デフォルトの閾値0で十分良い仕事がされる。
Reference LoGピッキングでは無視。
autopicking Picking threshold LoGピッキングでは無視。
autopicking Minimum inter-particle distance (A) LoGピッキングでは無視。
autopicking Maximum stddev noise LoGピッキングでは無視。
autopicking Minumum avg noise LoGピッキングでは無視。
autopicking Write FOM maps? No 参照像ベースのピッキングで使用する。
autopicking Read FOM maps? No 参照像ベースのピッキングで使用する。
autopicking Shrink factor 0 (RELION 3.0.5の場合)LoGピッキングではこの欄は無視され、常に0が使われる。0の場合、Laplacianタブで指定したMaximum resolution to considerの値に従い画像が縮小される。
autopicking Use GPU acceleration? No LoGピッキングは十分高速なのでGPU化されていない。
Helix LoGピッキングでは無視。


【Run!】を押せばLoGピッキングが実行されます。シングルMPIプロセスでおよそ40秒程度です。

【Display:】ボタンを押してcoords_suffix_autopickを選択するとRELION manual-picking GUIが立ち上がり、ピッキングの結果を確認できます。

rln-log4.png


Manual pickingと同じ要領で、追加で単粒子のピッキングや削除をすることができます。また、Laplacianタブの各パラメータを調整して再計算することもできます。しかし参照像あり自動ピッキングの参照像を作る目的であればこの程度のクオリティで多くの場合問題ないため、このまま先へ進みます。


また、【Display:】でout:logfile.pdfを選択すると、各画像でピッキングされた単粒子数(_rlnGroupNrParticles)のプロットとヒストグラム、各画像のピッキングの平均FOM(ピッキングスコアの指標, _rlnAutopickFigureOfMerit)のプロットとヒストグラムを確認できます。


LoGピッキングの結果生成されるSTARファイルについてはこちらをご覧ください→ LoGピッキングのSTARファイル

1.6 Particle extraction(単粒子の切り出し)

単粒子を切り出したい画像の座標ファイルを得られたら、Particle extractionジョブタイプにより単粒子を切り出し、そのメタデータを集めることができます。Particle extractionジョブタイプを選択し、以下の様にパラメータを指定します。あとから行う切り出しと区別するため、LoG_basedとうエイリアス名を指定しておきます。

rln-extract-log-1.png

rln-extract-log-2.png

LoGピッキング結果に基づくParticle extractionのパラメータ
タブ パラメータ名 説明
I/O micrograph STAR file CtfFind/job003/micrographs_ctf.star Browseボタンを使って指定する。Select/5mics/にあるmicrographs_selected.starを使うこともできる。いずれにせよ単粒子の座標は最初の5枚の画像についてしか用意していないため、どちらを使っても良い。単粒子の座標が指定されていない画像があると【Run!】した時にGUIの下部ウィンドウ(stderr)に赤字で警告が出るが、無視して良い。(※他のエラーや警告まで見逃さない様に注意)
I/O Input coordinates AutoPick/job006/coords_suffix_autopick.star Browseボタンを使って指定する。
I/O OR re-extract refined particles? No このオプションは2D classification、3D classification、3D auto-refineジョブタイプで生成される _data.star ファイルを用いて、STARファイル内に記録されている情報を元に単粒子を再切り出しするために使う。例えば、はじめに単粒子を縮小して切り出し、2Dまたは3Dの分類(classification)を行った後に、今度は元の縮小前のスケールで切り出し直したい、といった場合に有用。RELION3ではこの機能をさらに拡張した"re-center refined coordinates"というオプションが追加されており、3D classificationや3D auto-refineで用いる3D参照像内における任意の点(x, y, z)を指定できる。それを用いると、全ての単粒子のx, y座標の原点が調整され、その結果得られる3D再構成像はそこで指定したx, y, z座標が中心にくる。Focused refinementを行う際に有用。
I/O Manually set pixel size? No 入力STARファイル(micrograph STAR file)にピクセルサイズの情報が記録されていない場合はこれをYesにし、Pixel size[Å/pixel]を手入力する。
extract Particle box size (pix) 256 単粒子の切り出しサイズ [pixel]。256x256の大きさで切り出される。偶数でなければならない。
extract Invert contrast? Yes 白黒反転して単粒子が背景よりも白くなる様にする。
extract Normalize particles? Yes 画像は常に正規化する。正規化については下記参照。
extract Diameter background circle (pix) 200 単粒子の中心を中心とする直径200pixelの円の外側を背景領域と見なし、背景領域の輝度値の平均が0、標準偏差が1となる様に単粒子画像を正規化する。負の値を指定すると、Particle box sizeの75%の値が使われる。
extract Stddev for white dust removal -1 異常に明るい画素(があれば)を取り除くために使用。平均よりも標準偏差xこの値倍を超えて大きい輝度値を持つ画素は、ガウス分布からサンプリングした別の値で置き換えられる。5程度が典型的。今回のデータセットは異常画素が存在しないので、-1を指定して無効化する。
extract Stddev for black dust removal -1 異常に暗い画素(があれば)を取り除くために使用。平均よりも標準偏差xこの値倍を超えて小さい輝度値を持つ画素は、ガウス分布からサンプリングした別の値で置き換えられる。5程度が典型的。今回のデータセットは異常画素が存在しないので、-1を指定して無効化する。
extract Rescale particles? Yes 単粒子画像を縮小するか否か。縮小すると計算が高速化される。最初のクラス分類を高速化するために、解析の初期段階では単粒子画像を縮小するのが典型的。再構成像の分解能がナイキスト周波数に近づいてきたら、縮小なしで単粒子を再切り出しし、高分解能解析に進む。
extract Re-scaled size (pixels) 64 縮小後の単粒子画像のサイズ[pixel]。今回は256x256 → 64x64へと縮小する。


Helixタブは今回は無視し、RunningでシングルMPIプロセスにして【Run!】実行します。


単粒子画像はMRCスタックとして切り出し元画像ごとに1つのファイルにまとめて書き出されます(RELIONでは.mrcsの拡張子を持ちます)。保存場所はExtract/job007/Movies/になります。切り出し後は、特に問題が起きていないか確認するため、単粒子画像の様子をざっと確認しましょう。【Display】でout:particles.starを選択するとRelion display GUIが立ち上がるので、適当にパラメータを指定して【Display!】を押し、表示します。

rln-extract-log2.png


切り出した単粒子画像の平均画像を確認することもできます。画像ウィンドウ上で右クリックし、メニューから"Invert selection"を選択すると全ての単粒子画像が選択され(選択されると赤枠で表示されます)、さらに右クリックメニューから"Show average of selection"を選ぶと平均画像が表示されます。

rln-extract-log3.png


Particle extractionジョブで生成されるSTARファイルについてはこちらで説明します → Particle extractionジョブのSTARファイル

1.7 Making templates for auto-picking (自動ピッキングのための参照像作成)

参照像ベースの自動ピッキングで用いる参照像を得るために、2D classificationジョブタイプを利用します。

2D classificationジョブタイプを選択し、以下の通りパラメータを入力します。

rln-2dclass-log1.png

rln-2dclass-log2.png

rln-2dclass-log3.png

rln-2dclass-log4.png

rln-2dclass-log5.png


【Run!】を推した時にMPI系のエラーが発生した場合 → RELIONトラブルシューティング


自動ピッキングのための参照像作成パラメータ
タブ パラメータ名 説明
I/O Input images STAR file Extract/job007/particles.star 2D classificationにかけたい全ての単粒子画像とそのメタデータを記録したSTARファイル。CTF補正を行う場合はCTF情報も記録されていることが必要。
CTF Do CTF-correction? Yes ベイズ推定(MAP最適化)により最適なウィーナーフィルタが推定され、CTFの振幅および位相補正を実行する。
CTF Have data been phase-flipped? No 画像がCTFの位相反転済みであるか否か。RELION以外で前処理した画像を扱う場合のみ有用なオプション。
CTF Ignore CTFs until first peak? No Yesとした場合、CTFの振幅補正はCTFの第1ピーク以降の周波数帯でのみ行われる。CTF補正した結果低分解能が異常に強くなり、単粒子がごく少数のぼやけたクラスに集まってしまった場合などは、これをYesとすることが有用な場合あり。
Optimisation Number of classes 50 クライオ電顕の場合は少なくとも1クラスあたり100粒子あることが望ましい。負染色の場合はより少数、例えば1クラスあたり20-50粒子程度が目安。しかし、今回のように単粒子数が少ない場合(~1,500粒子)、この制限を緩和した方が単粒子が異なるクラスにうまく分散する傾向がある。単粒子数に関わらず、少なくとも50クラスに分けることは悪い考えではない。(※1クラスあたり100粒子にこだわろうとすると今回の場合15クラス程度にせねばならないが、それだとうまくクラス分類できない場合があるので、条件を緩和して50クラス程度にした方がいい、ということ)
Optimisation Regularisation parameter T 2 Tの定義は Sjors H. W. Scheres, JMB(2012)を参照。T=1はベイズ理論的に導出される最適化用の数式をそのまま用いることに相当する。T>1とすると、実験データ(画像)由来の情報により重みを置いて最適化計算を行う。クライオ電顕画像の2D classificationであればT=2〜3, 3D classificationであればT=3〜4が典型的な値である。負染色の場合はより小さめのTを使うのが良い。クラス平均像のS/Nが悪い場合には、Tを下げた方が良い。クラス平均像の分解能が低すぎる場合は、Tを上げると良い。高周波数のノイズに対しオーバーフィッティングしてしまうのを避けることが重要。
Optimisation Number of iterations 25 最適化は逐次計算により行われる。このパラメータで逐次計算の反復(イテレーション)回数を指定する。25以外の値を使うことは滅多にない。
Optimisation Use fast subsets (for large data sets)? No Yesとすると、最初の5イテレーションはランダムに選択された(クラス数×100)粒子のみで計算を実行し、次の5イテレーションは(クラス数×300)粒子で計算し、次の5イテレーションは全単粒子の30%で計算し、残りは全データで計算する。これはNiko Grigorieffらが開発したcisTEMの実装から着想を得たものである。数十万粒子あるような大きなデータセットで、計算の高速化のために大変有用な仕組みである。今回のような小さなデータセットでは不要。
Optimisation Mask diameter (A) 200 単粒子およびクラス平均像に対し適用するソフトな円形マスクの直径[Å](※ソフトなマスクとは、エッジがステップ状ではなく滑らかなマスクのこと)。溶媒領域および隣接単粒子由来のノイズを除去する役割。ノイズを出来るだけ除くためにマスクを小さくしつつも、単粒子およびクラスの情報をマスクで消してしまわないよう、単粒子の最大サイズよりも大きく取らねばならない。
Optimisation Mask individual particles with zeros? Yes Yesとした場合、個々の単粒子のマスクの外側部分が、フーリエ変換前に値ゼロで埋められる。これによりノイズが除去され、アラインメントとクラス分類の感度が向上する。一方で値ゼロで埋まったマスクの存在により偽相関が生じる。Noとした場合、マスクの外側はランダムなノイズで埋められ、したがってマスクによる偽相関は生じない。大きな分子の高分解能リファインメントでは溶媒をランダムノイズで埋めた方がうまくいく傾向があるが、それ以外のほとんどのケースでは値ゼロで埋めた方がうまくいく傾向がある。
Optimisation Limit resolution E-step to (A) -1 正の値を指定すると、EMアルゴリズムのEステップ(アラインメント計算)は、ここで指定した分解能[Å]までの情報しか使わなくなる。ノイズへのオーバーフィッティングを防ぐために有用な場合がある。今回は特にこの機能を使わないが、もし指定するとすれば10

〜15Å程度だろう。非常にノイズの多いデータなど、クラス分類が難しい場合はこれが役立つことがある。0以下の値を入れると無効化。

Sampling Perform image alignment? Yes 画像のアラインメント(回転、並進移動)を行うか否か。Noとするとクラス分類のみを実行する。極めて細かいfocused classificationを行う場合などにNoにしたりするが、うまくいくためにはSTARファイル内の各単粒子のオリエンテーションが十分最適化されていることが前提である。
Sampling In-plane angular sampling 6 面内回転角(2D画像の回転)のサンプリングレート[度]。細かくとると計算に時間がかかる。ほとんどの2Dリファインメントにおいては6度で十分である。自動サンプリングを使う場合には、この値は最初のイテレーションのみで使われ、それ以降は自動でサンプリングレートが上がっていく。
Sampling Offset search range (pix) 5 並進移動に関する確率計算は、この値[pixel]を直径とする円の内側だけで行われる。円の中心はイテレーションごと、喀痰粒子ごとに更新され、直前のイテレーションにおいて最適な並進移動位置に置かれる。自動サンプリングを使う場合には、この値は最初のイテレーションのみで使われ、それ以降は自動で調整される。
Sampling Offset search step (pix) 1 並進移動のサンプリングレート[pixel]。
Helix 螺旋型繊維状蛋白質複合体の解析用。今回は無視。
Compute Use parallel disc I/O? Yes Yesとすると、全てのMPIスレーブプロセスがそれぞれ自分用の単粒子を読み込む。高速なファイルシステムを使っている場合はこれをYesにすると良い。ファイルシステムが並列化に対応していない場合、複数のMPIノードから並列アクセスすることは不可能であり、その場合はNoとする。Noの場合、MPIマスタープロセスのみが単粒子の読み込みを行い、ネットワーク通信によりMPIスレーブへ単粒子を分け与える形となる。
Compute Number of pooled particles 30 MPIスレーブの中ではI/O高速化のため、単粒子画像はバッチ(※複数枚のまとまり)単位でやりとりされる。最低でも1スレッドあたり1粒子以上となるバッチサイズが必要である。Number of pooled particlesは1スレッドあたりに一度に渡される単粒子画像の枚数を指定する。1MPIプロセスあたり8スレッドの時、この値が30であれば、MPIスレーブのバッチサイズは 30×8 = 240 枚になる。ディスクアクセスが遅いシステムを使っているときは、これを大きめにすると効率が上がり高速化される。典型的には、GPUを使用する場合は10〜30枚程度、CPUのみを使う場合はもっと小さな値、例えば3などで良い。大きな値を使うとその分RAM使用量が増えることに留意。
Compute Pre-read all particles into RAM? Yes Yesとすると、全ての単粒子画像を一度にRAMに読み込む。ディスクアクセスが遅いシステムでは大変有用である。ただしRAM使用量に注意せねばならない。単粒子は倍精度浮動小数点数で読み込まれるため、Nを単粒子画像枚数、box_sizeを画像サイズ[pixel]とすると、(N × box_size × box_size × 4) / (1024 × 1024 × 1024) ギガバイトとなる (※倍精度なら4バイトではなく8バイトでは...?あるいは単精度の間違い?) 。parallel disc I/OがYesの場合、全てのMPIスレーブがそれぞれ全ての単粒子を読み込む。parallel disc I/OがNoの場合、MPIマスターのみが全ての単粒子を読み込み、それをネットワーク通信によりMPIスレーブへ分け与える形になる。
Compute Copy particles to scratch directory Pre-read all particles into RAMがNoの時のみ指定可能。全ての単粒子を一度に読み込めるほどのRAM容量はないが、SSDのような高速なストレージをマウントしている場合、一旦全ての単粒子をSSDへ一時的に移し(移し先のディレクトリをscratch directoryと呼ぶ)、そこから各MPIプロセスがデータ読み出しを行うようになる。
Compute Combine iterations through disc? No Yesの場合、各イテレーションごとに、各MPIスレーブが各々の計算結果ファイルをディスクへ書き出す。MPIマスターはそれらを読み込み、1つの大きなファイルへと結合し、ディスクへ書き出す。各MPIスレーブは結合された結果をそれぞれ読み込み、次のイテレーションの処理へ進む。こうするとネットワーク通信の負荷は軽減されるが、ディスクI/Oの負荷が大きくなる。Noとした場合、マスターとスレーブ間のデータやり取りはディスクを介さず、ネットワーク通信により実行される。ディスクI/Oとネットワーク通信のどちらが強いかによって、各実行環境ごとで適宜変更すること。
Compute Use GPU acceleration? Yes 計算をGPUで実行する。
Compute Which GPUs to use 0:1:2:3 このチュートリアルを実行する環境ごとで適宜変更すること。空欄にすればシステムが自動的に決めて実行される。ここで示すようにデバイスIDを手動で指定することも可能。複数のMPIプロセスを使う場合、各MPIプロセスごとに使用するGPUデバイスIDを変えることおが可能。各MPIプロセスに対し割り当てるGPUデバイスIDはコロンで区切って指定する。例えばMPIプロセス4つ(プロセス0〜プロセス3)で計算する場合に 0:1:0:1 と指定すると、プロセス0および2はGPUデバイス0, プロセス1および3はGPUデバイス1を使って計算する。
Running Number of MPI procs 5 MPIの全プロセス数。使用されるCPUコア数は、この値と次のNumber of threadsとの積になる。2D classification, 3D classification, 3D initial model, 3D auto-refineジョブタイプでは、MPIプロセスのうち1つはマスタープロセスとなり、それ自体は計算を行わずに他のスレーブプロセスへのジョブ送信を担当する。したがって、上記のように4GPUで計算したい場合には、MPIプロセスを5個立ち上げれば良い。
Running Number of threads 6 MPIプロセスあたりのスレッド数。スレッドの方がメモリ効率が良く、MPIプロセスの方がスケーリング性能が高い(※数を増やした時の並列化効率がスレッドよりも高い)。3D classificationや3D auto-refineはRAM消費が激しいため、スレッドを増やしてRAM使用の効率化を行うと良いだろう。2D classificationはそれほどRAM消費が高くないことから、スレッドを多くとる必要はあまり無い。


以上のパラメータに関し、注意点として、まずSamplingタブの内容はデフォルト値から変える必要は滅多にありません。今回もデフォルト値を使用しています。ほとんどのプロジェクトでは6度の面内回転サンプリングで十分です。一部、巨大な正20面体ウイルス粒子では、より細かいサンプリングにすると結果が良くなります。その場合、まずはデフォルトの6度サンプリングで25イテレーションを実行し、その後【Continue now】ボタンを使って、例えば2度のサンプリングにて追加で5イテレーション程度実行すると良いです(Number of iterations: 30とすれば追加で5回回せますNumber of iterationsはこれから走らせるイテレーション数ではなく、そのジョブの累積イテレーションを表します)。今回のデータセットに関しては、これは全く必要ないです。

また、【Continue now】ボタンは、何らかの理由により失敗したり停止したジョブを再開するのにも使えます。その場合、ジョブのパラメータは変えてはなりません。2D classification, 3D classificationおよび3D auto-refineジョブを再開するためには、続きを計算したいイテレーションで生成された *_optimizer.star ファイルをI/Oタブで指定することが必要です。


ComputeタブおよびRunningタブのパラメータは実行環境ごとに適切な設定をしてください。


2D classificationジョブはこれ以外にも実行するため、区別のためにLoG_basedというエイリアス名をつけておきましょう。【Run!】を押して計算が終了したら、【Display】でout: run_it025_model.starを選択し、クラス平均像を確認しましょう。

rln-class2d-log-result.png


Relion display GUIで Sort images on "rlnClassDistribution"としReverse sortにチェックを入れるか、またはSort images on "rlnAccuracyRotations"としReverse sortのチェックを外すかして【Display!】すれば、このように良いクラス平均像が左上に集まる傾向があり、結果の質を判断しやすくなります。


ちなみに、クラス平均像を左クリックすると赤枠でそのクラスが選択され、右クリックして"Show particles from selected classes"を選択すると、そのクラスに属する単粒子画像が表示されます。

rln-class2d-log-result2.png


また、本ジョブの結果生成されるスタックファイルについて、Class2D/LoG_based/ 以下に、各イテレーションごとのファイル (run_it<イテレーション番号>_* と、run_unmasked_classes.mrcsというファイルが生成されています。

run_it<イテレーション番号>_classes.mrcs が、各イテレーションにおける2D classificationの結果のクラス平均像を収めたスタックです。run_unmasked_classes.mrcsは最終イテレーション(25)のクラス平均像の、マスク無し版です(それでも何らかの大きなマスクがかかっているようには見えますが)。

2D classificationジョブで生成されるSTARファイルの内容ついてはこちらに書いていきます → 2D classificationジョブのSTARファイル

1.8 Selecting templates for auto-picking (自動ピッキングに使う参照像を選ぶ)

自動ピッキングの参照像として使う2Dクラス平均像の選択は、Subset selectionジョブタイプで行います。

Subset selectionジョブタイプを選択し、I/Oタブの"Select classed from model.star"に"Class2D/LoG_based/run_it025_model.star"を指定します。Subset selectionはこの後も使用するため、区別のためにtemplates4autopickというエイリアス名をつけておきましょう。

ここで他に注目すべきパラメータは次の2つです。

自動ピッキング用参照像選択時に注目すべきオプション
タブ パラメータ名 説明
Class options Re-center the class averages? Yes Yesとすると、2Dクラス平均像の自動センタリングが実行される。平均像がその重心が画像の中心になるようセンタリングされる。単粒子が背景に対し白くなっていることが必要。2D classificationを走らせるとクラス平均像の重心が動いてしまう傾向があるため、このセンタリング処理は多くの場合必要。特に自動ピッキングの参照像としてクラス平均像を使う場合はセンタリングは特に重要になる。参照像の中心がずれていると、それでピッキングされる単粒子の中心もずれてしまい、データセットに系統的な位置ずれが生じてしまうため。
Class options Regroup the particles? No このオプションは、Subset selectionにより画像1枚あたりの単粒子数が非常に少なくなってしまう場合に有用。画像1枚あたりの単粒子数が非常に少ない場合はノイズのパワースペクトル およびスケール因子(scale factors. ?)の推定が安定しなくなる。後者はデフォルトでは切り出し元画像それぞれで独立に計算される。このオプションでYesとすると、いくつかの切り出し元画像をまとめて一つのグループにし、グループあたりの単粒子数を増やせる。ClassificationジョブまたはAuto-refineジョブを実行する際、グループあたりの単粒子数が非常に少ない時はRELIONがその旨警告を出してくれるので、その時対処すれば良いだろう。


【Run!】を押し、RELION display GUIを立ち上げます。

RELION display GUIで適宜パラメータを調整してクラス平均像の一覧を表示します。このときrlnClassDistributionでreverse sortすると良いクラスが左上に集まる傾向があるため、参照像を選ぶのが楽になります。

rln-refselect-2.png


左クリックするとクラス平均像が選択され、赤枠で表示されます。クラス平均像をいくつか選択し終えたら、右クリックで表示されるメニューから"Save selected classes"を選択し、選択結果を保存します。


クラス平均像選択のポイントを以下に挙げます。

  • 色々な見た目のクラス平均像があるかと思いますが、それぞれの見た目の中で代表的と思われるクラス平均像をいくつか選択します。 (※見た目viewとは、面内回転psiの違いを無視した時の単粒子の見た目を意味します。)
  • 似たような見た目のクラス平均像を複数選択するのは避けましょう。
  • 悪いクラス平均像(ぼやけていたり、タンパク質らしくなかったり、ゴミのようにごちゃごちゃしたものなど)を選択するのは避けましょう。

ここでは4つのクラス平均像を選択しました。

rln-autopick-template.png


このジョブで生成されるSTARファイルに関しては特筆すべきところはありません。選択したクラス平均像のメタデータが class_averages.starに、選択したクラスに属する単粒子のメタデータが particles.starに保存されます。また、backup_selection.starには、50クラスあったうちのどのクラスが選択されたかがブール値のフラグ列で記録されています。

1.9 Auto-picking(自動ピッキング)

それでは、選択した2Dクラス平均像を参照像として、Auto-pickingジョブタイプにより参照像ベースの自動ピッキングに進みます。

しかし、全画像で自動ピッキングを行う前に、自動ピッキングのパラメータチューニングをしましょう。チューニング対象はautopickingタブのPicking threshold, Minimum inter-particle distance, Maximum stddev noise, Minimum avg noiseです。時間を節約するため、少数の画像のみを使用します。LoGピッキングを行なったのと同じ5枚の画像を使うことにします。

まずは以下のようにパラメータを入力してください。(※job番号が都合により異なっています)

rln-autopick-tune-1.png

rln-autopick-tune2.png

rln-autopick-tune3.png


以下にパラメータをまとめます。

自動ピッキングのパラメータチューニングのためのパラメータ
タブ パラメータ名 説明
I/O Input micrographs for autopick Select/5mics/micrographs_selected.star チューニング用に単粒子をピッキングしたい画像をリストアップしたSTARファイルを指定する。
I/O Pixel size in micrographs (A) -1 -1とするとピクセルサイズ[Å/pixel]は入力STARファイルから自動で読み込まれる。
I/O 2D references Select/templates4autopick/class_averages.star 自動ピッキングに使う参照像のSTARファイルを指定する。
I/O OR: provide a 3D reference? No 3Dの参照像を用いる時に使うオプション。
I/O OR: use Laplacian-of-Gaussian? No 今回はLoGピッキングではないのでNo。
Laplacian 今回は無視する。
References Lowpass filter references (A) 20 参照像にローパスフィルタをかけてからピッキングする。最終的に目標としている分解能よりも十分に低いカットオフでローパスをかけることが必要。これはいわゆる”Einstein-from-noise"を避けるため。高分解能の参照像を使ってしまうと高周波数ノイズに対するオーバーフィッティングが生じてしまう。
References Highpass filter (A) -1 正の値(例えば200Å)を指定すると、ピッキング前に画像に対しハイパスフィルターが施される。画像内に強いグレースケールの勾配がある場合に有効。
References Pixel size in references (A) 3.54 参照像のピクセルサイズ[Å/pixel]。負の値を指定すると、参照像のピクセルサイズは入力画像と同じ(0.885Å/pixel)と見なされる。しかし今回は参照像作成のための2D classificationにおいて、画像サイズを 256x256 pixel から 64x64 pixelへダウンサンプリングしてあるため、参照像のピクセルサイズは 0.885 x (256 / 64) = 3.54 [Å/pixel]である。これを手打ちする必要がある。
References Mask diameter (A) -1 負の値を与えると、マスクの直径[Å]は参照像作成時の2D classificationと同じ値が使われる。
References In-plane angular sampling (deg) 5 面内回転の全探索時の角度刻み[度]。ほとんどの場合5度で十分です。
References References have inverted contrast? Yes 画像内では単粒子は背景に対して黒、参照像は背景に対して白なので、参照像のコントラストは反転されている。
References Are References CTF corrected? Yes 2D classificationの際にCTF補正を有効にしたため、Yes。
References Ignore CTFs until first peak? No 参照像を作る際の2D classificationでこれをYesにした時のみ、Yesにする。
autopicking Picking threshold 0.8 参照像ありの自動ピッキングでは、FOM(画像の各画素でのピッキングスコア相当)のマップが作られる。ピーク検出アルゴリズムがFOMマップからこの閾値以上のピークを検出して単粒子検出する。
autopicking Minimum inter-particle distance (A) 200 2つの単粒子間で許される最小の距離[Å]。逐次クラスタリングアルゴリズムが、互いにこの距離以内にある単粒子を除去する。有用な値は多くの場合単粒子の直径の50〜60%程度。
autopicking Maximum stddev noise -1 カーボン膜領域や大型のコンタミがある領域をピッキングしてしまうのを防ぐのに役立つ。背景領域の輝度値の標準偏差がこの値よりも大きくなるピークは無視される。1.0〜1.2あたりがおそらく有用です。-1とするとこの機能はオフになる。
autopicking Minimum avg noise -999 カーボン膜領域や大型のコンタミがある領域をピッキングしてしまうのを防ぐのに役立つ。背景領域の輝度値の平均がこの値よりも小さくなるピークは無視される。-0.5から0あたりがおそらく有用です。-999とするとこの機能はオフになる。
autopicking Write FOM maps? Yes 計算されたFOMマップをディスクへ書き出す。
autopicking Read FOM maps? No FOMマップの計算はスキップし、すでに計算済みのFOMマップをディスクから読み込んで使う。
autopicking Shrink factor 0 0とすると、Lowpass filter references (A)で指定した値が画像の最大の分解能となるよう自動で画像が縮小される。これにより計算が高速化され、メモリも節約でき、チュートリアルを早く終えるのに役立つ。0から1の間の値を指定すると、元画像のサイズをその値倍したサイズへ縮小される。shrink factor = 1 すなわち縮小なしの場合と比較して若干ピッキングの精度が低下することに注意。
autopicking Use GPU acceleration? Yes 参照像ありの自動ピックングはGPU化されている。
Helix 今回は無視する。
Running Number of MPI procs 1 FOMマップを書き出す場合はシングルMPIプロセスでしか動作しない。


他のジョブと区別するため、optimise_paramsというエイリアス名をつけておきます。【Run!】を押してジョブを実行しましょう。GPUを使えば30秒程度で終わります。 

このジョブで最も計算コストが高いのは、確率に基づくFOMマップ(各回転角度の各参照像と、画像の各場所での相互相関係数に類似した量です)を計算するところです。FOMマップ計算後は、閾値と最小ピーク間距離をパラメータとするピーク検出処理がありますが、これはFOMの計算よりはるかに高速です。ピーク検出のパラメータはtry&errorで最適化する必要があるため、このジョブはReferenceで指定すればFOMマップをディスクへ書き出して保存することができます。参照像1枚につき、ピッキング対象の画像と同じサイズのファイルが2つ生成されます。ディスクI/Oに過負荷がかかるのを防ぐために、autopickingはFOMマップの書き出しを有効にした場合はシングルMPIプロセスでしか動作しません。

一度FOMマップをディスクへ書き出せば、それを使ってピッキング用のパラメータチューニング過程を早く済ますことができます。

ではチューニングに取り掛かりましょう。まず、先の設定によるピッキングの結果を確認します。【Display:】ボタンでcoords_suffix_autopickを選択します。Pickを押せば画像とピッキング結果が表示されます。

rln-autopick-tune4.png


画像の表示サイズやシグマコントラストのシグマ値などの表示設定は、最後に行ったManual pickingジョブの設定が使用されています。表示設定を変更するには、Manual pickingジョブタイプへ移動し、Displayタブで適当な値を打ち込んで、RELION GUI左上のJobsメニューから"Save job settings"を選択すれば良いです。そのとき、以下のようにColorsタブにもパラメータを設定しておきましょう。FOM値に基づいてピッキング結果が色付けされます(下記設定ではFOM = 0(悪い)が赤、FOM = 1(良い)が青となります。

rln-autopick-tune5.png.png


再度先ほどのAutoPickジョブに戻り(※Finished jobsウィンドウから選択すること)、【Display:】ボタンでcoords_suffix_autopickを選択し、ピッキング結果を表示してみます。

rln-autopick-tune6.png


どのピッキング結果も青色(FOM〜1)で表示されているので、ピッキングされた単粒子の質は良いです。ただしピッキングされている単粒子の数は少なく、閾値が高すぎであることがわかります。このように、ピッキング結果を確認して、閾値をあげるべきか下げるべきか、単粒子間の最小距離をより狭めるべきか広げるべきかなど判断していきます。


画像は開きっぱなしで良いので、再度AutoPickジョブ(AutoPick/optimise_params)のautopickingタブに戻ります(Finished jobsリストから計算済みのジョブを選択するようにします。一度Runしたジョブは【Run!】ボタンが【Continue!】ボタンになっているので、それが目印になります。)

今度はFOMの設定を変えます。

  • Write FOM maps? : No
  • Read FOM maps? : Yes

その上で、Picking threshold, Minimum inter-particle distance, Maximum stddev noise, Minimum avg noiseパラメータを新しい値に変更します。

これにより、【Continue!】ボタンを押すと、ジョブはFOMマップの計算を省略して、すでにディスクに書き出されたFOMマップを読み込んで、わずか数秒でピーク探索、ピッキング処理を行います。処理が終了したら、先ほど開きっぱなしにしておいて画像ウィンドウ上で右クリックし、"Reload coordinates"をすれば新しい条件で計算されたピッキング結果が表示されます。(Reload coordinatesを押したらピッキング結果が全て消えてしまった場合、画像上で一度右クリックすると、新しい結果が読み込まれます)

rln-autopick-tune7.png


この過程を繰り返して最適なパラメータを見つければ良いです。


※RELIONチュートリアルでは最終的に以下のパラメータになっています。その時の5micsのピッキング結果画像を示しておきます。偽陽性が多いような気もしますがこれで進めます。

  • Pickng threshold: 0.0
  • Minimum inter-particle distance (A) : 100
  • Maximum stddev noise: -1
  • Minimum avg noise: -999

rln-autopick-tune-result1.png rln-autopick-tune-result2.png rln-autopick-tune-result3.png rln-autopick-tune-result4.png rln-autopick-tune-result5.png


パラメータを十分最適化できたら、ジョブタイプブラウザから新規にAuto-pickingジョブを立ち上げます(【Continue!】ではなく【Run!】ボタンが表示されます)。各ウィンドウには最適化したパラメータが入力された状態になっているはずです。

最適化されたパラメータを用いて全ての画像で自動ピッキングを行うためには、以下のようにします。

  • I/Oタブで、Input micrographs for autopick を CtfFind/job003/micrographs_ctf.star にする。
  • autopickingタブで、以下のようにパラメータ指定する。
    • Write FOM maps? : No
    • Read FOM maps? : No
  • ジョブのエイリアス名を template_based にしておく。

今度はFOMマップの書き出しを行わないので、複数のMPIプロセスを使って並列計算できます。今回のデータならば、シングルMPIプロセスでGPUを使えば1分程度で終了します。

(※本記事執筆者の環境では、全部で9,726粒子がピッキングされました。)


【Continue!】を押した時の挙動について

注意事項として、【Continue!】ボタンを押した時の挙動は、FOMマップの読み込み、書き出しの設定に依存して変わります。FOMマップの書き出しまたは読み込みのいずれかを有効にした状態で【Continue!】を押すと、全ての入力画像に対し再度ピッキングが行われます。しかしFOMマップの読み込みも書き出しも無効にした状態で【Continue!】を押すと、まだピッキングが実行されていない画像のみ処理が行われます。この機能は、電顕でデータを取得する先からオンザフライに処理を実行する際などに有効です(13.3節で詳しく説明します)。全て綺麗にやり直したければ、先ほど述べたようにジョブタイプブラウザから新規にAutoPickジョブを立ち上げて実行すれば良いです。


ピッキング結果の確認と手動での微調整

全画像に対する自動ピッキングが終了したら、【Display:】ボタンから"coords_suffix_autopick.star"を選択し、結果を確認しましょう。

rln-190522-2246-1.png


※やや偽陽性も目立ちますが、安定して単粒子をピッキングできています。


ここで全ての画像を確認し、偽陽性があれば手動で全て消す人々もいます。例えば、カーボン膜のエッジや高コントラストのアーティファクトは単粒子と間違って拾われていることが多いです。これをやるには、【Display:】ボタンを押して各画像を表示し、偽陽性の単粒子を中ボタンクリックして削除していきます。中ボタンを押しっぱなしにして動き回ればたくさんの偽陽性を楽に消し去ることができます。作業が終わったら、各画像ごとに、右クリックしてSave STAR with coordinatesを必ず実行し、座標ファイルを更新するようにしましょう。

また、【Display!】からout: logfile.pdfを選択すると、画像あたりの単粒子数、平均FOM値のプロットとヒストグラムを確認できます。


中間ファイルを削除してディスクを節約

結果に満足したら、書き出し済みのFOMマップを削除してディスク使用量を削減することもできます。そのためには、FOMを消したいAutoPickジョブをFinished jobリストから選択し、【Job actions】ボタンから"Gentle clean"を選択して、FOMマップを削除します。


参照像あり自動ピッキングにより生成されるSTARファイルについて

AutoPick/template_based/Movies/ 以下に、各画像ごとに <画像ルート名>_autopick.star というファイルが生成されます。中身は以下のような感じです。

# RELION; version 3.0.5

data_

loop_
_rlnCoordinateX #1
_rlnCoordinateY #2
_rlnClassNumber #3
_rlnAutopickFigureOfMerit #4
_rlnAnglePsi #5
  293.494119   440.241179            2     3.190388   320.000000
 2968.805900  2077.035306            1     2.989841   315.000000
 3228.435313  1614.217657            2     2.681165    50.000000
 3442.911785  2088.323542            3     2.176928   145.000000
 3273.588255  2031.882365            2     2.101359   200.000000
 3183.282372  1839.982364            2     1.891942   130.000000
(以下略)

LoGピッキングのSTARファイルと同様のラベル列ですが、今回は参照像を用いていることから、_rlnClassNumberにはピッキングに使われた参照像の番号、_rlnAnglePsiにはピッキングされた時の面内回転角度[度]が記録されています。なぜか_rlnClassNumberと_rlnAutopickFigureOfMeritの順序が入れ替わっていますが。


単粒子の切り出し

最後に、全画像から単粒子の切り出しを行いましょう。ジョブタイプブラウザからParticle extractionジョブを選択し、I/OタブのInput coordinatesに"AutoPick/template_based/coords_suffix_autopick.star"を指定しましょう(Browseボタンを使って指定すると良いです)。また、エイリアス名としてtemplate_basedを使います。

それ以外のパラメータは、前回のParticle extractionジョブ(LoG_based)と全く同じで良いです。

これで【Run!】すれば、以降の解析で使う、単粒子データセットを得られます。

【Display!】で切り出された単粒子を確認できます。

rln-190522-224925.png


このとき生成されるSTARファイルも、LoGピッキングの後のParticle extractionジョブのSTARファイルと同様です。参照像ありでピッキングしたため、_rlnAnglePsi, _rlnClassNumberがピッキング時の結果の値になっています。

1.9.1 The shrink parameter (shrinkパラメータについて)

処理を高速化するため、RELIONでは --shrink というコマンドラインオプションを用いて画像サイズを変更することができます。

最もシンプルな使い方は --shrink 0 ですが、このオプションはもう少し細かい使い方ができます。

--shrink のデフォルト値は 1.0 であり、これは何もしないことを意味します。

  • shrink = 0 のとき
    • --lowpassオプションで指定した分解能が画像の最大の分解能となるよう、画像が縮小されます。単粒子解析においてはこれが推奨されます。(ただしヘリカルな試料を扱う場合には向きません)


  • 0 < shrink <= 1 のとき
    • この場合は --shrink の値は画像の縮小率とみなされます。すなわち、新しい画像のサイズは shrink x (元の画像のサイズ) となります。


  • shrink > 1 のとき
    • 画像のサイズを shrink - (shrink mod 2) に変更します(mod は剰余演算です)。これにより、画像サイズが自分の望みのサイズ(偶数)にすることができます。


--shrinkによる新しい画像サイズの最大の空間周波数が、--lowpassで指定した空間周波数よりも小さくなってしまう場合、画像の中に存在しない分解能領域でローパスをかけようとしていることになるため、RELIONはnon-fatalな警告を発します。RELIONのautopickerは多数のFFT処理を行うので、FFTの主な落とし穴を回避するために、現在のRELIONでは画像のサイズが自動でFFTに適した値に調整されます。AutoPickを実行した際の標準出力に、これに関連したメッセージと、それを上書きするための方法が表示されます。

※ AutoPickのrun.outを見ると、冒頭近くで "+ Will use micrographs scaled to 340 pixels as requested. The largest prime factor in FFTs is 19" のような行があるので、それのことだと思います。ただ、FFT関連の値の上書き方法に関連したメッセージに関しては見つけられませんでした。

1.10 Particle sorting (単粒子の並べ替え)

RELIONには、各単粒子と、それとアラインされたCTF補正済み参照像の間の差異の大きさに基づいて、単粒子を並べ替える機能が実装されています。並べ替えた結果を利用して、コンタミや異常な単粒子画像を捨て去ることができます。

並べ替え機能はParticle sortingジョブタイプで提供されており、単粒子とそれに対応する参照像のセットがあれば実行することができます。すなわち、Auto-picking, 2D classification, 3D classification, 3D auto-refineのいずれかを行なった後なら実行できます。ここでは先ほどのauto-pickingの結果に対し並べ替えを実行してみます。

まず、以下のようにパラメータを入力します。(※都合によりジョブ番号が異なっています。また、RELION3の初期版ではReferenceタブの名称がCTFタブになっているようです。)

rln-190523-102627.png

rln-190523-102727.png

単粒子並べ替えのパラメータ
タブ パラメータ名 説明
I/O Input particles to be sorted Extract/template_based/particles.star 今回はParticle extractionジョブで切り出された単粒子を並べ替えるので、これを指定。2D classification, 3D classification, 3D auto-refineジョブの結果で単粒子を並べ替える場合には、最後のイテレーションの_data.starファイルを指定する。
I/O References from model.star 入力単粒子が2D classification, 3D classfication, 3D auto-refineジョブ由来の場合は、ここに最終イテレーションの_model.starファイルを指定する。今回は空白。
I/O OR autopicking references Select/templates4autopick/class_averages.star 入力単粒子のピッキングに使用した参照像ファイルを指定する。
References Pixel size in references (A) -1 -1とすると、参照像のピクセスサイズ[Å/pixel]は単粒子画像のピクセルサイズと同じとみなされる。
References Are References CTF corrected? Yes 説明済み
References Ignore CTFs until first peak? No 説明済み


エイリアス名としてafter_autopickを指定しておきましょう。MPIプロセスを増やして並列化することもできますが、今回はシングルMPIプロセスでも数秒で処理が完了します。【Run!】を押し、ジョブを実行します。

ジョブが完了すると、Sort/after_autopick/ に particles_sort.star というSTARファイルが生成されています。これは入力したParticle extractionジョブのSTARファイルに対し、_rlnParticleSelectZScore という新たなラベルが追加されたものになっています。

Particle sortingにより新たに追加される値
ラベル名 説明
_rlnParticleSelectZScore (恐らく)参照像と単粒子画像の差異を表す指標のZスコア。Zスコアとは、その値が平均値から標準偏差の何倍離れているかを示す値。例えばZスコア4なら、平均値から4標準偏差離れた値で、外れ値であることが示唆される。Zスコアが大きいほど参照像との差異が大きいと考えれば良い。


【Display:】でout: particles_sort.starを選択してRelion display GUIを表示し、rlnParticleSelectZScoreでソートして単粒子画像を表示してみます。

rlnParticleSelectZScoreにより昇順(小さいものから順に)で表示すると、最初の100枚は以下のような感じです。

rln-190523-123555.png


rlnParticleSelectZScoreにより降順で表示すると、最初の100枚は以下のような感じです。

rln-190523-123717.png


rlnParticleSelectZScoreが大きいほど、コンタミなど蛋白質でないものが目立つことが確認できます。


ということで、このrlnParticleSelectZScoreを用いて、単粒子画像の選別をすることができます。これにはSubset selectionジョブタイプを使用します。

まずSubset selectionジョブタイプで、I/Oタブの"OR select from particles.star"に対し、rlnParticleSelectZScore付きのSTARファイルを指定します(Sort/after_autopick/particles_sort.star)。I/Oタブの他の項目は空欄にします。ジョブのエイリアス名として after_sorting と入力しておきましょう。

【Run!】を押すとRelion display GUIが立ち上がります。Sort images onにチェックを入れ、rlnParticleSelectZScoreを指定します。全ての単粒子を表示するために、Max. nr. imagesに-1を指定しておきます。左上からrlnParticleSelectZScoreが小さい順にソートされて表示されます。下の方に(rlnParticleSelectZScoreが大きい、すなわち参照像と似ていない)スクロールしていくと、あるあたりから氷のコンタミやよくわからない蛋白質以外の何かが目立ち始めるところがあると思います。その境界あたりに位置する単粒子を左クリックで選択し、右クリックして"Select all above"を選択すれば、そこから上の全ての単粒子が選択されます。もちろん個別に選択したり除外したりすることもできます。この状態で右クリックメニューから"Save STAR with selected images"を押せば、選択された単粒子のみで構成されるparticle STARファイルを書き出すことができます。これは、rlnParticleSelectZScoreが高い単粒子をデータから捨て去り、データセットをクリーンにする意味があります。

rln-190523-125640.png


(※本記事執筆者の場合、9,726粒子のうち、9,514粒子が選択されました。212粒子を捨てることになります。)


上記では捨て去る単粒子たちの境界を目視で決定しましたが、データセットが大きい場合などは、上記のように全ての単粒子を表示して目視で境界を決めるのは時間がかかりますし、メモリ消費も激しくなります。そこでRELION3以降は、STARファイル内の何らかのラベルの値で閾値処理することにより単粒子のサブセットを選択する機能が実装されています。本チュートリアルではこの方法は撮りませんが、rlnParticleZScoreを使って選択する場合について、その方法を説明しておきます。

まず、rlnParticleSelectZScoreの閾値を決めます。autopickの単粒子並べ替えを行なったジョブ(Sort/after_autopick)をFinished jobsパネルから選択し、【Display:】から out:logfile.pdfを選択します。すると下記のように、rlnParicleSelectZScoreのヒストグラムとプロットが表示されます。


ヒストグラム

rln-190523-130814.png


プロット(横軸は単粒子の番号)

rln-190523-130932.png


ヒストグラムは潰れていてわかりにくいですが、プロットの方では、ほとんどの単粒子がZスコア1あたりに収まっており、それ以外は外れ値(捨ててもいい画像)とみなして良さそうです。0.8を閾値とすることにします。

Subset selectionジョブタイプで、I/Oタブでparticles_sort.starを選択し、Subsetsタブに以下のように入力します。

rln-190523-132327.png

これで【Run!】すれば、メタデータrlnParicleSelectZScoreの値が0.8以下の単粒子のみが選択され、STARファイルに書き出されます。

2 Reference-free 2D class averaging (参照像無しでの2Dクラス平均化)

ほとんどの場合、参照像無しでの2Dクラス平均化を用いて、"悪い"単粒子を捨てる工程を行います。前節までで行なってきた通り、自動ピッキングの結果を手動で編集したりrlnParticleSelectZScoreによる並べ替えに基づくサブセット選択などにより、悪い単粒子を捨てるよう努力はしたわけですが、ほとんどの場合これは十分ではなく、目的物以外が写った画像が単粒子に紛れ込んでいます。そのような画像は2Dクラス平均にかけるとうまく平均化されないがために、見た目が悪く粒子数も少ないクラスに割り振られることが多いです。そうした見た目の悪いクラスに属する単粒子を捨てることにより、効率よくデータセットをクリーンアップできるのです。

2.1 Running the job

ほとんどのオプションは、参照像あり自動ピッキングのための参照像作りの際に行なった2D classfiicationと同じです。

一部、I/OタブとOptimisationタブのパラメータを以下の通り変更します。また、エイリアス名としてafter_sortingを付けておきましょう。

  • I/O : Input images STAR file
    • Select/after_sorting/particles.star
  • Optimisation : Number of classes
    • 100
      • 単粒子数が増えたのでクラス数も増やします。(概ね10,000粒子程度かと思いますので、100クラスなら単純計算で1クラスあたり100粒子を期待でき、適切なレンジです。)


(Sjorsらの環境では)4GPU、5MPIプロセス、6スレッド/MPIプロセスで20分程度かかります。コーヒーでも飲んで待ちましょう。

小話

ジョブが終了したら、Subset selectionジョブタイプに行き、I/Oタブの"Select classes from model.star"に、最後のイテレーションの Class2D/after_sorting/run_it025_model.star を指定します。ジョブのエイリアス名は class2d_aftersort としておきましょう。

【Run!】ボタンを押し、rlnClassDistributionでリバースソートするなどしてクラス平均像を表示します。CTF補正を行なっているので背景に対し蛋白質が白くなっています。

rln-190523-143834.png


上の方に、蛋白質の内部構造まで見えるきれいなクラスが固まっています。それらを選択します。

rln-190523-144037.png


(※上記9クラスだけ選択してみました。)

選択し終えたら、右クリックメニューから"Save selected classes"を選択し、選択結果を保存します。 (※入力が 9,514粒子に対し、5,816粒子が選択されました。残り2,698粒子は捨てることになります。)


"良い"クラス平均像とは

良いクラス平均像は、蛋白質の原子モデルにローパスフィルタをかけて投影したのと似たような見た目をしています。蛋白質ドメインの内部構造が見え、蛋白質周囲の溶媒領域は一様であるのが理想です。良い2Dクラス平均像が得られる場合、3Dマップも高い分解能で再構成できることが期待できます。蛋白質から溶媒領域に向かって放射状のストリークが見える場合、オーバーフィッティングが発生している兆候です。その場合は、2D classificationで Limit resolution E-step を使い、アラインメント時に使う分解能を制限して再計算してみるのが良いでしょう。


クラス選択の時の注意点 - "Einstein from noise"

自動ピッキングにおいて低い閾値を使用した場合、"Einstein-from-noise"クラスに注意しましょう(※参照像とのアラインメントにより、ノイズ画像から参照像とよく似た模様が浮かび上がってしまう現象。R.Henderson, PNAS(2013))。そのようなクラス平均像の特徴として、ピッキングに使用した参照像の低分解能のゴースト的なものの上に、高分解能のノイズが蓄積しています。それらは選択しないように注意してください。(※上で示した100個のクラス平均像の中には、そうしたクラスは無いように見えます。)


何度か2Dクラス分類と選別を繰り返す

2Dクラス分類を用いて単粒子を選別する工程を何度か繰り返すこともできます。つまり、2D classificationしてSubset selectionした後、もう一度(あるいは何度か) 2D classification & Subset selectionを行う、ということです。


2Dクラス分類の後は...

また、今回のチュートリアルではやりませんが、2D classificationを行なった後は、原則として単粒子の並び替え(Particle sorting)を行い、異常な単粒子を除去する工程を行いましょう。

2.2 Analysing the results in more defail

※チュートリアルをお急ぎの時は読み飛ばしてください。

※元チュートリアルでは、STARファイルの中身などについて解説されています。本記事ではこちらに書いてあります → 2D classificationジョブのSTARファイル

2.3 Making groups (グループ分けする)

※チュートリアルをやる上では必要ありません。お急ぎの時は読み飛ばしてください。

RELIONでは単粒子がグループ分けされます。各グループごとに以下のような統計量が計算されます。

  • ノイズのパワースペクトル。
  • 輝度値のスケール因子。(データセット内での氷の厚み、デフォーカス、コンタミ量のばらつきなどに由来するS/Nのばらつきを知ることができます。)

デフォルトでは、同じ画像から切り出された単粒子がグループ化されます。(切り出し元画像の枚数分だけグループができます)

画像あたりの単粒子切り出し数が十分大きければ、これで問題ありません。しかし、高倍率で撮影した場合や、試料が希薄な場合や、単粒子選別の結果単粒子が大きく減ってしまった場合などはグループあたりの単粒子数が少なくなってしまいます。すると輝度値のスケール因子(およびノイズのパワースペクトル )の推定が不安定となります。

一般的には、少なくとも10〜20粒子/グループあることが推奨されます。2D/3D classification (& subset selection)をすると単粒子が減りますので、この値に注意してください。


もし画像あたりの単粒子数が少なくなりすぎた場合には、複数の画像由来の単粒子を1つのグループへまとめてしまうことを推奨します。このためには、Subset selectionジョブタイプを使用します。

Subset selectionジョブタイプで、以下入力します。

  • I/Oタブの"Select classes from model.star" に *_model.starファイルを指定
  • Class optionsタブの"Regroup the particles"にYesを指定
  • Approximate nr of groups に適当なグループ数を入力

【Run!】すると単粒子が概ねApproximate nr of groupsで指定したグループ数に分けられます。(概ねというのは、指定した値とはやや異なるグループ数へと別れる場合もあるためです。RELIONにまかせましょう。)

今回のチュートリアルでは、画像あたり十分な数の単粒子があるため、この工程は不要です。


RELIONにおけるグループは、他のソフトウェアにおける"デフォーカスグループ"とは全く別物です。RELIONは常に単粒子ごとにCTF補正を行います。どのグループに単粒子が所属しているかはCTF補正と関係ありません。

3 De novo 3D model generation (De novoに3Dモデルを生成する)

RELIONでは、確率的勾配降下法(SGD; Stochastic Gradient Discent)を用いて、2Dの単粒子画像から de novo に3D初期モデルを生成できます(※事前情報なしで2Dから3Dを推定するということ)。RELION3.0以降では、cryoSPARCにおける実装にほぼ準拠した実装になっています。色々な投影方向の単粒子があり、かつ2Dクラス分類において良質なクラス平均像が得られている場合、このアルゴリズムによって、3D classificationや3D auto-refinementの初期モデルとして使用可能な低分解能3Dモデルを得られる確率が高いです。

RELION3.0以降は、実装が改善されたことにより、2Dクラス分類で得られたクラスの中からランダムに一部の単粒子を選び出して使うといった手続きは不要になりました。現在のアルゴリズムは大変頑健に動作するので、全単粒子を入力として使っても問題ありません。

※cryoSPARCのSGDによるde novo 3D構造決定のアルゴリズムの参考文献

  • Ali Punjani et al., Nature Methods (2017)
    • cryoSPARCの元論文です。
  • https://arunabh98.github.io/reports/cryosparc.pdf
    • ドキュメントの素性がよくわかりませんがcryoSPARCにおけるde novo 3D構造決定とリファインメントの概要が解説されています。

3.1 Running the job

ジョブタイプブラウザから3D initial modelジョブタイプを選択し、以下のようにパラメータを入力します。(※都合によりジョブ番号が違っています。また、ComputeタブとRunningタブの内容は各実行環境に合わせて変更してください。)

※SGDのイテレーションは初期(initial)、中期(in-between)、後期(final)の3期に別れます。それぞれ処理が少し異なります。

rln-init-1.png

rln-init-6.png

rln-init-2.png

rln-init-3.png

rln-init-4.png

rln-init-5.png


3D initial modelジョブのパラメータ
タブ パラメータ名 説明
I/O Input images STAR file Select/class2d_aftersort/particles.star 構造推定したい単粒子セットを記述したSTARファイル。数千〜1万粒子程度あれば十分。
CTF Do CTF-correction? Yes 説明済み
CTF Have data been phase-flipped? No 説明済み
CTF Ignore CTFs until first peak? No 説明済み
Optimisation Number of classes 1 推定する3D構造(クラス)の個数。SGDによるde novo構造決定では、一度に複数の3D構造を推定することが可能。SGDのイテレーションの初期では1種類の3D構造推定が行われ、それを出発点として徐々に異なる3D構造へと変化していく。2個以上を指定すると、データセット内に含まれているsub-optimalな単粒子を振り分けるための'sink'(流し)として働き、構造決定に役立つことがある。その場合はRunningタブのadditional argumentsに --sgd_skip_annel を指定し、イテレーション中期におけるアニーリング処理を無効化することも有用。今回は計算時間短縮のために1クラスだけ推定する。
Optimisation Mask diameter (A) 200 これまで通り。
Optimisation Flatten and enforce non-negative solvent? Yes (3Dで? )球形のマスク処理を施し、マスクの外側を非負の値にする。
Optimisation Symmetry C1 対象蛋白質の構造の空間群(対称性)。空間群が不明の場合は、C1でとりあえず再構成するのが良いだろう。また、高次の対称性を持つ場合は、本来の空間群でSGDするよりはC1でSGDした方が構造を得やすい。今回用いているデータは素晴らしいもので、本来の空間群D2でもうまくいく。しかし、最初にC1で3D構造推定し、その後D2でリファインメントする流れを説明するため、今回のチュートリアルではC1でまずSGDすることにする。
Optimisation Initial angular sampling 15 degrees イテレーション初期における投影角サンプリング間隔。オイラー角rot, tiltのサンプリング処理にHealPixライブラリを利用している都合上、サンプリング間隔は離散的ないくつかの値しか選択できない。低分解能で3D初期構造を得る際は、3D classificationやrefinementと比べて荒いサンプリング間隔で良い(15度など)。イテレーション中期と後期ではサンプリング間隔は単粒子のサイズと分解能に応じて自動調整される。
Optimisation Offset search range (pix) 6 2D classificationで説明済み。
Optimisation Offset search step (pix) 2 2D classificationで説明済み。
SGD Number of initial iterations 25 初期イテレーション数。分解能はInitial resolution (A)でカットされ、ミニバッチサイズは Initial mini-batch sizeで固定、複数3Dモデルは全て同一、の条件でSGD。50で多くの場合うまくいく。ダメなら増やす。
SGD Number of in-between iterations 100 中期イテレーション数。分解能カットオフはイテレーションに応じてInitial resolutionからFinal resolutionまで線形に増やす。ミニバッチサイズについても同様。複数3Dモデルがある場合、アニーリング処理により徐々に各モデルが変化するようにする。200で多くの場合うまくいく。複数の構造にうまく別れなかったり、正しい構造が出ない場合には増やすと良い。
SGD Number of final iterations 25 後期イテレーション数。分解能はFinal resolution (A)でカットされ、ミニバッチサイズは Final mini-batch sizeで固定、複数3Dモデルは各々異なるに任せる、の条件でSGD。50で多くの場合うまくいく。複数の構造にうまく別れなかった場合は増やすと良い。
SGD Write-out frequency (iter) 10 3Dモデルを何イテレーションごとに保存するか。
SGD Initial resolution (A) 35 初期イテレーションの分解能カットオフ[Å]
SGD Final resolution (A) 15 後期イテレーションの分解能カットオフ[Å]
SGD Initial mini-batch size 100 初期イテレーションの単粒子画像ミニバッチサイズ[枚]
SGD Final mini-batch size 500 初期イテレーションの単粒子画像ミニバッチサイズ[枚]
SGD SGD increased noise variance half-life -1 -1とすると無効。正の値を指定すると、ノイズの分散の一番最初の推定値は8倍され、処理が進むに連れて徐々に小さくなる。このオプションで指定した枚数の単粒子を処理した時にちょうど50%となるようにする。SGDがうまくいかない時に試すとうまくいく可能性がある。1000枚あたりにするとうまく行くようである。最初にノイズの分散を何倍するかは --sgd_sigma2fudge_ini オプションで変更できる。


典型的には、SGDタブの内容はとりあえず変更せずに走らせます。多くの場合はそれでうまくいきます。今回はチュートリアルを早く進めるため、イテレーション回数をデフォルトの半分にしています。


エイリアス名として symC1 を指定しておきましょう。RunningとComputeの設定は適宜変更して、【Run!】してください。4 GPU、5 MPIプロセス、6 スレッド/MPIプロセスの条件で15 分程度かかります。お茶でも飲んでゆっくりしてください。

3.2 Analysing the result

出力された3Dモデル(InitialModel/symC1/run_it150_class001.mrc)を、UCSF Chimeraなどの3Dビューアで見てみましょう。

$ chimera InitialModel/symC1/run_it150_class001.mrc


rln-190525-1.png


以下はPDB 6DRVをフィッティングした様子です。

rln-190525-2.png


可視化により対称性があることに気づいたら、対称性の軸が3D座標系のx, y, z軸と(RELIONの慣例にならって)揃える必要があります。RELION3では、これを助けるためのプログラムが新たに実装されています。コマンドラインで以下のコマンドを実行してみましょう。

$ relion_align_symmetry --i InitialModel/symC1/run_it150_class001.mrc --o InitialModel/symC1/run_it150_class001_alignD2.mrc --sym D2

※ 本記事執筆者の場合、出力は以下のようになりました。D2対称軸を探索した上で、それがx, y, z軸と揃うように回転を施してくれているようです。

 Reading map: InitialModel/symC1/run_it150_class001.mrc
 The input box size: 64
 Using the pixel size in the input image header: 3.54 A/px
 Center of mass: x= -0.363242 y= 1.31628 z= 1.18839
 Re-centred to the centre of the mass
 Downsampled to the working box size 64 px. This corresponds to 3.54 A/px.
 Generating 400 projections taken randomly from a uniform angular distribution.
 Searching globally ...
   8/   8 sec ............................................................~~(,_,">
 The best solution is ROT = 309.953 TILT = 121.807 PSI = 232.434

 Refining locally ...
   3/   3 sec ............................................................~~(,_,">
 The refined solution is ROT = 305.953 TILT = 121.807 PSI = 228.434

 Now rotating the original (full size) volume ...

 The aligned map has been written to InitialModel/symC1/run_it150_class001_alignD2.mrc


UCSF chimeraやrelion_displayで対称軸がx, y, z軸に揃えられたかどうか確認しましょう。

※ここではEMAN2のe2display.pyのXYZ projectionを使ってみました。

アラインメント前 (x, y, z軸方向への投影)

rln-align-pre-x.png rln-align-pre-y.png rln-align-pre-z.png


アラインメント後 (x, y, z軸方向への投影)

rln-align-after-x.png rln-align-after-y.png rln-align-after-z.png


実際に適切にアラインメントされていることが確認できたら、以下コマンドによりD2対称性を課す(対称性に従い平均化する)ことができます。

$ relion_image_handler --i InitialModel/symC1/run_it150_class001_alignD2.mrc --o InitialModel/symC1/run_it150_class001_symD2.mrc --sym D2


対称性を課す前と課した後で比較してみましょう。おおよその見た目は類似しているはずです。(※大きく違っているとしたら何かが間違っています。)

※ relion_displayで3Dモデルを開くと、z軸方向のスライスが並べて表示されます。

$ relion_display --i InitialModel/symC1/run_it150_class001_alignD2.mrc &
$ relion_display --i InitialModel/symC1/run_it150_class001_symD2.mrc &

※左が対称性を課す前、右が対称性を課した後です。

rln-190525-3.png

※それほど変化は見られません。間違った対称性を課していないと確認できました。


※3Dでも確認してみます。左の灰色がD2対称性を課す前、右の黄色が課した後です。

平均値で等高面表示

rln-190525-4.png

※対称性に基づく平均化がされたことで、溶媒部分のノイズが軽減されていることがわかります。


平均値 + 4.5σで等高面表示

rln-190525-5.png

※高い等高面で見ると、対称性を課す前は細かいところで非対称な密度分布が見られますが、対称性を課すことでそれが解消されます。


STARファイルについて

3D initial modelで生成されるSTARファイルは、項目としては2D classificationで生成されるSTARファイルとほとんど同じになります。

3D推定が行われたため、_data.starで _rlnAngleRot, _rlnAngleTiltなど、投影角を表すパラメータに値が割り当てられているのがわかると思います。

比較的大きく変わっているのは_sampling.starで、3D投影角のサンプリングに関する設定が追加されています。


_grad001.mrc について

※(たぶん)SGDでは3Dモデルに対する対数事後確率分布の勾配が計算され、その勾配を用いて対数事後確率分布が大きくなる方向に3Dモデルが更新されていくようです。その時の勾配をマップとして保存してあるのが_grad001.mrcなるmrcファイルかと思います。こういう方法で3Dモデルを最適化していくということとCTF補正を行うことがどう両立しているのか、よくわかりません。

4 Unsupervised 3D classification (教師無し3Dクラス分類)

どんなデータであっても、そこに含まれている蛋白質の構造は多様です。全く同一なんていうことはありません。問題は、それをどこまで許すか、ということです。RELIONでは、複数の参照像を用いた3Dリファインメントにより、教師無しでの3Dクラス分類を行うことができます。

4.1 Running the job

教師無し3Dクラス分類は3D classificationジョブタイプから実行することができます。以下のようにパラメータを入力します。(※ジョブの番号が都合により異なっています。)

rln-cls3d-1.png

rln-3dcls-2.png

rln-3dcls-3.png

rln-3dcls-4.png

rln-3dcls-5.png

rln-3dcls-6.png

rln-3dcls-7.png


3D classificationジョブのパラメータ
タブ パラメータ名 説明
I/O Input images STAR file Select/class2d_aftersort/particles.star 3Dクラス分類の対象とする単粒子セット。ここではInitialModel生成によりパラメータがリファインメントされた単粒子セット(run_it150_data.star)ではなく、InitialModel生成の対象になった単粒子セットを指定している。
I/O Reference map InitialModel/symC1/run_it150_class001_symD2.mrc D2対称性を課した3Dマップを参照像として指定。このマップはまだRelionのデータパイプラインにImportしていないため、Browseボタンでは現れてこない。ファイル名を先の通り手打ちするか、ImportジョブでNode type"3D reference(.mrc)"としてインポートし、その後ここでBrowseボタンを使いインポートするかせねばならない。また、今回のジョブはC1空間群としてリファインメントするので、参照像としてInitialModel/symC1/run_it150_class001.mrc を指定しても良い。ただ、ここで正しい対称性を課したマップを参照像にしておいた方が後々楽である。
I/O Reference mask (optional) 今回は空白。このオプションは、例えばリボソームを対象にしたプロジェクトで、largeサブユニットもしくはsmallサブユニットだけを切り出すようなマスクを使ってfocusedリファインメントするといった場合に使う。この欄が空白だと、Optimisationタブのparticle diameter (A)で指定した値を直径とするソフトな球形マスクが適用される。それが3D分類におけるバイアスが最も少ないマスクでもある。
Reference Ref. map is on absolute greyscale? Yes 参照像がこのデータセットを用いて再構成されたものであるなら、それはabsolute greyscaleである(Inputに指定した単粒子データセットとReferenceに指定した参照像とでグレースケールが一致している)。参照像がこのデータセット以外のデータを用いて用意されたものなら、absolute greyscaleでは無いので、Noとする。RELIONの確率モデルにはデータと参照像の自乗誤差項が含まれるため、参照像の輝度値のスケールがデータのそれと異なっているとおかしなことになる。(※例えばPDB構造を密度マップ化したものや、RELIONまたはXMIPP以外のソフトで再構成されたマップや、同一の蛋白質だが別のデータセットからRELIONで3D再構成されたマップなど。)
Reference Initial low-pass filter (A) 50 高分解能の初期構造は使うべきでは無い(リファインメントにバイアスを導入してしまうため)。Sjors Scheres, Methods in Enzymology(2010)で述べられている通り、初期構造にはできるだけ強いローパスをかけた方が良い。リボソームの場合(Sjorsらは) 70Åをよく使う。それより小さい蛋白質なら40Å〜60Åをよく使う。
Reference Symmetry C1 我々はこの蛋白質がD2対称であると知っているわけだが、最初の3Dクラス分類においては対称性を課さない方が基本的には良い。そうすることにより、対称性とマッチしない変な単粒子がゴミ溜め('sink')クラスに振り分けられやすくなる。
CTF Do CTF-correction? Yes 説明は省略
CTF Has reference been CTF-corrected? Yes SGDにおいてCTF補正有りで再構成したので、Yes
CTF Have data been phase-flipped? No 説明は省略
CTF Ignore CTFs until first peak? No 説明は省略
Optimisation Number of classes 4 3Dクラス数(k)。最初のイテレーションでデータがランダムにk個のサブセットに分割され、与えられた単一の初期参照像に対しそれぞれアラインメント、3D再構成されることでk個の異なる初期3Dモデルが生成され、クラス分類が進む。クラス数を大きくすればより対象の多様性に対処する能力が上がる。計算コスト(CPU時間および消費メモリ)はクラス数に比例して増える。
Optimisation Regularisation parameter T 4 説明は省略
Optimisation Number of iterations 25 説明は省略
Optimisation Use fast subsets (for large data sets) No 説明は省略
Optimisation Mask diamter (A) 200 説明は省略
Optimisation Mask individual particles with zeros? Yes 説明は省略
Optimisation Limit resolution E-step to (A) -1 説明は省略
Sampling Perform image alignment? Yes 説明は省略
Sampling Angular sampling interval 7.5 degrees 説明は省略
Sampling Offset search range (pix) 5 説明は省略
Sampling Offset search step (pix) 1 説明は省略
Sampling Perform local angular searches? No Yesにすると、網羅的な投影角探索は行わず、現状の投影角の周りのみで局所的に投影角推定を行う。(投影角度の事前確率分布が、現状の投影角を中心とし、Local angular search range の 1/3 の値を標準偏差としてもつガウス分布に設定される)


※すでに説明済みのパラメータについては説明を省略しました。多くは2D classificationのパラメータと被っています。

※Compute, Runningは各々の環境に合わせて設定してください。


Samplingタブについては、普通はデフォルトから何も変更しなくて大丈夫です。(正20面体ウイルス粒子のように、巨大で対称性の高い粒子の場合のみ、この時点で投影角サンプリングレートを3.7度にするなどします。)

MPIのプロセス数とスレッド数に関し、2D classificationで説明したように、3D classificationジョブはメモリ消費が2D classificationよりも増えるため、より多くのスレッド数を使うことが多いです。しかし今回のチュートリアルの場合、データの数がそれほど多く無いので、メモリ消費はあまり問題になりません。

ジョブのエイリアス名として、first_exhaustiveを指定しておきましょう(これが最初の3Dクラス分類で、かつ角度探索がexhaustive(網羅的)だからです)。

【Run!】を押してジョブを実行してください。(Sjorsらの環境では)4 GPU、5 MPIプロセス、6 スレッド/MPIプロセスで、およそ10分かかります。


クラス分類の結果を確認する際、例えばUCSF Chimeraなどで3Dマップを等高面表示するだけでなく、スライス画像シリーズとして見ることは非常に大切です。スライス画像シリーズで可視化すると、マップの中で構造が一定していない部分(ヘテロな部分)がぼやけた、あるいはストリーク状の領域として現れ、構造分類の良し悪しが把握しやすくなります。また、溶媒領域が一様であるか否かも見やすくなります。【Display】で out: run_it025_class<クラス番号>.mrc を選択し、スライス表示してみます。(※多分z軸に垂直な平面を、z座標が小さい方から順にみています)

class001

rln-190525-6.png


class002

rln-190525-7.png


class003

rln-190525-8.png


class004

rln-190525-9.png


※ class001, class002はとても綺麗です。どちらもそれぞれ似ているのであまり分かれている意味がある気がしませんが、特にぼやけたところもないので大きく違う色々な構造が入り混じっていることはなさそうです。他のクラスはゴミの集まりのように見えます。3D classificationには乱数ガチャ要素があるため、人によってどのクラスが良いかは変わります。

※ class002を拡大した図が以下です。綺麗です。

rln-190525-10.png


3Dでも確認してみましょう。最も綺麗な3Dマップに対し、他の3Dマップをアラインメントして表示すると良いです。リファインメントの途中でわずかにそれぞれ独立に回転してしまっているためです。UCSF Chimeraで、run_it025_classのclass001〜class004の.mrcを全て開き、Tools > Volume Data > Fit in Mapを使います。

rln-chimera-align-1.png


すると以下のFit in Mapウィンドウが出てきます。(本記事執筆者の場合は)class002が最も良いマップなので、class001, class003, class004を、それぞれclass002に対しFitします。

rln-chimera-align-2.png


Tools > Structure Comparison > Tile Structures すると、各モデルがタイル状に並べて表示されます。Favorites > Side Viewと押すとViewingパネルが表示されます。Rotationタブを開き、center of rotation methodをIndependentにすると、各モデルがそれぞれの中心の周りで回転してくれるので、みやすくなります。

rln-190525-11.png


※一番綺麗なclass002にPDB 6DRVをフィッティングしてみました。一部のαヘリックスがが筒として見え始め、βシートもシート状に見え始めています。

rln-190525-12.png


クラス選択

2D classificationの後と同じように、Subset selectionジョブを使って、良い3Dクラス(達)に属する単粒子のみを選択してサブデータセットを作ります。Subset selectionジョブタイプへ行き、I/OタブのSelect classes from model.starに Class3D/first_exhaustive/run_it025_model.starを指定します。ジョブのエイリアス名として、class3d_first_exhaustiveと指定しておきましょう。

【Run!】を押して画像を表示します。以下のように、左から順にclass001からclass004まで、3Dモデルの中央スライスが表示されます。選択したいクラスを左クリックで選択し(複数可。選択すると赤枠で表示。)、選択し終えたら右クリックで"Save selected classes"を押し、選択結果を保存します。

rln-190525-13.png


※本記事執筆者の環境では、class001とclass002を選択し、それに属する 5,355 粒子が保存されました。選択前の粒子数は 5,816 粒子だったので、461粒子を捨てたことになります。ヘテロな構造を見つけ出すというより、ゴミ落とし的に働いた感じです。

4.2 Analysing the results in more detail

チュートリアルで先を急ぐ場合は読み飛ばして構いません。

3D classificationジョブタイプで出力されるファイルは、2D classificationと基本的には同じです。(実は内部的には2D classificationと3D classificationは同じコードを用いています。)違いはどこにあるかというと、2D classificationでは各クラスの平均画像がまとめてmrcスタックとして出力されるのに対し、3D classificationではそれぞれのクラスの密度マップが個別のMRCファイルとして保存される点です(例: run_it025_class00?.mrc)。

以前と同様、所属する単粒子が少ない小さなクラスは、より強いローパスフィルタがかかっています。_model.starファイルの中のdata_model_class_N (N = 1, 2, ..., K)データブロックに、それぞれのクラスのspectral signal-to-noise ratio (SSNR)が記録されています。

この辺りで、STARファイルの中から任意のデータを抽出するのに便利な2つのスクリプトについて紹介しておきましょう。

コマンドラインで、以下を試してみてください。

$ relion_star_printtable Class3D/first_exhaustive/run_it025_model.star data_model_class_1 rlnResolution rlnSsnrMap

出力

0.000000 42504.825790
0.004414 5516.278366
0.008828 2994.762452
0.013242 1367.467508
0.017655 712.952953
0.022069 399.238037
0.026483 436.627611
0.030897 270.392675
0.035311 205.188045
0.039725 141.329182
0.044138 106.747491
0.048552 60.015935
0.052966 26.527252
0.057380 13.794756
0.061794 8.779093
0.066208 12.363512
0.070621 15.517685
0.075035 15.879017
0.079449 9.766969
0.083863 3.001588
0.088277 0.447204
0.092691 0.085640
0.097105 0.031358
0.101518 0.012671
0.105932 0.005392
0.110346 0.002789
0.114760 0.001623
0.119174 8.085094e-04
0.123588 4.006819e-04
0.128001 2.385425e-04
0.132415 1.600752e-04
0.136829 0.000000
0.141243 0.000000


上記コマンドにより、class001について、空間周波数(rlnResolution [1/Å]とMAP推定されたSSNR(rlnSsnrMap)が出力されました。これを自分の好みのプロットプログラムに渡してプロットすることができます。もしgnuplot(有名なデータプロットツール)がインストールされていれば、以下コマンドを打つと、

$ relion_star_plottable Class3D/first_exhaustive/run_it025_model.star data_model_class_1 rlnResolution rlnSsnrMap

空間周波数に対するSNRをプロットすることができます。

(※なぜかよくわかりませんが私の環境ではうまくプロットできませんでした。代わりにmatplotlibでプロットしたものを載せておきます。空間周波数に対してSNRが指数関数的に減少していることがわかります。)

rln-ssnr-1.png


また、25イテレーションで計算が収束したか否か確認するには、以下コマンドで確認できます。

$ grep _rlnChangesOptimalClasses Class3D/first_exhaustive/run_it???_optimiser.star

※ ???はシェルのワイルドカードです。任意の3文字(数字含む)にマッチします。 ※ エイリアス名がfirst_exhaustive_try2になっていますが実験の都合なので気にしないでください。

Class3D/first_exhaustive_try2/run_it000_optimiser.star:_rlnChangesOptimalClasses                             9.999999e+06
Class3D/first_exhaustive_try2/run_it001_optimiser.star:_rlnChangesOptimalClasses                                 0.987105
Class3D/first_exhaustive_try2/run_it002_optimiser.star:_rlnChangesOptimalClasses                                 0.688618
Class3D/first_exhaustive_try2/run_it003_optimiser.star:_rlnChangesOptimalClasses                                 0.437758
Class3D/first_exhaustive_try2/run_it004_optimiser.star:_rlnChangesOptimalClasses                                 0.359869
Class3D/first_exhaustive_try2/run_it005_optimiser.star:_rlnChangesOptimalClasses                                 0.231774
Class3D/first_exhaustive_try2/run_it006_optimiser.star:_rlnChangesOptimalClasses                                 0.234525
Class3D/first_exhaustive_try2/run_it007_optimiser.star:_rlnChangesOptimalClasses                                 0.251032
Class3D/first_exhaustive_try2/run_it008_optimiser.star:_rlnChangesOptimalClasses                                 0.313618
Class3D/first_exhaustive_try2/run_it009_optimiser.star:_rlnChangesOptimalClasses                                 0.325997
Class3D/first_exhaustive_try2/run_it010_optimiser.star:_rlnChangesOptimalClasses                                 0.394085
Class3D/first_exhaustive_try2/run_it011_optimiser.star:_rlnChangesOptimalClasses                                 0.452717
Class3D/first_exhaustive_try2/run_it012_optimiser.star:_rlnChangesOptimalClasses                                 0.442572
Class3D/first_exhaustive_try2/run_it013_optimiser.star:_rlnChangesOptimalClasses                                 0.432944
Class3D/first_exhaustive_try2/run_it014_optimiser.star:_rlnChangesOptimalClasses                                 0.481431
Class3D/first_exhaustive_try2/run_it015_optimiser.star:_rlnChangesOptimalClasses                                 0.391162
Class3D/first_exhaustive_try2/run_it016_optimiser.star:_rlnChangesOptimalClasses                                 0.428301
Class3D/first_exhaustive_try2/run_it017_optimiser.star:_rlnChangesOptimalClasses                                 0.456843
Class3D/first_exhaustive_try2/run_it018_optimiser.star:_rlnChangesOptimalClasses                                 0.404058
Class3D/first_exhaustive_try2/run_it019_optimiser.star:_rlnChangesOptimalClasses                                 0.369498
Class3D/first_exhaustive_try2/run_it020_optimiser.star:_rlnChangesOptimalClasses                                 0.427098
Class3D/first_exhaustive_try2/run_it021_optimiser.star:_rlnChangesOptimalClasses                                 0.565853
Class3D/first_exhaustive_try2/run_it022_optimiser.star:_rlnChangesOptimalClasses                                 0.448246
Class3D/first_exhaustive_try2/run_it023_optimiser.star:_rlnChangesOptimalClasses                                 0.439133
Class3D/first_exhaustive_try2/run_it024_optimiser.star:_rlnChangesOptimalClasses                                 0.445667
Class3D/first_exhaustive_try2/run_it025_optimiser.star:_rlnChangesOptimalClasses                                 0.525791

※ rlnChangesOptimalClassesの意味ですが、relion_refine --print_metadata_labels でメタデータラベルを表示すると、"The number of particles that changed their optimal clsas assignment in the last iteration"と説明されています。"The number"とありますが、上の出力を見ると多分 "fraction"の間違いのような気がします。各イテレーションで所属クラスの割り当てが変化した単粒子の割合を表す数だと思います。計算が収束しているならばこれは 0 に近づいて欲しい気がしますが、おそらく今回は、ほとんどの粒子がclass001とclass002に集中しており、かつclass001とclass002が非常によく似ているため、単粒子がそれらを行ったり来たりし、0.5近傍の値で落ち着いているのではないでしょうか。よくわかりませんが。


3D classificationのSTARファイル

※ こちらで説明します → 3D classificationのSTARファイル

5 High-resolution 3D refinement(高分解能3Dリファインメント)

十分に構造が均質な単粒子サブセットを得られたら、3D auto-refineジョブタイプを使って高分解能化を目指した自動リファインメントを行えます。この処理では、フーリエシェル相関(Fourier Shell Correlation, FSC)を計算する際にいわゆるgold-standard法を使います。gold-standard法では単粒子データセットを半分に分け、それぞれで独立に3D再構成を行い、出来上がった2つの3D構造間でFSCを計算します。これによりオーバーフィッティングを防ぐことができます[Sjors Scheres and Shaoxia Chen, Nature methods(2012)]。投影角度推定の精度を推定する手法[Sjors Scheres, JSB(2012)]と併せて、3D auto-refineはリファインメントが収束したか否かを自動で判定します。ゆえにこの処理はユーザーが入力すべき情報(ユーザーによる主観が入ります)が大変少なくなっています。また、多くのデータセットにおいて素晴らしい密度マップの再構成に成功してきました。この処理はチューニングすべきパラメータがほとんどないので、たいていは1度だけ実行すれば十分というのも利点です。

とはいえ、高分解能リファインメントを始める前に、現在の単粒子セットをより縮小率の小さい状態で切り出しておく必要があります。そのために、まずParticle extractionジョブタイプを実行します。I/Oタブおよびextractタブに以下のようにパラメータを入力します。(※都合によりジョブ番号が異なっています。)

rln-extract4autorefine-1.png

rln-190531-194659.png


3D auto-refineのための単粒子再切り出しパラメータ
タブ パラメータ名 説明
I/O micrograph STAR file CtfFind/job003/micrographs_ctf.star 切り出し元画像とそのCTF情報を記録したSTARファイル
I/O Input coordinates Refine particles STAR fileに指定したSTARファイルの情報を用いて切り出すため、ここは空欄にする。
I/O OR re-extract refined particles? Yes 3D classificationにより投影パラメータをリファインメントされた単粒子を切り出す場合はYesとする
I/O Refine particles STAR file Select/class3d_first_exhaustive/particles.star 3D classification後にSelectジョブタイプで選びだしたクラスに所属する単粒子を記録したSTARファイルを指定する。
I/O OR: re-center refined coordinates? Yes 3D classificationによるリファインメントの結果、単粒子画像の回転中心座標(rlnOriginX, rlnOriginY)は単粒子画像の中心からずれた場所に来ている。このオプションをYesとすることにより、単粒子を切り出したときに回転中心座標が(0, 0)となるよう、単粒子の元画像上での座標(rlnCoordinateX, rlnCoordinateY)を修正(re-center)する。(Recenter on - X, Y, Zを0, 0, 0としたときの挙動。)
I/O Recenter on - X, Y, Z (pix) 0 0 0 上参照。
I/O Manually set pixel size? No Noとするとmicrograph STAR fileに記載のピクセルサイズが使われる。
extract Particle box size (pix) 360 これまでの切り出しでは256 pixelとしていたが、ここではde-localiseされたCTFシグナルをより良くモデル化できるよう、より大きいサイズで切り出す。後々行うCTFリファインメントのために必要。(※蛋白質は結像時にCTF(正確にはその逆フーリエ変換のPSF) で畳み込まれているため、蛋白質の情報は溶媒部分にまで染み出している(de-localise)。それをカットしてしまわないよう、より大きな画像サイズで切り出す、ということではないかと思います。)
extract Rescale particles? Yes 大きなサイズで切り出すとはいえ、あまり大きすぎる画像でチュートリアルを進めるのも大変なので、少しダウンサンプリングする。ダウンサンプリング後の画像サイズは下のRe-scaled size (pixel)で指定し、今回は256 pixelとする。それによりピクセルサイズは 360 [pixel] * 0.885 [Å/pixel] / 256 [pixel] = 1.244 [Å/pixel]となる。するとナイキスト周波数は 2.5 Åとなり、今回使っているごく小さなデータセットでは十分な値である。
extract Re-scaled size (pixels) 256 上参照。


エイリアス名を best3dclass_bigbox とし、【Run!】してください。

上記に加え、ここまでで得られている中で最も良い3Dマップを 256 pixelのボックスサイズとなるようリスケールしておきます。以下コマンドで実行します。

(※私の場合 Class3D/first_exhaustive/run_it_025_class002.mrc が一番いいマップなので、それを使いました。どのクラスが一番いいかは人によって異なります。)

relion_image_handler --i Class3D/first_exhaustive/run_it025_class002.mrc --angpix 3.54 \
     --rescale_angpix 1.244 --o Class3D/first_exhaustive/run_it025_class002_box256.mrc --new_box 256

もともとの3Dマップのピクセルサイズは 3.54 Å/pixelで、それを 1.244 Å/pixel、256 x 256 x 256 [pixel^3]の3Dマップへと変換しています。

5.1 Running the auto-refine job

3D auto-refineジョブを実行しましょう。以下のようにパラメータを入力します。(※都合によりジョブ番号は異なっています。)

rln-autorefine-01.png

rln-autorefine-02.png

rln-autorefine-03.png

rln-autorefine-04.png

rln-autorefine-05.png


ほとんど3D classificationジョブタイプのパラメータと同じです。一部説明を要するところだけ以下に挙げます。

3D auto-refineパラメータ(一部)
タブ パラメータ名 説明
I/O Input images STAR file Extract/best3dclass_bigbox/particles.star 先程リスケールして切り出した単粒子のSTARファイル。※公式のチュートリアルでは Extract/best3dclass_fullsize/particles.star になってますが誤植ですね。
I/O Reference map Class3D/first_exhaustive/run_it025_class002_box256.mrc 先程リスケールした3D参照像(※人によってクラス番号は違うと思います。) このファイルはRELIONのパイプラインにImportしていないためBrowseからは見えない。ファイルパスを手打ちするか、ImportジョブタイプでインポートしてからBrowseで指定する必要がある。
Reference Ref. map is on absolute greyscale? No 今回のプロジェクトで精成した3DマップなのでYesを指定しそうになりそうだが、Noとする。because of the different normalisation of down-scaled images, the rescaled map is no longer on the correct absolute grey scale. (※ ? 縮小率が異なること自体は問題ではなく、ボックスサイズを広く取ることで溶媒領域が広がることが問題、ではないか?(溶媒領域の輝度値をもとに画像正規化をするため) よくわからない。) Noとすると、最初のイテレーションではMAP最適化ではなく、グレースケールの違いに影響されない相互相関計算に基づく方法が用いられ(※? 相互相関係数が最大となる投影パラメータを単粒子に割り当てる、いわゆる古典的な単粒子解析を行うということ?)、正しいグレースケールの3D構造を得られる。それをInitial low-pass filterで指定した値でローパス処理し、それを初期構造として以降のイテレーションでMAP最適化を行う。
Reference Initial low-pass filter (A) 50 高周波数領域にバイアスを持ち込むことを防ぐためと、to maintain the “gold-standard” of completely independent refinements at resolutions higher than the initial one. (※? よくわからない。gold-standardのハーフセット同士をできるだけ独立なものにするため?)
Reference Symmetry D2 高分解能リファインメントを狙っているので対称性を課す。D2対称性を課すことにより、実行的には単粒子の数を4倍したことになる。(※D2対称は非対称単位が4つあり、対称性を課すことによりそれら4つを平均化することになるから。)
Optimisation Use solvent-flattened FSCs? No FSCを計算する際に溶媒平滑化を行うか否か。溶媒平滑化を行うには3Dマスクが必要で、今はまだ無いので、ここではNoとする。)
Auto-sampling Local searches from auto-sampling 1.8 degrees イテレーション後期では局所的な投影パラメータ推定のみが行われる。その際の初期サンプリングレートをここで指定する。イテレーションが進むと自動でサンプリングレートが上がっていく。


CTFタブ、Optimisationタブ、Auto-samplingタブは先に実行した3D classificationジョブとほぼ同じです。Samplingタブにある投影パラメータのサンプリングレート(※Initial angular sampling, Initial offset range, initial offset step)は最初の数イテレーションでしか使われません。それ以降は、リファインメントが収束するまで自動的にサンプリングレートが細かくなっていきます。8面体や20面体対称性よりも対称性が低い場合は(つまりほとんどの場合は)、デフォルトの7.5度のサンプリングレートで十分です。また、イテレーション後期では局所的な投影パラメータ推定が行われます。その時のサンプリングレートの初期値もデフォルトの 1.8度で十分です。対称性の高い粒子のリファインメントを行うときのみ、それぞれ3.7度と0.9度を指定するのが典型的です。

MPIプロセス数について

MPIプロセス数は、まずマスタープロセスのために1つ必要です(スレーブプロセスの管理のみを行う)。gold-standard FSC計算のため2つのハーフセットを独立に再構成しますが、それぞれのハーフセットごとにMPIスレーブが必要です。よって全部で奇数個のMPIプロセスで実行するのが最も効率がよいです。最小MPIプロセス数は 3 になります(マスター1つ、ハーフセット1を計算するスレーブ1つ、ハーフセット2を計算するスレーブ1つ、合計3つ)。

メモリ消費について

一番最後のイテレーションでは、ナイキスト周波数を含む全ての情報を使って最適化計算するため、メモリ消費が著しく増加します。今回使っているデータセットのボックスサイズよりもより大きいボックスサイズで計算する場合は、計算機が持つ全てのコアを使うだけのスレッド数を指定したほうがいいでしょう。


ジョブのエイリアス名を first3dref とし、【Run!】でジョブを実行してください。

5.2 Analysing the results

出力されるファイルは3D classificationジョブと大まかには同じです。しかし3D auto-refineでは、各イテレーションごとにrun_it0??_half?_model.star および run_it0??_half?_class001.mrc というファイルがそれぞれ2つ書き出されます(※ ?は任意の一文字にマッチするワイルドカードです。また、出力されるファイルはこれらだけではありません)。gold-standardのハーフセットごとに1つです。

※一応、イテレーション19の出力ファイルを下に列挙します。half1およびhalf2が、それぞれ各ハーフセット由来のファイルになります。

$ ls run_it019*

run_it019_data.star                    run_it019_half2_class001_angdist.bild
run_it019_half1_class001.mrc           run_it019_half2_model.star
run_it019_half1_class001_angdist.bild  run_it019_optimiser.star
run_it019_half1_model.star             run_it019_sampling.star
run_it019_half2_class001.mrc


イテレーションにより最適化が収束すると、単一の run_model.star ファイル及び run_class001.mrc ファイルが書き出されます(_it0??やhalf?といった文字列はファイル名に含まれていません)。最終イテレーションでは2つのハーフセットから再構成された3D構造がマージされるので、最終イテレーションでは分解能が著しく向上します。また、最終イテレーションではプログラムがナイキスト周波数までのすべての情報を使用するため、メモリとCPUの消費量も増えます。


ひとつ言及しておきますと、投影パラメータのサンプリングレートが自動的に細かくなっていくことは、3D auto-refine処理のポイントの1つです。Sjors Scheres, J. Struct. Biol(2012) のアルゴリズムに従い、各イテレーションでのS/N比に基づいて、投影角度および並進位置ズレ量推定の精度が見積もられます。プログラムは必要以上に細かい投影角度および並進位置ズレのサンプリングは行いません(結果の向上につながらないからです)。推定された投影角度および並進位置ズレの推定の精度、その時用いられたサンプリングレート、その時の分解能の推定結果は、全て _optimiser.star および _model.star ファイルに記録されます。また、RELIONの標準出力にも一部の情報が出力されます。例えば以下を試してみましょう。

$ grep Auto Refine3D/first3dref/run.out


※出力は私の場合以下のようになりました。分解能は 3.9 Åまで行きました。確かに一番最後のイテレーションで分解能が 7.4Å → 3.9Åへ跳ね上がっています。

 Auto-refine: Iteration= 1
 Auto-refine: Resolution= 53.1 (no gain for 0 iter)
 Auto-refine: Changes in angles= 999 degrees; and in offsets= 999 pixels (no gain for 0 iter)
 Auto-refine: Iteration= 2
 Auto-refine: Resolution= 53.1 (no gain for 1 iter)
 Auto-refine: Changes in angles= 19.3751 degrees; and in offsets= 2.64108 pixels (no gain for 0 iter)
 Auto-refine: Estimated accuracy angles= 12.79 degrees; offsets= 3.835 pixels
 Auto-refine: WARNING: The angular accuracy is worse than 10 degrees, so basically you cannot align your particles (yet)!
 Auto-refine: WARNING: You probably need not worry if the accuracy improves during the next few iterations.
 Auto-refine: WARNING: However, if the problem persists it may lead to spurious FSC curves, so be wary of inflated resolution estimates...
 Auto-refine: WARNING: Sometimes it is better to tune resolution yourself by adjusting T in a 3D-classification with a single class.
 Auto-refine: Angular step= 7.5 degrees; local searches= false
 Auto-refine: Offset search range= 5 pixels; offset step= 1 pixels
 Auto-refine: Iteration= 3
 Auto-refine: Resolution= 13.275 (no gain for 0 iter)
 Auto-refine: Changes in angles= 14.6653 degrees; and in offsets= 1.95652 pixels (no gain for 0 iter)
 Auto-refine: Estimated accuracy angles= 6.915 degrees; offsets= 2.49 pixels
 Auto-refine: Angular step= 7.5 degrees; local searches= false
 Auto-refine: Offset search range= 5 pixels; offset step= 1 pixels
 Auto-refine: Iteration= 4
 Auto-refine: Resolution= 9.95625 (no gain for 0 iter)
 Auto-refine: Changes in angles= 17.8327 degrees; and in offsets= 2.04644 pixels (no gain for 1 iter)
 Auto-refine: Estimated accuracy angles= 2.379 degrees; offsets= 0.965 pixels
 Auto-refine: Angular step= 7.5 degrees; local searches= false
 Auto-refine: Offset search range= 5 pixels; offset step= 1 pixels
 Auto-refine: Iteration= 5
 Auto-refine: Resolution= 9.95625 (no gain for 1 iter)
 Auto-refine: Changes in angles= 5.36794 degrees; and in offsets= 0.951966 pixels (no gain for 0 iter)
 Auto-refine: Estimated accuracy angles= 2.021 degrees; offsets= 0.816 pixels
 Auto-refine: Angular step= 7.5 degrees; local searches= false
 Auto-refine: Offset search range= 5 pixels; offset step= 1 pixels
 Auto-refine: Iteration= 6
 Auto-refine: Resolution= 9.95625 (no gain for 2 iter)
 Auto-refine: Changes in angles= 5.20938 degrees; and in offsets= 0.747524 pixels (no gain for 0 iter)
 Auto-refine: Estimated accuracy angles= 2.06 degrees; offsets= 0.839 pixels
 Auto-refine: Angular step= 7.5 degrees; local searches= false
 Auto-refine: Offset search range= 5 pixels; offset step= 1 pixels
 Auto-refine: Iteration= 7
 Auto-refine: Resolution= 9.65454 (no gain for 0 iter)
 Auto-refine: Changes in angles= 5.25729 degrees; and in offsets= 0.735559 pixels (no gain for 1 iter)
 Auto-refine: Estimated accuracy angles= 1.976 degrees; offsets= 0.806 pixels
 Auto-refine: Angular step= 7.5 degrees; local searches= false
 Auto-refine: Offset search range= 5 pixels; offset step= 1 pixels
 Auto-refine: Iteration= 8
 Auto-refine: Resolution= 9.65454 (no gain for 1 iter)
 Auto-refine: Changes in angles= 4.89833 degrees; and in offsets= 0.715826 pixels (no gain for 0 iter)
 Auto-refine: Estimated accuracy angles= 1.885 degrees; offsets= 0.811 pixels
 Auto-refine: Angular step= 7.5 degrees; local searches= false
 Auto-refine: Offset search range= 5 pixels; offset step= 1 pixels
 Auto-refine: Iteration= 9
 Auto-refine: Resolution= 9.65454 (no gain for 2 iter)
 Auto-refine: Changes in angles= 5.29107 degrees; and in offsets= 0.736048 pixels (no gain for 1 iter)
 Auto-refine: Estimated accuracy angles= 2.06 degrees; offsets= 0.835 pixels
 Auto-refine: Angular step= 3.75 degrees; local searches= false
 Auto-refine: Offset search range= 3.68024 pixels; offset step= 0.62625 pixels
 Auto-refine: Iteration= 10
 Auto-refine: Resolution= 8.61081 (no gain for 0 iter)
 Auto-refine: Changes in angles= 3.90416 degrees; and in offsets= 0.581689 pixels (no gain for 0 iter)
 Auto-refine: Estimated accuracy angles= 1.346 degrees; offsets= 0.597 pixels
 Auto-refine: Angular step= 3.75 degrees; local searches= false
 Auto-refine: Offset search range= 3.68024 pixels; offset step= 0.62625 pixels
 Auto-refine: Iteration= 11
 Auto-refine: Resolution= 8.38421 (no gain for 0 iter)
 Auto-refine: Changes in angles= 2.47521 degrees; and in offsets= 0.39645 pixels (no gain for 1 iter)
 Auto-refine: Estimated accuracy angles= 1.195 degrees; offsets= 0.571 pixels
 Auto-refine: Angular step= 3.75 degrees; local searches= false
 Auto-refine: Offset search range= 3.68024 pixels; offset step= 0.62625 pixels
 Auto-refine: Iteration= 12
 Auto-refine: Resolution= 8.38421 (no gain for 1 iter)
 Auto-refine: Changes in angles= 2.30022 degrees; and in offsets= 0.352758 pixels (no gain for 2 iter)
 Auto-refine: Estimated accuracy angles= 1.191 degrees; offsets= 0.554 pixels
 Auto-refine: Angular step= 1.875 degrees; local searches= true
 Auto-refine: Offset search range= 1.76379 pixels; offset step= 0.4155 pixels
 Auto-refine: Iteration= 13
 Auto-refine: Resolution= 7.58571 (no gain for 0 iter)
 Auto-refine: Changes in angles= 1.64242 degrees; and in offsets= 0.281536 pixels (no gain for 0 iter)
 Auto-refine: Estimated accuracy angles= 0.934 degrees; offsets= 0.471 pixels
 Auto-refine: Angular step= 1.875 degrees; local searches= true
 Auto-refine: Offset search range= 1.76379 pixels; offset step= 0.4155 pixels
 Auto-refine: Iteration= 14
 Auto-refine: Resolution= 7.58571 (no gain for 1 iter)
 Auto-refine: Changes in angles= 1.15861 degrees; and in offsets= 0.219828 pixels (no gain for 1 iter)
 Auto-refine: Estimated accuracy angles= 0.877 degrees; offsets= 0.428 pixels
 Auto-refine: Angular step= 0.9375 degrees; local searches= true
 Auto-refine: Offset search range= 1.09914 pixels; offset step= 0.321 pixels
 Auto-refine: Iteration= 15
 Auto-refine: Resolution= 7.58571 (no gain for 1 iter)
 Auto-refine: Changes in angles= 0.758367 degrees; and in offsets= 0.169519 pixels (no gain for 0 iter)
 Auto-refine: Estimated accuracy angles= 0.772 degrees; offsets= 0.4 pixels
 Auto-refine: Angular step= 0.9375 degrees; local searches= true
 Auto-refine: Offset search range= 1.09914 pixels; offset step= 0.321 pixels
 Auto-refine: Iteration= 16
 Auto-refine: Resolution= 7.4093 (no gain for 0 iter)
 Auto-refine: Changes in angles= 0.549974 degrees; and in offsets= 0.135431 pixels (no gain for 1 iter)
 Auto-refine: Estimated accuracy angles= 0.747 degrees; offsets= 0.399 pixels
 Auto-refine: Angular step= 0.9375 degrees; local searches= true
 Auto-refine: Offset search range= 1.09914 pixels; offset step= 0.321 pixels
 Auto-refine: Iteration= 17
 Auto-refine: Resolution= 7.58571 (no gain for 1 iter)
 Auto-refine: Changes in angles= 0.553223 degrees; and in offsets= 0.144504 pixels (no gain for 2 iter)
 Auto-refine: Estimated accuracy angles= 0.759 degrees; offsets= 0.4 pixels
 Auto-refine: Angular step= 0.46875 degrees; local searches= true
 Auto-refine: Offset search range= 0.9 pixels; offset step= 0.3 pixels
 Auto-refine: Iteration= 18
 Auto-refine: Resolution= 7.4093 (no gain for 0 iter)
 Auto-refine: Changes in angles= 0.367747 degrees; and in offsets= 0.140687 pixels (no gain for 0 iter)
 Auto-refine: Estimated accuracy angles= 0.7 degrees; offsets= 0.388 pixels
 Auto-refine: Angular step= 0.46875 degrees; local searches= true
 Auto-refine: Offset search range= 0.9 pixels; offset step= 0.3 pixels
 Auto-refine: Iteration= 19
 Auto-refine: Resolution= 7.4093 (no gain for 1 iter)
 Auto-refine: Changes in angles= 0.288471 degrees; and in offsets= 0.138931 pixels (no gain for 1 iter)
 Auto-refine: Estimated accuracy angles= 0.673 degrees; offsets= 0.385 pixels
 Auto-refine: Angular step= 0.46875 degrees; local searches= true
 Auto-refine: Offset search range= 0.9 pixels; offset step= 0.3 pixels
 Auto-refine: Iteration= 20
 Auto-refine: Resolution= 7.4093 (no gain for 2 iter)
 Auto-refine: Changes in angles= 0.273581 degrees; and in offsets= 0.148914 pixels (no gain for 2 iter)
 Auto-refine: Refinement has converged, entering last iteration where two halves will be combined...
 Auto-refine: The last iteration will use data to Nyquist frequency, which may take more CPU and RAM.
 Auto-refine: Estimated accuracy angles= 0.677 degrees; offsets= 0.389 pixels
 Auto-refine: Angular step= 0.46875 degrees; local searches= true
 Auto-refine: Offset search range= 0.9 pixels; offset step= 0.3 pixels
 Auto-refine: Refinement has converged, stopping now...
 Auto-refine: + Final reconstruction from all particles is saved as: Refine3D/job036/run_class001.mrc
 Auto-refine: + Final model parameters are stored in: Refine3D/job036/run_model.star
 Auto-refine: + Final data parameters are stored in: Refine3D/job036/run_data.star
 Auto-refine: + Final resolution (without masking) is: 3.93333
 Auto-refine: + But you may want to run relion_postprocess to mask the unfil.mrc maps and calculate a higher resolution FSC


※ちなみに run_class001.mrc に原子モデルをはめるとこんな感じです。D2対称性があるとはいえ、5,333粒子からこのクオリティはすごいですね。

rln-autorefine-06.png


※ _angdist.bild というファイルをChimeraなどで開くと、(たぶん)投影角度の分布を示す棒グラフ的なものが見れます。D2対称性があるため球のうちの1/4のみ表示なのではないかと思います(たぶん)。若干てっぺんに分布の密なところがあるので微妙にpreferred orientation(なのか単粒子選別したときにバイアスを入れてしまったのかpickingで拾いにくかったのかわりませんが)が発生しているのが見て取れるような気がします。たぶんtop viewが相対的に多いのだと思います。

rln-autorefine-07.png


3D auto-refineのSTARファイルについてはこちら → 3D auto-refineのSTARファイル

6 Mask creation & Postprocessing (マスク作成 & ポストプロセッシング)

3D auto-refineを行った後は、密度マップのシャープニングを行う必要があります。また、auto-refineの中ではgold-standard FSCの計算時にマスクを施していないマップを使っていました('Use solvent-flattened FSCs'オプションを有効にしていない場合)。それにより、リファインメントの過程で分解能は過小評価されています。溶媒領域のノイズがFSCを下げてしまうためです。RELIONには"post-processing"と呼ばれる処理があり、これはB-factorによるシャープニングと、マスクを施した状態でのFSC計算を行います(Shaoxia Chen et al., Ultramicroscopy(2013))。これを行うために、まずは3Dボックスの中のどこまでが蛋白質でどこからが溶媒領域かを定義するマスクを作成する必要があります。それにはMask Creationジョブタイプを使用します。

6.1 Making a mask (マスクの作成)

Mask Creationジョブタイプを選択し、以下のようにパラメータを入力します。(※都合によりジョブ番号が異なっています。)

rln-mask-01.png

rln-mask-02.png


Mask creationのパラメータ
タブ パラメータ名 説明
I/O Input 3D map Refine3D/first3dref/run_class001.mrc マスクのもととなる密度マップ。
Mask Lowpass filter map (A) 15 滑らかな溶媒マスクを作るのには15Åをカットオフとするローパスフィルタが概ね適している。Input 3D mapにこの値をカットオフとするローパスを適用後、マップが2値化される。
Mask Pixel size (A) 1.244 マップのピクセルサイズ [Å/pixel]を指定する。ローパスフィルタ計算のために使用。
Mask Initial binarisation threshold 0.005 ローパスを施されたマップは、この値を閾値として2値化される。この値は、ローパスフィルタを施したマップを例えばChimeraなどのビューアで表示したとき、蛋白質以外の領域にノイズ由来の密度が一つも生じなくなる等高面の値とすべき。表示等高面の値を上下に振ってみて、ちょうどノイズが見えなくなる値を探すこと。ローパスフィルタを施したマップは、relion_image_handlerプログラムで --lowpass 15 --angpix 1.244 とすれば得られる。大体の場合、0.01 ~ 0.04程度の値。
Mask Extend binary map this many pixels 0 上記までの手続きで2値化されたマスクを得た後、マップの中の白い領域(2値化の結果値1となった領域)を全方位に向かってこの値 [pixel]だけ膨張する。2値化マスクを緩めにするために使用する。
Mask Add a soft-edge of this many pixels 6 (膨張処理の後)マスクのエッジをコサイン型のソフトなエッジにする。ソフトエッジの幅[pixel]をこのオプションで指定する。(Postprocessingにおいて)FSCに対するマスクの影響を計測してFSCを補正する処理は、シャープすぎるマスクに対して極めて敏感なので、エッジをソフトにすることは重要である。マスクの生成は比較的高速な処理なので、いろいろなパラメータでマスクを作成し、一番分解能(FSC)がよくなるところを探すのが良い。


Helixタブは無視します。エイリアス名は first3dref としておきましょう。RELION 3.0以降、この処理はマルチスレッドで高速化することが可能です。【Run!】して実行します。

処理が終了したら、【Display:】ボタンから out: mask.mrc を選択すると、出来上がったソフトマスクのzスライスを表示できます。(以下はσコントラスト0、黒値0、白値1で表示しています。)

rln-mask-03.png


または、UCSF Chimeraで確認することもできます(こちらのほうがおすすめです)。3D auto-refineの結果マップと一緒にマスクを表示し、マスクが3D auto-refineの結果マップを包含していて、かつマスク内に溶媒領域が取り残されすぎていないことを確認しましょう。

$ chimera Refine3D/first3dref/run_half1_class001_unfil.mrc MaskCreate/first3dref/mask.mrc &


以下は、3D auto-refineの最終イテレーションのhalf1の密度マップを灰色で示し、ソフトマスクを黄色のメッシュで示した図です。密度マップの方は、溶媒領域のノイズがギリギリ見えてこない程度の等高面で表示しています。

マスクの表示等高面 = 1

rln-mask-04.png


マスクの表示等高面 = 0.0001

rln-mask-05.png


マスクの等高面1で微妙に蛋白質が削れているところもありますがごく僅かですし、マスクの等高面0では蛋白質領域がすべて包含されているため、これで良い事にします。また、マスクは溶媒領域に密度が観察されませんので、適切な閾値でマスクを作成できたことがわかります。


適切なマスクが得られるまで、新しいパラメータでこのジョブを実行し続けることができます。その場合、古いマスクは mask.mrc.old という名前でバックアップされます。


6.2 Postprocessing (ポストプロセッシング)

Post-processingジョブタイプを選択し、以下のようにパラメータを入力します。(※都合によりジョブ番号が異なっています。)


rln-postprocess-01.png

rln-postprocess-02.png

rln-postprocess-03.png

3D initial modelジョブのパラメータ
タブ パラメータ名 説明
I/O One of the 2 unfiltered half-maps Refine3D/first3dref/run_half1_class001_unfil.mrc ポストプロセッシングしたいイテレーションのハーフマップの片方を指定。(もう片方は勝手にロードされる)
I/O Solvent mask MaskCreate/first3dref/mask.mrc Mask creationで作成したソフトマスク
I/O Calibrated pixel size (A) 1.244 ピクセルサイズ[Å/pixel]を修正したければここで入力。原子モデルを構築し始めたときなどに、ピクセルサイズが実際には数%ずれていることに気づいた、といったことはしばしば起こりうる。RELIONの中では、ここまでの処理はすべて辻褄があっているため、マップを正しいピクセルサイズでリファインメントし直したりクラス分類をやり直したりする必要は無い。ただここで正しいピクセルサイズを入力し、それをもとにマップの修正と最終的な分解能推定を行えば良い。今回は特にずれは存在しないため、現在のマップのピクセルサイズをそのまま入力する。
Sharpen MTF of the detector (STAR file) mtf_k2_300kV.star 使用した検出器のModulation Transfer Function(MTF)を記述したファイルを与える。Falcon2@300kV, DE-20@300kV, K2-summit@300kVのMTFならばRelionのウェブサイトにSTAR形式でアップされている。今回はK2-summitのデータのため、K2-summitのMTFを mtf_k2_300kV.star という名前のファイルにコピペして保存し、それをここで指定すれば良い。もしMTFが不明なら空欄でよく、その場合はB-factor推定で検出器のMTFごと推定、補正される。多少不正確になるが、大きくマップの質を変えてしまうことは無いだろう。(※試してみました→MTF無しでPostProcessing)
Sharpen Estimate B-factor automatically? Yes この処理ではRosenthal and Henderson, JMB(2003)の方法に基づいてB-factorを自動推定する。そのため分解能が10Åを大きく超えている必要がある。もし10Å分解能に届いていない場合は推定はうまくいかないのでここはNoとし、自分でB-factorの値をアドホックに入力する必要がある。
Sharpen Lowest resolution for auto-B fit (A) 10 普通は変えなくて良い。B-factor推定時、ギニエプロットの線形フィッティングに含める最小の分解能(Å)を指定する。10Åからあまり離れた値は使うべきではない。
Sharpen Use your own B-factor? No 今回は自動推定するのでここはNo
Filter Skip FSC-weighting? No Noとすると、マスクの影響を補正したgold-standard FSC曲線に従い、3Dマップにローパスフィルタを施す。FSCで計測されるのは(ローカルマスクを施したりしない限り)マップの全体的な分解能なので、もしマップ全体の分解能よりも局所的に良い分解能を持つ領域を分析、可視化したい場合はこれをYesとし、Ad-hoc low-pass filterで好きなローパスフィルタを指定し施すこともできる。


エイリアス名は first3dref とし、【Run!】してください。計算はすぐ終わるので性能の良い計算機を使う必要はありません。

計算が終わったら【Display:】でポストプロセッシング後のマップ(マスクなし: postprocess.mrc、マスクあり: postprocess_masked.mrc)のZスライス、FSC曲線とギニエプロットが描かれたPDFファイルを閲覧できます。また、PostProcess/first3dref/postprocess.mrcを Chimeraで表示すると、3D auto-refineで収束した3Dマップよりもαヘリックスが見やすくなっていることがわかると思います。

※ポストプロセッシング後のマップに原子モデルをフィッティングした様子を以下に示します。

rln-postprocess-04.png


分解能の計測はShaoxia Chen et al., Ultramicroscopy(2013)で報告した位相ランダム化処理に基づいています。【Display:】で out: logfile.pdf を開くと、1ページ目に4つのFSC曲線が描かれています。

  • 黒がマスクの影響を補正後のFSC(rlnFourierShellCorrelationCorrected)
    • ※分解能はこのFSCカーブで評価します。
  • 緑がマスクをかけていない状態のFSC(rlnFourierShellCorrelationUnmaskedMaps)
    • ※溶媒領域のノイズの影響で過小評価されているFSCです。
  • 青がマスクをかけた状態のFSC(rlnFourierShellCorrelationMaskedMaps)
    • ※マスクにより溶媒領域のノイズの影響は軽減されていますが、2つのハーフマップに同一のマスクを施したことにより偽相関が生じ、FSCが過大評価されます。マスクのエッジがシャープすぎたり、マスクが細かい情報を持ちすぎたりしていると、より強い偽相関が生じます。
  • 赤が位相ランダム化を施したマップのFSC(rlnCorrectedFourierShellCorrelationPhaseRandomizedMaskedMaps)
    • ※ポストプロセッシングジョブの標準出力ファイルに "+ randomize phases beyond: 9.65042 Angstroms"などという記述が見られると思います。赤のFSCカーブは、それよりも高分解能の領域の位相をランダムな値で置き換えた後でマスクを施したマップ間で算出したFSCです。位相ランダム化された分解能では本来相関が0に近くなるべきですが、マスクで偽相関が発生すると位相ランダム化された領域でもFSCが正の値を持ちます。この位相ランダム化FSCで青のFSCカーブを補正することにより、マスクによる偽相関を差っ引いた真のFSC(黒のカーブ)を得られます(Shaoxia Chen et al., Ultramicroscopy (2013))。

です。(これらは PostProcess/first3dref/postprocess.star の data_fscデータブロックにカッコ内をラベル名として記録されています。)

rln-postprocess-05.png


※マスクで溶媒ノイズ領域を削ってFSCを計測したことにより、分解能が 3.9 Å → 3.1 Åに向上しています。(向上したというか、もともと蛋白質領域近傍の分解能は3.1Åで、その領域の分解能をより正確に測れたと言う方が正確でしょうか)


位相ランダム化したマップのFSC(赤)が、ポストプロセッシングで推定された分解能の近辺で 0 に近い値をとっているか確認してください。もしそうでない場合、マスクのエッジがシャープすぎたり、細かい情報を含みすぎです。その場合はMask creationステップに戻り、より強いローパスフィルタをかけたり、より広くソフトなエッジにしたマスクを作り直し、再度ポストプロセッシングを行いましょう。

Actually, you may want to make the mask you just made a bit softer, so perhaps try an edge of 5 pixels as well? (※???先程6 pixelのコサインエッジを使ったのに、5 pixelでsofterとは??)

Mask creationで異なるパラメータでマスクを再度作りたい場合、Finished jobsからすでに終了したジョブを選択してパラメータを変更して【Coutinue!】すると前回の結果が上書きされます。またはジョブタイプブラウザからMask creationジョブを新規に立ち上げて新しいパラメータで【Run!】すると、前回の結果を上書きすることなく新しいマスクを作成できます。お好みで。

7 CTF and beamtilt refinement(CTFおよびビーム傾斜のリファインメント)

RELION 3.0より、単粒子画像ごとのデフォーカス量の推定と、データセット全体のビーム傾斜の推定を行うプログラムが追加されました。これはCTF refinementジョブタイプとして実装されています。この処理により、比較的少ない計算コストでさらなる分解能向上が見込めます。この処理は前の3D auto-refineジョブおよびそのPost-processingジョブの結果を出発点として走らせます。

7.1 Running the job

CTF refinementジョブタイプを開き、以下のようにパラメータを入力します。(※都合によりジョブ番号が異なっています。)

rln-ctfrefine-01.png

rln-ctfrefine-02.png


CTF refinementジョブのパラメータ
タブ パラメータ名 説明
I/O Particles (from Refine3D) Refine3D/first3dref/run_data.star 3D refinementにより最適化された単粒子のSTARファイルを指定。
I/O Postprocess STAR file PostProcess/first3dref/postprocess.star PostProcessジョブにより生成されたSTARファイル。そのジョブにおいて使用されたマスクが3D refinementのハーフセットの3D再構成マップに施され、FSCが計算される。CTF refinement時はそのFSCに基づいてフーリエ空間が重み付けされる。
Fit Minimum resolution for fits (A) 30 デフォルト値のままにしておく。ビーム傾斜推定時に使われる最小分解能[Å]を指定。
Fit Perform CTF parameter fitting? Yes 今回は単粒子毎にCTFパラメータを最適化したい。それにより、氷の層が水平でないことや、単粒子が氷の層内の異なる高さに分布することによるデフォーカスずれを補正したい。
Fit Range for defocus fit (A) 2000 普通はこのデフォルト値で上手くいく。入力STARファイルのデフォルト値を初期値とし、±この値[Å]の範囲の中でフィッティングを行う。
Fit Fit per-micrograph astigmatism? No Yesとすると元画像単位で非点収差のリファインメントを行うが、それには大量の単粒子画像と高いS/N比が必要。今回はNo。
Fit Fit per-particle astigmatism? No Yesとすると単粒子画像単位で非点収差のリファインメントを行うが、非常に大きい単粒子でかつS/Nが素晴らしく良い場合など、とても強いデータでしかうまくいかない。今回はNo。
Fit Fit per-particle phase-shift? No Yesとすると元画像単位でCTFの位相シフト(振幅コントラスト由来)のリファインメントを行う。ボルタ位相版を使用したデータで役に立つかもしれないが、元画像あたりの単粒子数が多く、かつS/Nが良いデータセットでないと難しい。今回はNo。
Fit Perform beamtilt estimation? Yes データセット全体に対し、X方向およびY方向のビーム傾斜量[mrad]を推定する。


※ビームシフト/イメージシフトによりグリッドを動かさず多数の穴を取るなどした場合、あるいはデータ収集中に照射系のアラインメントをやり直した場合などは、特定の画像グループ毎に個別にビーム傾斜推定したいかと思います。_rlnBeamTiltClassなるラベルをSTARファイルに追加してあげるとグループ単位で推定してくれるみたいです、どうやら → https://bitbucket.org/scheres/relion-3.0_beta/issues/3/beam-tilt-estimation-group-micrographs (これはβ版時代のIssueらしく、今のバージョンでどうなってるかまでは確認していません。)


このプログラムはCPUでのみ動作します。(※本家の環境では)4 MPI、各8スレッドの実行でだいたい1分ほどです。

7.2 Analysing the results

【Display:】から out:logfile.pdf を選択すると、各単粒子のリファインメント後のデフォーカス値を確認できます。氷の層が水平でない画像が見つかりますか?

※以下は結果の一部です。横軸および縦軸はそれぞれ元画像のx座標、y座標[pixel]を表していて、その画像内の各単粒子のデフォーカス量が色でプロットされています(青が最小デフォーカス、オレンジが最大デフォーカス?)。この画像は(たぶん)デフォーカス値の最大と最小の差が 188 Å ということなので、氷の層がごく薄いか、もしくは単粒子が氷の層の中の(z方向に)狭い範囲に分布していることがわかります。また、わずかですがxが大きくなるに連れてデフォーカスが大きくなる傾向がありますので、わずかに氷の層が傾いているかもしれません。

rln-ctfrefine-03.png


ビーム傾斜量の推定値は CtfRefine/job024/beamtilt_0.txt に記録されています。(※beamtilt_0.txtと書いてありますが beamtilt_iter-fit_class_0.txtおよびbeamtilt_lin-fit_class_0.txtのことでしょうか?あるいは本家チュートリアルはRELION 3.0のものなので、古いバージョンの挙動でしょうか。)

※ iter-fitとlin-fitはどちらもビーム傾斜量の推定値が記録されていますが、iter-fitの方が最終的な推定結果だと思われます。ソースコードを見るとまずlin-fitが推定され、それを初期値とした最適化計算によりiter-fitが推定されているようにみえ、かつCTF refinementジョブで生成されるparticle starファイル particles_ctf_refine.starの rlnBeamTiltX, rlnBeamTiltYはiter-fitに記載の推定値とおなじになっています。

※ beamtilt_iter-fit_class_0.txt の中身は、以下のようになっています。ビーム傾斜の推定値はX方向に 0.02 mrad, Y方向に -0.2 mrad です。

$ cat CtfRefine/job049/beamtilt_iter-fit_class_0.txt

beamtilt_x = 0.0236249
beamtilt_y = -0.214688


また、このビーム傾斜推定に使用された位相ずれを表した画像は以下のコマンドで確認できます。

(※CTF補正済みの参照像の位相と比較したときのズレと思います。ビーム傾斜はコマ収差による位相ズレをモデルでフィッティングすることにより推定されています。位相ずれが最も大きくなる方向がビーム傾斜の方向です。)

$ relion_display --i CtfRefine/job049/beamtilt_delta-phase_per-pixel_class_0.mrc --scale 3 &

rln-ctfrefine-04.png


※若干ですが、Y軸方向の下側ほど白く(位相ずれが正に大きい?)、Y軸方向の上側ほど黒く(位相ずれが負に大きい?)なっています。X方向にはそのような傾向は見られません。ビーム傾斜の推定値と対応しています。


フィッティングされたビーム傾斜(位相ずれ)モデルは以下のコマンドで確認できます。(※本家チュートリアルでは lin-fit を表示していますが、最終的なフィッティング結果は iter-fitの方な気がするので、ここではiter-fitを示しておきます。)

$ relion_display --i CtfRefine/job049/beamtilt_delta-phase_iter-fit_class_0.mrc --scale 3 &

rln-ctfrefine-05.png


もしビーム傾斜が大きい場合、これらの画像は円の片側が他よりも暗く、反対側が他よりも明るくなります。今回の場合は推定されたビーム傾斜量は X方向に 0 mrad, Y方向に -0.2 mrad でした。今回のデータの分解能では、これはあまり問題にならないくらい小さな値です。ビーム傾斜の推定無しで Ctf refinementをやり直してもいいでしょう。ここでは特にそうしたことはせず、このまま先へ進むことにします。

8 Bayesian polishing(ベイジアンポリッシング)

3.0がリリースされた現在、ビームによって引き起こされた運動補正をベイジアン法をRELIONも実行します。この方法は規則正しくされた見込みの最適化が狙いです。強いることなしに空間的干渉性の行為と時間的に滑らかな動きが限定する前の見込みで粒子軌道のそれぞれの仮説のセットを連想することを許します。前の期間の平滑度は動きを観察した統計学で説明する3つのパラメータを必要とします。この粒子のデータセットのための最も良い収量の動きを追跡する前の見積もりをすることは最初のtraining modeのプログラムを実行することができます。

8.1 Running in training made(トレーニングモードの実行)

平行な16の筋道を使うことでこのjobはコンピュータで1時間と15分かかります。もしこれを待ちたくないのであれば8.2へ進んであらかじめ計算された結果からsigma-valuesを使うことができます。もし、このjobを実行したいのであれば、job-typeをBaysian polishingを選択して以下のように設定してください。

●Micrographs (from MotinCorr): MotinCorr/own/corrected_micrographs.star

(このMotioncorrection jobはRELION3.0かそれ以上のもので実行されていることは重要です。運動補正について必要なメタデータが書かれていないのでRELION 2.1 もしくはそれ以下ではMotion correction jobは動きません。)

●Perticles (from Refine3D or CtfRefine):CtfRefine/job024/particles_ctf_refine.star

(これらの粒子はポリッシュされます。)

●Postprocess STAR file Postprocess/first3dref/postprocess.star

(このjobからのマスクとFSCカーブはポリッシュの手順で使われます。)

●First movie frame:1
●Last movie frame:-1

(それらの映画から最初もしくは最後のフレームを捨てる人もいます。RELIONでベイジアンポリッシングを行うときこれは勧めないことを注意してください。B-ファクターの映画フレームの重さは光る粒子のシグナルとノイズの割合で自動的に最適化するので、全ての映画フレームを含むための最善です。) [Train]タブにセットして以下のように設定してください。

●Train optional parameters? Yes
●Fraction of Fourier pixcels for testing:0.5

(デフォルトのままにしてください。)

●Use this many particles:5000

(それはほとんどすべてのものです。さらに多くの粒子、多くのRAMでこのプログラムを行うことを留意してください。もし、メモリーが尽きたら、少量の粒子でトレーニングをしてみてください。5000よりもとても少ない粒子を使うのにおすすめしないです。) Polishタブを確かめてセット以下のように設定してください。

●Perform particle polishing? No

このプログラムのトレーニングステップは、MPIと並行しないことに留意してください。それゆえに、1つのMPIの過程のみを使う確認をしてください。プログラムのスピードを上げるため16の道筋でプログラムを実行しました。まだ、計算には1時間以上かかります。trainのエイリアスを使いました。

8.2 Running in Polishing mode(ポリッシングモードの実行)

一度トレーニングステップが終わってしまえば、そのプログラムはPolish/train/opt-params.txtと呼ばれるテキストファイルに書き出されます。粒子をポリッシュするためにこれらのパラメータを使うにはjob-typeメニューの左にある新しいBayesian polishingを選んでクリックしてください。以前と同じ[I/O]タブのパラメータを維持して[Train]タブをトレーニングオフに変換したことを確認してください。その時、[Polish]タブをセットして以下のように設定してください。

●Perform particle polishing? Yes
●Optimised parameter file: Polish/train/opt-params.txt
●OR use your own parameters? No
●Minimum resolution for B-factor fit (A):20
●Maximum resolution for B-factor fit (A):-1

(最後の2つのパラメータはデフォルトのままにしてください。) 代わりにもしトレーニングセットをスキップすると決めたならば、実行することで得たsigma-parametersでPolishタブを満たすことができます。

●Perform particle polishing? Yes
●Optimised parameter file: 

(下記に従って手に入れた最適なパラメータを使うのでここは空のままにしてください。)

●OR use  your own parameters? Yes
●Sigma for velocity (A/dose) 0.462
●Sigma for divergence (A) 1935
●Sigma for accleration (A/dose) 2.46
●Minimum resolution for B-factor fit (A):20
●Maximum resolution for B-factor fit (A):-1

(最後の2つのパラメータはデフォルトのままにしてください。) このプログラム部分はMPI並行化されています。3MPIを使う過程でそれぞれ16道筋、2分で実行し終えます。ポリッシュのエイリアスを使いました。

8.3 Analysing the results (結果の分析)

Bayesian polishing jobは、shuny starと呼ばれる粒子をポリッシュしたSTARfileとlogfileのPDFに出力します。後者はスケールのプロットを含み、放射線障害の重さにため、B-ファクターが使われました。加えて、粒子の精密化のプロットは、全ての顕微鏡写真のすべて含んだ粒子を追跡しています。このデータセットのプロットを見てください。そのステージの流れが少し出現します。ほとんどすべての粒子は映画の間、上の右から左の下へ動きます。ポリッシュ後、粒子のシグナルとノイズの割合は改良されていました。そして、それは新しい3D auto-refine jobと一致するPost-processing jobを提示すべきです。[I/O]タブで下記のオプションを使うため、光る粒子で3Dauto-refineを実行してください。

●Reference mask (optinal): MaskCreate/first3dref/mask.mrc

(これは最初のPost-processing jobのために作ったマスクです。このオプションを使うことで、マスクの外側のすべてのピクセルを0にセットします。これは参考の中でノイズが減少し、このようにさらに良い配向指定を導き、このように再構築します。) [Optmisation]タブ上で下記のように設定してください。

●Use solvent-flattened FSCs? Yes

(精密化は毎相互作用でFSCカーブの基準で溶媒補正を使います。このオプションを良く使うことは、Post-processing job-typeで使われるようです。粒子の箱の中身の相対的に小さい体積をタンパク質が占めるときこのオプションは特に便利です。例えば、分子をとても引き伸ばすかマスクを使って小さい部分で精密化して焦点を合わせるときです。3D精密化でFSCを計算することのデフォルトの方法は精密化の間、分解能は下回っているマスキングしているハーフマップ以外である。これは、毎相互作用でFSCカーブの溶媒補正の位相ランダム化を計算することによって正しく設定され、これは通常注目すべき分解能の改善を導きます。) 以前計算した結果を見ることができるので最終的に2.9Åに伸びている分解能を得ることができます。3GBのデータのために悪くないということですか?

8.4 When and how to run CTF refinement and Bayesian polishing(いつどのようにしてDTF精密化とベイジアンポリッシングを実行するか)

粒子ごとに焦点をぼかしてビームチルトを見積もったベイジアンポリッシングとCTF精密化はともに分解能の再構築を改善するかもしれません。これは最初にどちらを適用すべきなのかという疑問が生じます。この例では、粒子ごとに焦点をずらした値とビームチルトの精密化を最初にしたが、最初にベイジアンポリッシングも行うことができます。両方法とも高い分解能モデルから恩恵を被るので相互作用の手順が有益であるかもしれません。例えば、ベイジアンポリッシングの後、CTFの精密化を繰り返すことができます。通常、初めに最も大きい問題といくつかの試みに取り組むことがおそらく最善です。そして、エラーは場合によっては、必要です。さらにトレーニングが矛盾した結果を生成するケースがいくつか見られます。すなわち、とても異なるシグマ値の収量で複数実行します。しかしながら、よくポリッシングのために使われた実際のシグマ値が光っている粒子を再精密化した後、分解能のマップのためにあまり重要ではないと述べました。それゆえに、トレーニングは計算的に高いのでデフォルトのパラメータで直接ポリッシングを実行した方がいいでしょう。(σvel=0.2;σdiv=5000;σacc=2)すなわち特定のデータセットでは、トレーニングなしがいいでしょう。

9 Locak-resolution estimation(ローカル解像度推定)

後処理プログラム[Post-processing]から推定される分解能は広域的な推定となります。しかしながら、たった1つの数では高分子複合体の再構築でしばしば観測される分解能が場所ごとで異なることを描写できません。Alp Kucukelbir and Hemant Tagare はマップ全体の解像度の変化を推定する効果的なプログラムを記述しています。[7] RELIONは[Local resolution]job-type を通してプログラムのラッパーを実装しています。代わりに全体のマップ周りを動く柔らかい球状のマスクでpost-processingのような手順を実行するために選ぶことができます。以下が例です。

9.1 Running the job(jobの実行)

[I/O]タブ上で以下のように設定してください。

●One of the two unifiltered half-maps:Refine3d/polished/run_half1_class001_unfil.mrc
●User-provided solvent mask:MaskCreate/first3dref/mask.mrc
●Calibrated pixel size: 1.244

(時には、正しいピクセルサイズと思ったモデルを作り始めたときに実際は何%かずれていることを発見します。RELIONの中身はこの点までのすべては、首尾一貫していました。なので、マップの精密化や/とデータの分類化は必要がありませんでした。する必要があるすべては、正しい地図と最終的な解像度推定のためにここに正しいピクセルサイズを提供することです。) ResMapをセットして、Use ResMap?をNoとしてRelionタブをセットして以下のように設定してください。

●Use Relion? Yes
●User-provided B-factor:-30

(この値は、ローカリーフィルターとマップを鋭くする計算をすることも慣れています。おそらく、Post-processing jobの間で機械的に測定した値に近い値を使いたいでしょう。)

●MTF of the detector (STAR file):mtf_k2_300kV.star

(Post-processing jobと同じです。)

9.2 Analysing the results(結果の分析)

エイリアスとしてポリッシュを使いました。8MPIprocessで実行したときこのjobは7分かかります。出力はLocalRes/polished/relion_locres.mrcと呼ばれます。ローカル解像度によりPostprocess/polished/postprocess.mrcマップに色を付けるUCSFChimeraで使われたかもしれないファイルです。Tools->Volume data->Surface colorを使ってby volume data valueを選んでresmapファイルを拾い読みます。特有のRELIONのオプションは1つのマップの質における全体の差異をするlocally-filteredとsharpend mapの出力を追加します。このマップは、local/Res/polished/relion_locres_filtered.mrcに保存されUCSFChimeraで直接ビジュアライズすることができます。(オプションで前のローカル解像度によって色を付けます。)

10 Checking the handedness(利き手のチェック)

マップの慎重な検査は利き手が間違ったことを指し示します。例えば、α-ヘリックスは間違った方法の向きを変えるからです。顕微鏡使用ステージのチルティングなしでデータセットから絶対的な利き手を決心するのは不可能であることを覚えておいてください。したがって、3D initial model job-typeのSGPアルゴリズムは、逆の手である50%の可能性があります。あらかじめ計算された結果では、これは本当に本当です。ある方法で以下の通りに後処理されたマップの利き手を動かしてください。 relion_image_handler --i PostProcess/polished/postprocess.mrc --o PostProcess/polished/postprocess_invert.mrc --invert_hand 同じコマンドもまた他のどのマップで実行することができます。手が間違っていると画像処理手順の前に気が付いたら他の手に変更してください。RELIONは自身で両手が区別されない問題を解決することができませんが、気が付くとすぐに手を動かすことはより便利である場合があります。一度正しい手にするとマップでUCSFChimeraに読み込みβ-ガラクトシダーゼの原子モデルで上書きしたいとします。PDBID 5a1aを用いたPDBから直鎖上の原子モデルをフェッチすることができたはずです。

11 まとめ

11.1 フローチャートの作成

どのようにして最終的な再構築をするか気になるでしょう。Finished jobリストから実行した最後の動作を選択し、Job actionsボタンからMake flowchartオプションを試してください。これを実行するためにはシステムにLATEXTikZパッケージが必要となります。最初の項は正確な作業名のないフローチャートの全体像です。これは、出版目的に役立つでしょう(おそらくお気に入りのvector-based designプログラムの編集後に)。フローチャートの全体像の後に、最初の詳細なフローチャートにどのようにして終了するかが示されています。10ステップ以上からなるフローチャートは複数の構成要素に分けられることに注意してください。作業中に複数枝分かれすることがあるでしょう。それ故に、最終作業のフローチャートに従って、それぞれの分岐についてフローチャートがあるでしょう。リンクをクリックすることでPDFファイルから合致する位置を得ることができます。

11.2 ディレクトリのクリーンアップ

ディスクの空間を確保するために、RELIONは作業ディレクトリをクリーンアップするオプションを持っています。これには2つのモードがあります。gentleクリーンは作業ディレクトリから中間ファイルのみを削除します。harshクリーンは作業からのインプットを必要とする新たな作業を開始するのに必要なファイルを削除します。例えば、harshクリーンはMotion Corr作業から平均化された、あるいはParticle extraction作業から少量のスタックを抽出された顕微鏡写真を削除します。一方、gentleクリーンは2D classification, 3D classificationあるいは3D auto-refine作業の中間反復からすべてのファイルを削除します。Job actionsボタンから個別に作業をクリーンアップすることができます。あるいは、GUIのトップメニューから下がった、Jobsからすべての作業をクリーンアップすることができます。私たちは計算前の結果に頭を悩ませていたプロジェクトディレクトリのtarballを作る前にそのメニューからGently clean all jobsオプションを使いました。長期間保存されていたデータを取る前にプロジェクトディレクトリをクリーンアップしたいでしょう。

11.3 質問や引用

以上です!このチュートリアルを楽しんで、そして役に立つことを願います。もしRELIONについて質問があれば、初めにRELION WikiFAQCCPEMメーリングリストを確認してください。それが助けにならなければ、CCPEMリストを利用して質問してください。お願いですから、決してSjorsに直接メールを送信しないでください。質問全てに返答できるわけではありません。もしRELIONが利用者の研究に有益であると思ったら私たちの論文を引用し、あなたの研究仲間たちに伝えてください。

11.4 参考文献

RELIONの手順の改善を支持した理論の詳細は以下の論文で述べられています。

  • S.H.W. Scheres (2012) “RELION: Implementation of a Bayesian approach to cryo-EM structure determination” J. Struc. Biol., 180, 519-530.
  • S.H.W. Scheres (2012) “A Bayesian view on cryo-EM structure determination” J. Mol. Biol., 415, 406-418.

RELIONの使い方に関する包括的な概観は以下の論文に述べられています。

  • S.H.W. Scheres (2016) “Processing of structurally heterogeneous cryo-EM data in RELION” Meth. Enzym., 579, 125-157.


12 付録 A:インストールの注意

12.1 MPIのインストール

MPI(message passing interface)がインストールされた計算機クラスター、もしくはNVIDIA GPU搭載のマルチコアデスクトップマシンが必要になることに注意してください。RELIONをコンパイルするには、mpi-develパッケージが必要です。おそらく存在する種類(MPICH, LAM-MPI, etc)またはバージョンはそれほど重要ではありません。まだシステム内にmpi-develをインストールしてなければ、openMPIをインストールすることをお勧めします。

12.2 CUDAのインストール

NVIDIAの比較的新しいGPU(compute capability 3.5+以上)を持っている場合、オートピッキング、分類、精密化作業を高速に実行できます。GPU-accelerationサポートを含むRELIONをコンパイルするには、CUDAをインストールする必要があります。このチュートリアルの準備にはCUDA-8.0を使用しました。NVIDIAのWebサイトからダウンロードしてください。

12.3 RELIONのインストール

RELIONはオープンソースのソフトウェアです。the RELION wikiから無料でダウンロードして、手順に従ってインストールしてください。もし、job submission system(Sun Grid EngineやPBS/TORQUE etc)に詳しくない場合は、インストール手順の説明にあるように、qsub.cshスクリプトの設定に関して、システム管理者にお尋ねください。分散メモリ並列化のためのMPIと、共用メモリの並列化のためのpthreadsの両方を使用する、いわゆるハイブリッド並列計算を実行したいときは注意してください。ジョブキューイングシステムはこれを可能とするためにいくつかの調整が必要な場合があります。再度、システム管理者にお尋ねください。

12.4 モーション補正ソフトウェアのインストール

RELION-3.0は全フレーム顕微鏡写真ムービーアライメントに使用される、UCSFプログラムMOTIONCOR2へのラッパーを提供します。Dabid Agardのページからプログラムをダウンロードし、インストールの手順に従ってください。また、RELIONの保有する、MOTIONCOR2の実行を使うかもしれません。なので、もし、UCSF実行をインストールするのに問題が発生しても心配しないでください。バージョン3.0において、Nico grigorieff’s groupからUNBLURへのラッパーはGUIから中止されていることに注意してください。

12.5 CTF推定ソフトウェアのインストール

CTF推定はRELIONに含まれていません。代わりに、RELIONはAlexis RohouやNiko Grigorieff’s CTFFIND4へのラッパーを提供します。Niko’s CTFFIND websiteからダウンロードし、手順に従ってインストールしてください。また、もし、パソコンにNVIDIA graphics card(GPU)が入っているならKai Zhang’s GCTFを使うこともできます。その際は、LMBに関するKaiのWebサイトからダウンロードしてください。

12.6 RESMAPのインストール

局所分解能推定はRELIONが所有するポストプロセッシングプログラムの中、あるいはAlp KucukelbirのRESMAPへのラッパーを通して実行される。AlpのRESMAP Webサイトからダウンロードし、手順に従ってインストールしてください。



13 付録B:RELIONの利用

13.1 GUIについて

13.1.1 pipeline approach

RELION-3.0では、一つのプログラムの出力を次のプログラムの入力として順次繋いでいくような情報の受け渡し(パイプライン)を円滑に実現できるようなGUIを設計しました。詳細は2016 Proceedings of the CCP-EM Spring Symposium[3]で公開されています。  RELION-3.0を使うにあたって、まずプロジェクト(すなわち構造決定したい分子)ごとにディレクトリを作成することを推奨します。ここではこのディレクトリの呼称をプロジェクトディレクトリとします。RELIONのGUIは「relion」とコマンドを打つことで起動させることができますが、この動作は常にプロジェクトディレクトリに移動してから行うことを念頭に置いてください。

このソフトのGUIでは、すべてのジョブや、あるジョブによる出力がどのように他のジョブの入力に使われたのかといった履歴を保存することができ、この機能によってパイプラインが形成されます。各種ジョブにはそれぞれ独自の出力ディレクトリが与えられます。ここではたとえば あるジョブの出力のためにClass2D/ というディレクトリができるとします。こういったジョブディレクトリの中でさらに新しいジョブが実行されると、その新しいジョブには一連の番号が与えられます。たとえばClass2Dのもとで作られたとあるジョブにはClass2D/job010という番号が振り分けられます。さらに、ジョブディレクトリのもとでは、Class2D/job010/run といったように出力名が定められます。このソフトでは、各ジョブに意味のある名前が付けられるような仕組みを備えるために、ファイルシステム上の個々のジョブディレクトリへの記号リンクとして実装された"エイリアス"システムを使用しています。 パイプラインに関するすべての情報はdefault_pipeline.starという名前のファイルに保存されていますが、基本的にはこのファイルを見る必要はありません。もしdefault_pipeline.starファイルが壊れてしまった場合は、最後に実行したジョブディレクトリで保存されたバックアップから復元することができます。

13.1.2 The upper half: jobtype-browser and parameter-panel(上半分:jobタイプブラウザとパラメータパネル)

このソフトのGUIは上半分と下半分に分かれています。上半分の左側にはジョブの種類が縦リスト上に並べられており、行いたいジョブを選ぶことができます。たとえば2D classificationを行いたいときは、上半分の左側のブラウザで選択すればいいわけです。上半分の右側には複数のタブがあるパネルがあり、それぞれのジョブタイプへのパラメータを入力できます。GUIの左上には機能のおおよその概要が分かる3つのメニューが書かれています。 【Schedule】ボタンでは、後でジョブを実行するときのスケジューリングを行うことができます。それに対して【Run now!】ボタンは今すぐに実行したいときに使います。【Schedule】ボタンは、繰り返し実行できるよう完全自動化された「パイプライン」を用意したいときに役立ちます。たとえば、データが収集されていくにつれてそれが瞬時に反映できるようなパイプラインを実現したいときに有効です。詳細は13.3を参照してください。 GUIの左側のjobtype-browserをクリックすると、新しい動作(【Run now!】ボタン付き)が右側のパラメータパネルにロードされます。

13.1.3 The lower half: job-lists and stdout/stderr windows

GUIの下半分には、まだ実行中のジョブ(【Running jobs】)、すでに完了しているジョブ(【Finished jobs】)、または後で実行するようにスケジュールされたジョブ(【Scheduled jobs】)のリストがあります。これらのリスト内のジョブをクリックすると、そのジョブのパラメータがパラメータパネルにロードされ、【Run now!】ボタンの色が変わり、【continue now!】に変わります。それをクリックすると、新しい出力jobディレクトリが作成されることはありませんが、パラメータパネルで指定されたパラメータに従ってjobが続行されます。【2D classification】、【3D classification】、【3D auto-refine】のジョブでは、_optimiser.starというファイルが必要になり、ファイル名には、run_ct23のような継続された繰り返しのファイル名が付きます。他のjobタイプでは、それらが以前に実行されるまで、その時点から継続することができます。【モーション補正】、【CTF推定】、【自動ピッキング】、【粒子抽出】は、これまでに行われていなかった顕微鏡写真でのみ実行されます。【Input to this job】(このjobへの入力)と【Output from this job】(このjobからの出力)は、リンクjobをまとめてリストし、プロジェクト履歴内を前後にブラウズするために使用できます。

GUIの下半分の下部には、選択された(実行済みまたは実行中の)jobの標準出力(stdout)と標準エラー出力(stderr)がそれぞれ黒と赤のテキストで表示されます。標準エラー出力は理想的には何もない状態でなければなりません。これらのテキスト表示は、jobリスト内のjobをクリックするたびに更新されます。stdoutまたはstderrのディスプレイをダブルクリックすると、スクロールがしやすくなるようテキスト全体を含むポップアップウィンドウが開きます。

13.1.4 The Display button(ディスプレイボタン)

runボタンとscheduleボタンの下にある【Display:】ボタンでは、各ジョブにおける最重要たる入力ファイルと出力ファイルを視覚的に表示することができます。GUIの下半分にあるjobリストのうちのひとつのjobを選択した後に【Display:】ボタンをクリックすると、選択されたjobの入力と出力がすべて(例えば、パーティクル、顕微鏡写真、座標、 PDFファイルなど)がポップアップメニューとして表示されます。中間体のファイルを表示するといったような一般的な機能はGUI左上の3つのメニューのうちFile>Displayを選択することで使うことができます。

2.1.5 The Job action button(jobアクションボタン)

【Job action】ボタンは、選択した(実行中、完了後、またはスケジュールされた)ジョブの設定を含む小さなメニューを開きます。ここで、note.txtというファイルにアクセスすることができます。このファイルは個々のジョブディレクトリに保存されており、ユーザーのコメントを格納したり、jobのエイリアスを変更したり、jobを終えたものをマークすることができたり、job履歴のフローチャートを作ったり(LATEXとTikZパッケージがシステムにインストールされている場合は、第11章を参照してください)、ディスク容量を節約するためにjobを削除またはクリーンアップ(13.1.6参照)したりすることができます。

2.1.6 Clean-up to save disk space(ディスク領域を節約するためのクリーンアップ)

ジョブを削除すると、ジョブディレクトリ全体がプロジェクトディレクトリからTrash/というディレクトリに移動します。RELIONの左上にある【File】メニューの【Trash】フォルダを空にすると、完全に削除され空き領域が増えます。これを行うまでは、左上の[job]メニューから対応する設定を使用して、ジョブを復元することができます。

ディスク容量を節約するため、ジョブを「clean」にすると、refine jobのすべての中間的な繰り返しに対して書き出されたファイルのような中間ファイルをごみ箱フォルダに移動します。これには2つのクリーニングオプションがあり、片方は【gentle clean】という方法です。この方法では、他のジョブへの入力として使用できるすべてのファイルをそのまま残します。もう一方は【harsh clean】という完全なクリーニング方法です。この方法では全てのファイルを消すことになるため、より多くの容量を確保できます。その中でも【Motion correction】、【Particle extraction】、【Movie refinement】、【Particle polishing】のジョブタイプから選択できる粒子スタック、または顕微鏡写真を含むディレクトリに関しては、とりわけ多くの容量を確保できます。

 また、RELIONの左上にある【job】メニューの対応するオプション(【Gently clean all jobs】あるいは、【Harshly clean all jobs】)を使用して、プロジェクト内のすべてのディレクトリをワンクリックで消去することもできます。 その際、特定のディレクトリを削除対象から除外したい場合には、NO_HARSH_CLEANというファイルをその中に置いてください。

 例えば、保護したい粒子データがjob098にある場合、次のコマンドを実行してください。

 $ touch Polish/job098/NO_HARSH_CLEAN

※ touchコマンドは、指定したファイルが存在しない場合には空のファイルを作成するプログラムです。存在している場合には、ファイルの修正時刻を変更します。


13.2 それぞれの環境におけるけさんを最適化

13.2.1 GPU-acceleration(GPUによる高速化)

Erik Lindahl(ストックホルム)のグループのDari KimaniusとBjoern Forsbergは、RELIONの計算コストの高い部分を移植して、GPUを使うようになりました。彼らはこの作業を行うためにNVIDIAのCUDA-librariesを使用したため、RELION内のGPUアクセレータはNVIDIAカードでしか動作しません。NAVIDIAのGPUはcompute campacity 3.5以上である必要があります。単精度と倍精度の両方のカードが動作するため、高価な倍精度カードに限らず安価なゲーミングカードを使用することもできます。2つの異なるRELIONプログラムがGPUにより高速化されています。それは、relion_autopick(自動ピッキング)とrelion_refine(2次元分類、3次元分類、3次元自動精密化job)です。これらのプログラムのシーケンシャル版とMPI版の両方が加速されました。

13.2.2 Disk access(ディスクへのアクセス)

GPUアクセレータによって画像処理速度があるかに向上したため、ハードディスクへのアクセスがますます妨げられるようになりました。RELION-2.0のGUIには、データセットとコンピュータセットアップのディスクアクセスを最適化するためのいくつかのオプションが用意されています。2次元分類、3次元初期モデル、3次元分類、そして3次元自動精密化のためにパラレルディスクI/Oの使用【use parallel disc I/O】を選択できます。Yesに選択したとき、すべてのMPIプロセスがハードディスクから同時に粒子を読み込みます。そうでないときは、マスターだけがイメージを読み取り、ネットワークを通してこれらをスレーブに送信します。fhgfsのglusterのような並列ファイルシステムは、並列ディスクI/Oで優れています。NFSは並行してたくさんのスレーブを読み込むと中断してしまうでしょう。     メモリ中にプールされた粒子の数【number of pooled particles】を設定することもできます。粒子はMPIスレーブによって個々のバッチに加工されます。各々のバッチの間で、粒子画像のスタックは、一度だけ開閉され、ディスクアクセス時間を改善します。1つのバッチのすべての粒子画像は一斉にメモリに読み込まれます。これらのバッチのサイズは、使われるスレッドあたり、少なくとも1つです。このパーラメータは、スレッドごとにバッチ内で一緒に読み込まれる粒子の数を制御します。もし、3にセットし、1つにつき8つのスレッドが使われるとすると、3x8 = 24個の粒子のバッチが一緒に読み込まれます。これによりディスクアクセス、特にメタデータ処理が問題であるディスクアクセスシステムのパフォーマンスが向上するでしょう。それは、控えめなRAM使用を与えます。もし、比較的小さなデータセット(及び/または大量のRNAを持つコンピュータ)を持つとすると、計算のはじめに全ての粒子をRNAに先行読み取りすることができます【pre-reqd all particles into RAM】。これによりディスクアクセスが比較的遅いシステムの計算を大幅に高速化することが出来るでしょう。しかしながら、使用可能なRNAの量には注意が必要です。粒子はフロート精度で読み込まれるため、N粒子をRNAにお読み込むためには(N×boxsize×boxsize×4)/(1024×1024×1024)GBが必要です。200ピクセルのボックスサイズの100,000個の粒子は15GBとなり、または、400ピクセルのボックスサイズで同じ粒子の数では60GBとなる。データセットが大きすぎてRAMに先読みできないが、各計算ノードに同じ名前でマウントされたローカルの高速ディスク(ソリッドステートドライブなど)(スクラッチディスクと呼ぶ)がある場合、各MPIスレーブはすべての粒子を 計算を開始する前にローカルディスクを使用することができます。これは、粒子をスクラッチディレクトリにコピーすることによって行われます【Copy particles to scratch directory.】。もし、多様なスレーブが同じノード上で実行されると、最初のスレーブだけが粒子にコピーされます。もし、全体のデータセットを保持するのにローカルディスクが小さすぎると、これらの粒子はもはやスクラッチディスクに当てはまらないで、元々の位置から読み取られるでしょう。Relion_volatileと呼ばれるAサブディレクトリは指定されたディレクトリ名の内部につくられます。例えば、/ssdを指定すると、/ssd/relion_volatileという名のディレクトリが作成されます。もし、/ssd/relion_volatileという名のディレクトリが既に存在してたとすると、粒子をコピーする前に消去されます。その後、プログラムはすべての入力粒子を、このディレクトリ内の単一な大きさのスタックにコピーします。もしこのjobが正確に終了すると、/ssd/relion_volatileのディレクトリは再び消去されるでしょう。もし、このjobが終了する前に故障すると、自分で除去したいと思うことでしょう。このプログラムは、すべての人が書き込むことの許可を持って/ssd/relion_volatileを作成します。それにより、スクラッチディレクトリとしてユーザーネーム無しで/ssdを使用することを選択することが出来ます。その方が、各計算ノード上の1人のユーザーが常に1つのjobしか実行しないという条件で、ローカルディスクは、jobがクラッシュし、ユーザーがスクラッチディスクのクリーニングを忘れたときでも、ジャンクでいっぱいになる危険性はありません。最後に、ディスクを通して繰り返しを組み合わせるオプション【combine iterations through disc】があります。Yesに設定すると、すべてのMPIスレーブはすべての繰り返しの最後に、蓄積された結果をひとつの大きなファイルに書き出します。MPIマスターはこれらすべてのファイルを読み込み、それらをすべて結合し、結合された状態で、ひとつの新しいファイルに書き出します。すべてのMPIスレーブは、その後、その結合された結果を読み込みます。これにより、ネットワークの負荷が軽減されますが、ディスクI / Oの負荷が増加します。これは、期待値ステップ(E-step)の終わりに達する進捗バー(マウスがチーズに達する)とその後の最大化ステップ(M-Step)の開始にかかる時間に影響します。この時間は、非常に効率的なシステム設定に依存します。このオプションは当初、LMBでの古いクラスタ上でのネットワークカードのバグを回避するために実装されました。 今日では、精密化が高解像度に達すると非常に遅くなる傾向があるので、このオプションは使用しない方が望ましいです。


13.3 Scheduling jobs(jobのスケジュール化)

【Schedule】ボタンから、今後のjobの準備をできます。スケジュール化された動作により期待される出力ファイルは、既に【Browse】ボタンから他のjobへの入力として利用できるため、Scheduled jobsに複数の連続した動作の「パイプライン」を構築することが可能です。GUIの左上にある【Autorun】メニューからスケジュール化されたjobを実行することができます。そのプログラムは何回スケジュール化された動作を実行したいかを訊ねてくるでしょう。1より大きい値を指定すると、すべてのスケジュールしたjobが繰り返し実行されます。この設定にした場合、スケジュール化されたjobが一周毎に、何分経過させるべきかという質問が出てきます。(もしスケジュールしたjobの実行時間がこの時間よりも短い場合、このプログラムは、次にスケジュールしてあるjobの実行を待ちます。)jobタイプが、【Import】、【Motion correction】、【CTF estimation】、【Auto-picking】、【Particle extraction】、【Movie refinment】である場合には、以前の反復で完了していない顕微鏡写真のみを処理します。他のすべてのjobタイプでは、スケジュール化されたjobの各反復後、最初から計算を再スタートするでしょう。

13.4 On-the-fly processing(撮影しながら処理を実行)

これから紹介するRELION-3.0はスケジュール化されたjobをパイプラインに追加することと、これらをコマンドラインから実行うることを可能にします。これによって進化を可能にされたRelion_it.pyと呼ばれるスクリプトは、RELIONに組み込まれた/scriptsのディレクトリの中で見つけることが出来ます。このスクリプトは次の手順で自動的に実行される。   ●*.mrcsか*.mrcか*/tiffのすべてのファイルの入力は、動画と平均化された顕微鏡写真(*.mrcのみ)のどちらにもなり得ることができます。大きなファイルをコピーするのにはたくさんの時間がかかり、まだコピーが終わっていない動画や顕微鏡写真を入力するのは好ましくないため、このスクリプトはファイルを一時的なファイル名でコピーし、コピーが完全に終了したときにのみ、ファイルを適切なファイル名に変更するべきでしょう。   ●入力されたファイルが動画の場合:Motion correction

  ● CTF estimation
  ● Auto-pickingはRELION-3.0の新しい特徴である:reference-freeからピックされLaplacian-of Gaussian(LoG)フィルターに基づいています

 ●Particle extraction

 この手順はたくさん繰り返され、そして、各反復の間に追加された新しい画像だけが加工されるでしょう。一方では、スクリプトがいくつの分子が引き抜かれたかを同様に観察しています。ユーザー指定の数の粒子が抽出されると、スクリプトは、ユーザー指定の数の粒子のバッチで、2D分類/3D分類jobも起動します。3D参照がまだ利用できない場合、スクリプトは3D初期モデルjobを使用して自動的に生成することもできます。必要に応じて、スクリプトはすべての顕微鏡写真をもう一度通過します。そのため、最初の分子から適切な3D参照が取得されているので、これを使用してテンプレートベースの自動ピッキングjobの投影を生成します。おそらく最初のパスで使用されたものより大きいBoxのサイズを有するピッキングランから抽出された粒子は、その後2D分類/3D分類jobに再度使用することができます。セットアップ全体では、収集中の画像をその場で簡単に処理することができます。また、プロジェクトディレクトリでGUIを開くことにより、プロセスを追跡することができます。顕微鏡でのデータ収集中にすでに2次元クラスの平均を調べることは、データを最初に取得する価値があるかどうか、または、方位分布が可能か、または、氷の厚い領域で収集する必要があるかどうかを判断するために役立ちます。オンザフライ3D再構成は、例えば、ある補因子が関心のある複合体に結合しているかどうかを確認するために有用であり得る。パイソンで書かれているので、スクリプトを修正するのは比較的簡単です。そして、私たちはこのスクリプトが多くのユーザーのニーズと設定を反映するために多くの異なる方法で修正されるかもしれないと思われます。例えば、取得されているデータセットの達成可能な解像度を評価するために3D自動更新jobを非常に簡単に追加できます。


13.5 Heerical reconstruction(ヘリカル再建)

Scheresグループの博士課程学生であるShaoda Heは、ヘリカルアセンブリ処理のためのワークフローを実装しています。これには、【Auto-picking】、【Partcle extraction】、【2Dclassification】、【3Dclassification】、【3D auto-refine】、【Particle polishing】、および【Mask create】といったjobタイプのパラメータパネルへの追加タブが含まれます。ヘリ    カル再建を処理するための個別のチュートリアルはありません。一般的な原則は、このチュートリアルで取り上げている単粒子分析と同じです。したがって、ヘリカル処理にRELIONを使用するつもりのユーザーは、このチュートリアルを最初に実行することをお勧めします。これらオプションのより詳細な説明は、ユーザーがRELION Wikiの関連するページ、またはShaodaの論文から参照できます。

13.6 Sub-tomogram averaging(サブ断層写真の平均化)

MRC-LMBのLoweグループの元博士号であるtanmay Bharatの助けを借りて実装されたサブ断層撮影の平均化。一般的な概念の多くは単一粒子解析と同じであり、RELIONで顕微鏡画像の平均化を実行するユーザーは、このチュートリアルを先に進めることをお勧めします。サブ断層画像平均化の手順の詳細な説明については、ユーザーはRELION wikiの対応するページ、またはTanmayの論文を参照してください。なお、サブ断層像平均化パイプラインをより使いやすく、よりよくするための努力を続けています。この作業は、MRC-LMBでも、John Briggsのグループと密接に協力して行われています。一方、サブ断層像平均化パイプラインは単一粒子のものよりもかなり合理化されておらず、多くの場合、RELIONパイプラインの外側でスクリプトを作成する準備をする必要があります。

13.7 Interaction with other programs(他のプログラムとの相互)

原則として、RELIONは異なるプログラムによって抽出されたパーティクルを使用できますが、これは推奨される手順ではありません。多くのプログラムは粒子自身を変更します。位相反転、バンドパスまたはウィーナーフィルタリング、マスキングなどを介して行われます。これらはすべて、その後のRELIONでの使用には最適とは言えません。さらに、必要なメタデータすべてを正しくフォーマットされたRELIONタイプのSTARファイルにまとめると、エラーが発生する可能性があります。 RELIONでパーティクルを再抽出するのは簡単で非常に高速であるため、以下に概説する手順は多くの場合RELIONへのはるかに簡単な(そしてより良い)経路です。また、RELIONを囲むラッパーの実装もいくつか報告されました(例えばEMAN2、SCIPION、APPIONなど)。他の人がこれらのラッパーを書くとき私たちは役に立つように努めますが、私たちはそれらを完全に制御することができず、彼らの最終製品が最良の方法でRELIONを使うかどうか知りません。したがって、これらのラッパーで得られた結果に疑問がある場合は、このチュートリアルで概説している手順に従うことをお勧めします。