RELION3チュートリアル

提供: Eospedia
移動: 案内検索

前書き

  • 本ドキュメントは、RELION3のチュートリアルを日本語訳し、図や補足説明などを追加して充実化したものです(Sjors Scheres博士の許可を得ています)。
  • 誤訳等がある場合には、ご連絡いただけますと幸いです。少しずつ、良いものにしたいと思います。
  • 書き進めてみて気づきましたが、どこまでが本家チュートリアルに書いてあることでどこからが注釈、補足なのか、分かりにくいと思います。あまり良くない気がしているので、随時切り分けを進めます。
    • 少なくとも ※ 記号で始まるところは、日本語訳に伴って追加した注釈、補足です。


実行環境

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のプロジェクトディレクトリの中に動画を置きたくない場合は、動画を格納したディレクトリのシンボリックリンクをプロジェクトディレクトリ内に作成すれば良いです。

$ ln -s xxxx Movies/yyyy

画像ファイル名の拡張子は .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


※de novo 3D構造決定は50%の確率で真の3D構造の鏡像を出してきます。今回はたまたま正しい鏡像で出たためPDBがうまくハマっていますが、そうでない場合はうまくはまらない部分が出てきます。PDBフィッティングはこのチュートリアルの目的ではないので(なんとなくきれいな図を載せたくてやってみただけです)、鏡像によりPDBがハマらなくても気にせず、先へ進んで問題ありません。鏡像の問題(handedness)に関しては #10 Checking the handedness を参照ください。


可視化により対称性があることに気づいたら、対称性の軸が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


※de novo 3D構造決定で真の構造の鏡像が出てしまった場合、それを初期構造として使ったリファインメント(3D classification, 3D auto-refine)でも同じく鏡像が出てしまいます。今回はたまたま正しい鏡像で出たためPDBがうまくハマっていますが、そうでない場合はうまくはまらない部分が出てきます。PDBフィッティングはこのチュートリアルの目的ではないので(なんとなくきれいな図を載せたくてやってみただけです)、鏡像によりPDBがハマらなくても気にせず、先へ進んで問題ありません。鏡像の問題(handedness)に関しては #10 Checking the handedness を参照ください。


クラス選択

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粒子を捨てたことになります。ヘテロな構造を見つけ出すというより、ゴミ落とし的に働いた感じです。

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


※de novo 3D構造決定で真の構造の鏡像が出てしまった場合、それを初期構造として使ったリファインメント(3D classification, 3D auto-refine)でも同じく鏡像が出てしまいます。今回はたまたま正しい鏡像で出たためPDBがうまくハマっていますが、そうでない場合はうまくはまらない部分が出てきます。PDBフィッティングはこのチュートリアルの目的ではないので(なんとなくきれいな図を載せたくてやってみただけです)、鏡像によりPDBがハマらなくても気にせず、先へ進んで問題ありません。鏡像の問題(handedness)に関しては #10 Checking the handedness を参照ください。


※ _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

※de novo 3D構造決定で真の構造の鏡像が出てしまった場合、それを初期構造として使ったリファインメント(3D classification, 3D auto-refine)でも同じく鏡像が出てしまいます。今回はたまたま正しい鏡像で出たためPDBがうまくハマっていますが、そうでない場合はうまくはまらない部分が出てきます。PDBフィッティングはこのチュートリアルの目的ではないので(なんとなくきれいな図を載せたくてやってみただけです)、鏡像によりPDBがハマらなくても気にせず、先へ進んで問題ありません。鏡像の問題(handedness)に関しては #10 Checking the handedness を参照ください。


分解能の計測は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が大きくなるに連れてデフォーカスが大きくなる傾向がありますので、わずかに氷の層が傾いているかもしれません。

※ごくまれに外れ値的なデフォーカス量を持つ単粒子がいたりします。そのような単粒子はSubset selectionで除外するなどした方がいいと思います。

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をやり直してもいいでしょう。ここでは特にそうしたことはせず、このまま先へ進むことにします。


STARファイル

Ctf refinementでは、particles_ctf_refine.star というparticle STARファイルが生成されます。CTF refiementにより単粒子毎にデフォーカス量(rlnDefocusU, rlnDefocusV)が更新されており、またビーム傾斜量( rlnBeamTiltX, rlnBeamTiltY)のカラムが追加されています。

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

RELION3.0より、ビーム誘起性の試料ドリフトをベイズ推定により補正するアルゴリズムがBayesian polishingジョブタイプとして実装されています。これは正則化尤度を最大化(※事後確率最大化)するような単粒子の軌跡の集合を推定することにより行います。それにより、個々の単粒子の軌跡に対し、空間的にcoherent(※画像中の近いところで全然違うドリフトになったりしない、ということ?)かつ時間軸方向に滑らかであるという事前情報を、ガチガチな拘束条件を課したりすることなく持たせることができます。滑らかさを表現する事前情報は3つあります(ドリフトに関する統計量)。今回のデータセットで起きているドリフトの性質を最もよく表現する事前情報を推定するため、まずはBayesian polishingジョブを訓練モードで実行します。訓練モードにより事前情報を推定できたら、再度プログラムを走らせて(ポリッシングモード)データセット内のすべての単粒子の軌跡を推定し、適切に重み付けされたドリフト補正動画フレームを生成することができます。

※こちらがBayesian polishingの論文です → http://journals.iucr.org/m/issues/2019/01/00/fq5003/index.html

8.1 Running in training mode(訓練モードでの実行)

16スレッドで並列実行すると、1時間15分かかります(※本家の環境で)。これを待つのが嫌であれば、訓練モードは飛ばして次節のポリッシングモードに進み、すでに計算済みのσ値(※訓練モードで推定したいのはこれ)で実行すれば良いです。計算済みσ値は次節に書いてあります。

このジョブを自分で実行する場合は、Bayesian polishingジョブタイプを選択し、以下のようにパラメータを入力します。

rln-polish-train-01.png

rln-polish-train-02.png

rln-polish-train-03.png


Bayesian polishingジョブ訓練モードのパラメータ
タブ パラメータ名 説明
I/O Micrographs (from MotionCorr) MotionCorr/job002/corrected_micrographs.star Motion correctionジョブで生成されたSTARファイル。RELION3以上で実行したMotion correctionジョブでなければならない。それ以前のバージョンでは、Bayesian polishingに必要なメタデータが書き出されていない。(※ RELION3でUCSF版のMotionCorr2使った場合はどうなるのだ?)
I/O Particles (from Refine3D or CtfRefine) CtfRefine/job024/particles_ctf_refine.star ここで指定した単粒子データセットがポリッシングにかけられる。
I/O Postprocess STAR file PostProcess/first3dref/postprocess.star ポストプロセッシングにより生成されたSTARファイル。そのポストプロセッシングで使われたマスクと、計測されたFSC曲線が、ポリッシング処理で使われる。
I/O First movie frame 1 ドリフトのフィッティングに使用する最初のフレーム番号。動画の最初の方のフレームや最後の方のフレームを捨てる人もいるが、Bayesian polishingを行う場合はそれは推奨されない。フレームごとのB-factorによる重み付けによりポリッシュ後の単粒子のS/Nが自動的に最適化されるので、すべての動画フレームを含めるのが良い。(※ ....?)
I/O Last movie frame -1 -1は最終フレーム。説明は上記参照。
Train Train optimal parameters? Yes Yesとすると、relion_motion_refineが入力データの一部(Use this many particlesで指定)を使って、事前情報(3種類のσ値)の推定を実行する。
Train Fraction of Fourier pixels for testing 0.5 ここではとりあえずデフォルト値を使用。This fraction of Fourier pixels (at higher resolution) will be used for evaluation of the parameters (test set), whereas the rest (at lower resolution) will be used for parameter estimation itself (work set). (※ ....? フーリエ空間のピクセル全体のうち、高分解能側のこの割合分だけのピクセルを、推定したパラメータ(3つのσ値?)の性能評価のために使用。パラメータの推定自体は残りの低分解能側の情報を使って実行する?)
Train Use this many particles 5000 パラメータ推定に使用する単粒子数。多ければ多いほど正確になるが、メモリ使用量と計算時間が増える。メモリが足りなくなった場合は単粒子数を減らして学習する。ただし一般に 5,000 粒子をかなり下回るような数での実行は推奨されない。いずれにせよ今回は5,000粒子ほどしか無いので、とりあえず5,000粒子で実行する。
Polish Perform particle polishing? No まずは学習モードでの実行だけ。


このプログラムの学習モードはMPI並列化はされていません。したがって単一MPIプロセスで実行するようにしてください(※スレッド並列化は可能です)。 エイリアス名を train とし、【Run!】を押して実行してください。

(※本家の環境で)高速化のために16スレッド並列で実行しましたが、それでも1時間以上かかりました。

(※すごく時間かかります。メモリ消費ですが、少なくとも 22 GBは行きます。)


※計算が終了すると、標準出力に以下のような表示がされました。本家チュートリアルでは s_vel = 0.462, s_div = 1935, s_acc = 2.46 なので、結構違う推定結果になっている感がありますが、そういうものでしょうか?

good parameters: --s_vel 0.3795 --s_div 1230 --s_acc 3.42

written to Polish/job050/opt_params.txt

※以下は relion_motion_refine --help して表示されるヘルプおよびBayesian polishingの文献から推測(憶測?)した説明です

  • s_vel( σV )はドリフトの速さの標準偏差 [Å/dose]
    • 単位はこういうことですかね? → [Å/(e-2)]
    • なんだかわかりにくいですが、フレームあたりで発生しているドリフト量(フレーム内での平均ドリフト速さ)の標準偏差[Å/frame]がまずあって、それを単位照射量あたりで正規化したものと考えればいいですかね?
  • s_div( σD )は、ドリフト速度の相関長(空間方向)[Å]を表す量
    • 文献の式(2)を見るに σD が大きいほど、距離の近い単粒子同士のドリフト速度の相関が高くなる、すなわちおんなじような動き方をするみたいな感じでしょうか。
  • s_acc ( σA )はドリフトの加速度の標準偏差 [Å/dose]
    • 加速度なのに単位が速さと一緒で..? ですが、隣接フレーム間のドリフト速度の変化の大きさの標準偏差[Å/frame]があって、それを単位照射量あたりで正規化したと考えばいいですか?

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

学習モードが終了すると、Polish/train/opt_params.txt というファイルが書き出されます。そこに記録されたパラメータを用いてポリッシング処理を行うためには、ジョブタイプブラウザからBayesian polishingジョブタイプを新規選択します。(※学習モードで使用したパラメータがそのまま残っているはずです。)

I/Oタブは学習モードと同じままにし、TrainタブおよびPolishタブに以下のようにパラメータを入力します。

rln-polish-01.png

rln-polish-02.png


Bayesian polishingジョブ ポリッシングモードのパラメータ
タブ パラメータ名 説明
Train Train optimal parameters? No すでに学習済みなのでNo
Polish Perform particle polishing? Yes relion_motion_refineが以下のパラメータのもとで単粒子毎のドリフト軌跡推定処理を実行し(ポリッシング)、ポリッシュされた単粒子画像を出力する。
Polish Optimised paramter file Polish/train/opt_params.txt 学習モードで生成された、事前情報パラメータを記録したテキストファイル。
Polish OR use your own parameters? No 自分で事前情報パラメータを入力する場合はYesとする。その場合上記のOptimised parameter fileは無視される。
Polish Minimum resolution for B-factor fit (A) 20 B-factorのフィッティングで使用される最小の分解能 [Å]。今回はデフォルト値を使用。
Polish Maximum resolution for B-factor fit (A) -1 B-factorのフィッティングで使用される最大の分解能 [Å]。負の値を指定すると、FSC曲線(Postprocessingの)から自動で決定される。(※具体的には?)。今回はデフォルト値を使用。


もし学習モードをスキップした場合は、PolishタブでOR use your own parameters? を Yes にし、予め計算してある(※本家の方で)σパラメータを入力してください。

  • Sigma for velocity (A/dose) : 0.462
  • Sigma for divergence (A) : 1935
  • Sigma for acceleration (A/dose) : 2.46


このプログラムはMPI並列化されています。エイリアス名を polish とし、【Run!】を押して実行してください。

3 MPIプロセス、16 スレッド/プロセスで、2分で終わりました(※本家では)。

8.3 Analysing the results

Bayesian polishingジョブタイプは、ポリッシングされた単粒子の情報を記録したSTARファイル shiny.star と、PDFログファイル (logfile.pdf)を出力します。後者は電子線照射損傷に基づく重み付けのためのスケール(※ ??)およびB-factorのプロット、全画像の全単粒子のリファインメント後のドリフト軌跡プロットを含んでいます。

【Display:】でout:logfile.pdfを選択し、ログファイルを見てみます。

ドリフト軌跡

※例として1枚だけ示しています。青丸が各単粒子の最初のフレームにおける位置で、黒いにょろにょろがドリフトの軌跡です(30倍に拡大されています)。 rln-polish-03.png


単粒子のドリフト軌跡のプロットを見ると、今回のデータセットでは、試料ステージが少しドリフトしていたようです。ほとんどすべての単粒子が動画の中で右上から左下へ向けて動いているからです。


フレームごとのB-factor推定

※ドリフト推定後にフレームのアラインメントを行った後、パーティクルポリッシングの文献に記載の方法で、参照像との比較に基づきフレームごとのB-factorが推定されます。B-factorは符号で混乱しそうになりますが、多分下のプロットのB-factorが exp(B*(k**2)) 的な感じで各フレームのフーリエ変換に掛け算され、ローパスフィルタ的に働くのではないかと思います。電子線損傷が蓄積しているあとの方のフレームほどB-factorが小さくなっているためより強いローパスフィルタとなり、つまり後の方のフレームほど高周波数成分が弱く抑えられる形になるかと思います。違っていたらご指摘ください。

rln-polish-04.png


※もう一つ、"Polishing scale-factors"なるプロットも見れます。何を読み取れば良いのかよくわかりません。。

8.4 3D refinement after polishing

※ 本家チュートリアルでは8.3節に軽めに描かれていたところですが、ここでは勝手に別の章として書いています。

ポリッシングにより単粒子画像のS/Nが向上したので、再度3D auto-refineと Post-processingを実行しましょう。

3D auto-refineジョブタイプを選択し、以下の通りパラメータを入力します(載せていない分については前回の 3D auto-refineと同じです。)

rln-refine-polish-01.png

rln-refine-polish-02.png

rln-refine-polish-03.png


一部説明が必要なパラメータだけ以下に示します。

ポリッシング後の3D auto-refineのパラメータ(一部)
タブ パラメータ名 説明
I/O Input images STAR file Polish/polish/shiny.star ポリッシングされた単粒子を指定。
I/O Reference map Refine3D/first3dref/run_class001.mrc 前回の3D auto-refineの結果マップを指定。
I/O Reference mask (optional) MaskCreate/job037/mask.mrc ポストプロセッシングで使用したマスクと同じ。このオプションを使うことにより参照像にマスクが施され、マスクの外側がゼロ値で埋められる。参照像のノイズを減らせるため投影パラメータ推定の精度が向上し、結果として再構成の分解能も向上する。
Reference Ref. map is on absolute greyscale? Yes 今回はYesとして良いはず。本家チュートリアルでは特に言及されていないのでNoでもいいかもしれない。RELION GUIのはてなマークで表示されるヘルプメッセージでも、"Noにしても結果が悪くなることはないです"的なことも書いてある。
Optimisation Use solvent-flattened FSCs? Yes Yesとすることで、各イテレーションでgold-standard FSCを計算する際、マスク(Reference maskで指定)をかけた状態でFSCが計算される。Post-processingとよく似た方法が使われる(※よく似た、ということは、全く同じではない?何がどう違う?)。このオプションは、目的蛋白質が小型で単粒子ボックスの中の比較的小さい領域しか占めていない場合に特に有用である(例えば、細長い分子や、マスクを使ってごく小さい領域のみリファインメントするときなど)。3D auto-refineでFSCを計算する際のデフォルトの挙動は、マスク無しでハーフマップ同士のFSCを計算することである。このオプションにより溶媒ノイズの影響をなくしたFSCを計算できるため、分解能が向上する。


エイリアス名は polish とし、【Run!】で実行してください。

※私の環境での標準出力を載せておきます → Polish後の3D auto-refineの標準出力


ポリッシングされた単粒子を使ったこと、参照像に溶媒ノイズをカットするマスクをかけたことで、分解能が 3.1 Å (Postprocess) → 2.9 Åに向上しました。


B-factorでシャープニングしたマップを得るべく、ポストプロセッシングに進みます。I/Oタブの One of the 2 unfiltered half-maps だけ、先程の Refine3D/polish/run_half1_class001_unfil.mrc を指定します。他は前回のポストプロセッシングと同じで良いでしょう。

出力は以下のようになりました。

== Reading input half-reconstructions: 
  + half1-map:                     Refine3D/job052/run_half1_class001_unfil.mrc
  + half2-map:                     Refine3D/job052/run_half2_class001_unfil.mrc
== Using a user-provided mask ... 
  + input mask:                    MaskCreate/job037/mask.mrc
== Masking input maps ...
  + randomize phases beyond:       9.65042 Angstroms
== Dividing map by the MTF of the detector ...
  + mtf STAR-file:                 mtf_k2_300kV.star
== Applying sqrt(2*FSC/(FSC+1)) weighting (as in Rosenthal & Henderson, 2003) ...
== Fitting straight line through Guinier plot to find B-factor ...
  + fit from resolution:           10
  + fit until resolution:          0
  + slope of fit:                  -3.48274
  + intercept of fit:              -12.6341
  + correlation of fit:            0.302465
  + apply b-factor of:             -13.931
== Low-pass filtering final map ... 
  + filter frequency:              2.86904
== Writing output files ...
  + Processed map:                 PostProcess/job055/postprocess.mrc
  + Processed masked map:          PostProcess/job055/postprocess_masked.mrc
  + Metadata file:                 PostProcess/job055/postprocess.star
  + FINAL RESOLUTION:              2.86904


※同じマップ、マスクを使っているので分解能は特に上がったりはしません。(なぜか微妙に違いますが)

※ポストプロセッシング前にマスクを再設計し、分解能のちょっとした向上を狙いに行くのもアリです。


マップは以下のようになりました。

rln-refine-polish-04.png


3 GB のデータから再構成したにしては、悪くない結果ですね!

8.4 When and how to run CTF refinement and Bayesian polishing

ベイジアンポリッシングもCTFリファインメントも、3D再構成の分解能を向上させるでしょう。ここで一つ疑問が出てきます。果たしてどちらを先にやるべきかと。今回の例では我々はまずCTFリファインメントにより単粒子毎のデフォーカス値およびデータセット全体のビーム傾斜を最適化しましたが、しかし先にベイジアンポリッシングをやることもできたでしょう。どちらの処理もより高分解能のモデルがあればよりよい結果を返すので、入れ代わり立ち代わりこれらのジョブを繰り返すと良いかもしれません。今回であれば、ベイジアンポリッシングをした跡に再度CTFリファインメントを行うなど。一般的には、より大きな問題がある方から先に実行するのが良いですが、いずれにせよtrial & errorでいろいろやってみるしか無いでしょう。(※例えばドリフトがとても激しいデータならベイジアンポリッシングを先にやってドリフト補正したほうがいいかもしれませんし、試料の屈曲や氷の厚みによりデフォーカスばらつきが大きい、ビームシフトでデータを取ったためにビーム傾斜が大きいデータなどならCTFリファインメントを先にやったほうがいいかもしれません。)


もう一つ言っておくと、ベイジアンポリッシングの訓練モードの推定値は再現性が取れない場合があります。つまり、複数回訓練モードを走らせたとき、推定されるσ値が大きく異なるなど。しかし、ポリッシングで使用されるσ値の値は最終的な結果の分解能にあまり影響を与えない場合が多いということもわかっています。それに加えて訓練モードは計算に時間がかかるという事情もあるので、訓練モードは実行せずに、デフォルトのσ値(σV = 0.2; σD = 5000; σA = 2) を使って直接ポリッシングモードで計算してしまうというのも有りでしょう。

9 Local-resolution estimation (局所分解能の推定)

ポストプロセッシングで推定された分解能は、マスク内の3D構造全体のグローバルな推定結果です。しかし巨大分子複合体の3D構造の中では、場所によって分解能が異なることがしばしば有り、その場合はただ1つの数字だけでは分解能を表現できません。Alp KucukelbirとHemant Tagereは、3D構造内における分解能のバリエーションを推定するプログラムを開発しました(Alp Kucukelbir et al., Nature Methods (2014))。RELIONでは、Local resolutionジョブタイプからそのプログラムを使用することができます。あるいは、RELION独自のアルゴリズムとして、3Dマップ内を移動しながら各所でソフトな球形マスクを用いてポストプロセッシング的な処理により局所的な分解能を測っていくこともできます。今回は後者の方法を使ってみます。

9.1 Running the job

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

rln-locres-01.png

rln-locres-02.png

rln-locres-03.png

Local resolutionジョブタイプのパラメータ
タブ パラメータ名 説明
I/O One of the 2 unfiltered half-maps Refine3D/polished/run_half1_class001_unfil.mrc 局所分解能を計測したい3Dマップのハーフセットの片方を指定。
I/O Calibrated pixel size (A) 1.244 説明はpostprocessingの節を参照。
ResMap Use ResMap? No RELION実装版の局所分解能推定を使うのでここはNo
Relion Use Relion? Yes RELION実装版を使用
Relion User-provided B-factor -30 局所的分解能に基づいて局所的に周波数フィルタリングおよびシャープニングを施したマップを作成する際に使用されるB-factorの値。ポストプロセッシングジョブでの推定値(postprocess.starの data_generalブロック _rlnBfactorUsedForSharpening)に近い値を入れるのがおそらくよい。マップをシャープニングしすぎるとノイズを観察、解釈してしまうことになるので重々注意。
Relion MTF of the detector (STAR file) mtf_k2_300kV.star ポストプロセッシングのときに使用したのと同じMTFファイル。

(※本家チュートリアル(ver 3.0)ではI/Oに"User-provided solvent mask"を指定していますが、RELION 3.0.5 ではそのようなオプションはありません(ResMapの方にはありますが。))


エイリアス名として polished を指定しておきましょう。

RELION版の局所分解能推定はMPI並列実行可能です(ResMapは無理)。8 MPIプロセス並列で 7 分程度かかります。結構時間がかかるのでできたら並列実行したほうが良いです。【Run!】を押し、実行してください。

9.2 Analysing the results

RELION版は2つのmrcファイルを出力します。

  • relion_locres.mrc
    • 入力した3Dマップと同じサイズで、各ボクセルには局所分解能の推定値がÅ単位で入っています。
  • relion_locres_filtered.mrc
    • (RELION版のみ)局所分解能に応じてフィルタリング・シャープニングされた密度マップです。

UCSF Chimeraを使って、密度マップを局所分解能に応じて色付けすることができます。ポストプロセッシングの結果マップ(Postprocess/polished/postprocess.mrc)を色付けしてもいいですし、relion_locres_filtered.mrcを色付けしてみても良いです。

まず、UCSF Chimeraで色付けしたいマップを開きます。ここではrelion_locres_filtered.mrcを開きます。次に、メニューバーからTools > Volume data > Surface colorを選択します。

rln-locres-04.png


すると密度マップの等高面の色付けに関するコントロールウィンドウが開きます。以下に貼り付けたキャプチャ画面の通り設定すれば概ね同じような見た目になると思います。細かいオプションはウィンドウ下部の Options を押すと表示されます。Colorボタンを押せば色づけされます。(ボクセルに局所分解能の推定値を格納した relion_locres.mrcを使って、relion_locres_filtered.mrcの等高面を色付けしています。)

rln-locres-05.png


上記では 4 Å分解能を赤、2.5 Å分解能を青、その中間を白で表示してみました。溶媒、溶媒-空気界面、他の分子などと接する蛋白質外側表面は局所分解能が低めで、蛋白質内側へ行くほど局所分解能が高いことがわかります。よくあるパターンです。

10 Checking the handedness

マップを良く見てみましょう。handednessが間違っていたりしませんか(※別の言い方をすれば、本当の構造の鏡像が出たりしてないですか)。例えば、アルファヘリックスをよく見ると、右巻きではなく左巻きになってしまっているなど。単粒子解析では、データセットだけからhandednessを決定することはできません(人間がヘリックスの巻方向を見たり、既知の原子モデルを当てはめたり、あるいは試料ステージを傾けた写真を撮ったりしない限り)。それ故、3D initial modelジョブタイプのSGDアルゴリズムで再構成した像は、50%の確率でhandednessが間違っています(※今回のチュートリアルではたまたま正しいhandednessが出ました。)ポストプロセッシング後のマップのhandednessを反転するには、以下のコマンドを実行すればよいです。

$ relion_image_handler --i PostProcess/polished/postprocess.mrc  --o PostProcess/polished/postprocess_invert.mrc --invert_hand


このコマンドは任意のマップに対して使えます。画像処理の早い段階でhandednessが間違っていることに気づいたら、その段階で上記コマンドにより修正した方が良いでしょう。RELIONで再構成をしている分にはどちらのhandでも問題ありませんが、気づいた段階で早めに直した方がいろいろ楽でしょう(※)。

正しいhandであることが確認できたら、βガラクトシダーゼの原子モデルをUCSF Chimeraでマップに重ねて表示してみましょう。PDB-ID 5a1a をfetchするのが手っ取り早いです。(※この日本語チュートリアルではかなり初期の段階から先走ってしまいましたが。)


※参照像のhandednessを反転してリファインメントすることについて

※本家チュートリアルには無い項目ですが気になったので書きのこしておきます。

例えば 3D classificationの後にマップのhandednessが間違っていると気づいたとして、その後の3D auto-refineの前にマップのhandednessを反転しそれを参照像として使ったとしても、単粒子の投影パラメータの方は3D classificationにおいてhandedness反転前のマップに対して最適化されているので、そこをスタート地点にしてリファインメントしてもうまくいかないんじゃない?と思うかもしれません。

そう思って実際に実験してみました。

  • Class3D/first_exhaustive/run_it025_class002_box256.mrc のhandednessを反転したマップを作成
  • それを参照像として 3D auto-refineを実行 (参照像以外のパラメータはすべてRefine3D/first3drefと同一)
  • 3D auto-refineの経過は以下のようになった。
 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= 48.3435 degrees; and in offsets= 2.6562 pixels (no gain for 0 iter) 
 Auto-refine: Estimated accuracy angles= 13.23 degrees; offsets= 3.865 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.4101 degrees; and in offsets= 1.94711 pixels (no gain for 0 iter) 
 Auto-refine: Estimated accuracy angles= 7.045 degrees; offsets= 2.467 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.4818 degrees; and in offsets= 2.01552 pixels (no gain for 1 iter) 
 Auto-refine: Estimated accuracy angles= 2.57 degrees; offsets= 0.993 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.65454 (no gain for 0 iter) 
 Auto-refine: Changes in angles= 5.62055 degrees; and in offsets= 1.00326 pixels (no gain for 0 iter) 
 Auto-refine: Estimated accuracy angles= 1.958 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 1 iter) 
 Auto-refine: Changes in angles= 5.10292 degrees; and in offsets= 0.76261 pixels (no gain for 0 iter) 
 Auto-refine: Estimated accuracy angles= 1.961 degrees; offsets= 0.815 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.30461 degrees; and in offsets= 0.758337 pixels (no gain for 1 iter) 
 Auto-refine: Estimated accuracy angles= 1.97 degrees; offsets= 0.819 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.90904 degrees; and in offsets= 0.706569 pixels (no gain for 0 iter) 
 Auto-refine: Estimated accuracy angles= 2.044 degrees; offsets= 0.818 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.23688 degrees; and in offsets= 0.732301 pixels (no gain for 1 iter) 
 Auto-refine: Estimated accuracy angles= 1.95 degrees; offsets= 0.821 pixels
 Auto-refine: Angular step= 3.75 degrees; local searches= false
 Auto-refine: Offset search range= 3.66151 pixels; offset step= 0.61575 pixels
 Auto-refine: Iteration= 10
 Auto-refine: Resolution= 8.61081 (no gain for 0 iter) 
 Auto-refine: Changes in angles= 3.90627 degrees; and in offsets= 0.580177 pixels (no gain for 0 iter) 
 Auto-refine: Estimated accuracy angles= 1.299 degrees; offsets= 0.587 pixels
 Auto-refine: Angular step= 3.75 degrees; local searches= false
 Auto-refine: Offset search range= 3.66151 pixels; offset step= 0.61575 pixels
 Auto-refine: Iteration= 11
 Auto-refine: Resolution= 8.38421 (no gain for 0 iter) 
 Auto-refine: Changes in angles= 2.25653 degrees; and in offsets= 0.377236 pixels (no gain for 1 iter) 
 Auto-refine: Estimated accuracy angles= 1.219 degrees; offsets= 0.569 pixels
 Auto-refine: Angular step= 3.75 degrees; local searches= false
 Auto-refine: Offset search range= 3.66151 pixels; offset step= 0.61575 pixels
 Auto-refine: Iteration= 12
 Auto-refine: Resolution= 8.38421 (no gain for 1 iter) 
 Auto-refine: Changes in angles= 2.30402 degrees; and in offsets= 0.354229 pixels (no gain for 2 iter) 
 Auto-refine: Estimated accuracy angles= 1.178 degrees; offsets= 0.558 pixels
 Auto-refine: Angular step= 1.875 degrees; local searches= true
 Auto-refine: Offset search range= 1.77115 pixels; offset step= 0.4185 pixels
 Auto-refine: Iteration= 13
 Auto-refine: Resolution= 7.58571 (no gain for 0 iter) 
 Auto-refine: Changes in angles= 1.66164 degrees; and in offsets= 0.288858 pixels (no gain for 0 iter) 
 Auto-refine: Estimated accuracy angles= 0.935 degrees; offsets= 0.473 pixels
 Auto-refine: Angular step= 1.875 degrees; local searches= true
 Auto-refine: Offset search range= 1.77115 pixels; offset step= 0.4185 pixels
 Auto-refine: Iteration= 14
 Auto-refine: Resolution= 7.77073 (no gain for 1 iter) 
 Auto-refine: Changes in angles= 1.17514 degrees; and in offsets= 0.22471 pixels (no gain for 1 iter) 
 Auto-refine: Estimated accuracy angles= 0.851 degrees; offsets= 0.432 pixels
 Auto-refine: Angular step= 0.9375 degrees; local searches= true
 Auto-refine: Offset search range= 1.12355 pixels; offset step= 0.324 pixels
 Auto-refine: Iteration= 15
 Auto-refine: Resolution= 7.58571 (no gain for 0 iter) 
 Auto-refine: Changes in angles= 0.763975 degrees; and in offsets= 0.163416 pixels (no gain for 0 iter) 
 Auto-refine: Estimated accuracy angles= 0.754 degrees; offsets= 0.4 pixels
 Auto-refine: Angular step= 0.9375 degrees; local searches= true
 Auto-refine: Offset search range= 1.12355 pixels; offset step= 0.324 pixels
 Auto-refine: Iteration= 16
 Auto-refine: Resolution= 7.58571 (no gain for 1 iter) 
 Auto-refine: Changes in angles= 0.576587 degrees; and in offsets= 0.162731 pixels (no gain for 1 iter) 
 Auto-refine: Estimated accuracy angles= 0.749 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= 17
 Auto-refine: Resolution= 7.58571 (no gain for 1 iter) 
 Auto-refine: Changes in angles= 0.357051 degrees; and in offsets= 0.127339 pixels (no gain for 0 iter) 
 Auto-refine: Estimated accuracy angles= 0.704 degrees; offsets= 0.392 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.58571 (no gain for 2 iter) 
 Auto-refine: Changes in angles= 0.269964 degrees; and in offsets= 0.096948 pixels (no gain for 1 iter) 
 Auto-refine: Estimated accuracy angles= 0.7 degrees; offsets= 0.39 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.58571 (no gain for 3 iter) 
 Auto-refine: Changes in angles= 0.277262 degrees; and in offsets= 0.137002 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.702 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/job058/run_class001.mrc
 Auto-refine: + Final model parameters are stored in: Refine3D/job058/run_model.star
 Auto-refine: + Final data parameters are stored in: Refine3D/job058/run_data.star
 Auto-refine: + Final resolution (without masking) is: 3.88537
 Auto-refine: + But you may want to run relion_postprocess to mask the unfil.mrc maps and calculate a higher resolution FSC


iter2のChanges in anglesは、正しいhandednessの参照像を用いた場合(Refine3D/first3dref)は 19.3751度でしたが、上記では 48.3435 度となっており、割と激しく角度が動いているのかなという気はします。が、それ以降の経過はどちらのhandednessの参照像を使った場合でも概ね同じであり、むしろ違うhandednessでやったほうが分解能が少し上がってしまいました(左巻きαヘリックスも確認できます)。

最適化後のparticle starファイルである run_data.star を比較すると、rlnAnglePsiが180度反転しており、他はほぼ同じであると確認できました。

ということで、今回の場合は、すでに最適化済みのparticle starファイルを入力として、参照像のhandednessを反転してからrefinementすると、反転後のhandednessで3D再構成されると確認できました(単粒子の面内回転psiは180度反転)。最初の方のイテレーションで投影角度の全探索が行われる(ですよね?)のでこういうことが可能なのだと思います。

全探索に任せるのが気持ち悪ければ、参照像のhandを反転することに加え、入力単粒子starファイルのrlnAnglePsiを [-180, 180]の範囲で180度反転させればいいと思います。

なにか勘違いしていればご指摘ください。

11 Wrapping up (終わりに)

11.1 Making a flowchart(フローチャートの作成)

どのような手順を経て最終的な3D再構成に至ったか、整理したいでしょう。Finished jobsリストから最後に実行したジョブを選択し、【Job actions】ボタンで"Making flowhart"オプションを実行すると、そこへ至るまでのフローチャートを出力できます(PDF)。これが機能するためにはシステムに LaTeXおよびTikZパッケージがインストールされている必要があります。

※Polish/polishの場合について示します。

rln-flowchart-01.png


最初のページには、ジョブ番号やエイリアス名を表示していない概要フローチャートがあります。論文へ載せる場合などに役に立つと思います(好みのベクトルデータ描画ソフトとで編集したりした後で)。

rln-flowchart-02.png


概要フローチャートの後には、選択したジョブへ至るまでの詳細フローチャートがあります。10ステップよりも長いフローチャートの場合は適宜分割して表示されます。また、ワークフローの中にはブランチがあると思います(※紫背景でジョブタイプ名が赤枠で囲ってあるブロック)。後のページにそれぞれのブランチの詳細ワークフローが続きます。ブランチのジョブタイプ名をクリックすると、そのブランチのワークフローのページまで飛べます。

rln-flowchart-03.png


※私の環境では、必要なパッケージはインストールされていましたが、なぜかMaking flowchartするとRELION GUIがフリーズしてフローチャートが出力されませんでした。latexコマンドでflowchart.dviを生成するところまではうまく行っていたので、dvipdf flowchart.dvi とコマンドを手打ちしたところ、flowchart.pdf を生成できました。

11.2 Cleaning up your directories (ディレクトリのお掃除)

ディスク領域を節約するため、RELIONには終了したジョブのディレクトリをお掃除する機能が備わっています。お掃除には2つのモードがあります。

  • 'gentle'
    • 選択したジョブのディレクトリ内の中間ファイルのみ削除します。
    • 例えば2D classification, 3D classification, 3D auto-refineジョブの場合、途中のイテレーションのデータのみが削除されます。
  • 'harsh'
    • 選択したジョブのディレクトリ内で、そのジョブに続けて次のジョブを実行する際に必要になるような重要なファイルも含めて削除します。
    • 例えばMotionCorrジョブの場合はドリフト補正後の平均画像、Particle extractionジョブの場合は切り出された単粒子画像のスタックファイルといった具合に、そのジョブの成果物ファイルも含めて消されます。

【Job actions】ボタンを使って個々のジョブを掃除することもできますし、上部メニューバーのJobsプルダウンメニューから 'Gently clean all jobs'または'Harshly clean all jobs'を選択してすべてのジョブを一括でお掃除することもできます。RELIONチュートリアルとともに配布されている計算済みデータを格納したプロジェクトディレクトリのtarボール(※relion30_tutorial_precalculated_results.tar.gz)は、'Gently clean all jobs'によるお掃除後に作成されたものです。一通りプロジェクトが終了した後で、データを長期間とっておきたいといった場合には、'Gently clean all jobs'を検討してもいいでしょう。

11.3 Asking questions and citing us

以上でチュートリアルは終了です!このチュートリアルが楽しく役に立つものであったことを願います。RELIONに関して何が質問があれば、まずはRELION WikiのFAQとCCPEMのメーリングリストをチェックしてください。もし答えが見つからない場合は、CCPEMメーリングリストに質問を投稿してください。そしてこれは心からのお願いですが、RELIONに関して直接Sjorsにメールするのはやめてください。メールが来てもSjorsはお返事できません。

RELIONがあなたの研究のお役に立った場合は、我々の論文を引用してください。また、RELIONが役に立つものであったことを周りの人にも教えてあげてください。

11.4 Further reading (RELIONをもっと深く知るために)

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.


※お疲れ様でした!

rln-final-map.png

12 Appendix A: notes on installation

13 Appendix B: using RELION