Human Connectome Projectが提供しているpipelineを用いた解析も行っています。
私はMacで解析を行っているので、このページの情報はすべてMacでのものです。この内容を更新しているときのHCP pipelineのバージョンは4.1.3です。頻繁にバージョンアップをするので、内容が古くなっている可能性があることをご承知おきください。
HCP pipelineはシェルスクリプト群でできていて、Githubからダウンロードできます。基本的には、Exampleフォルダにあるスクリプトを自分の研究データとコンピュータ環境に合わせて変更して使用します。そのほとんどはパスの設定です。相対パス指定はせずに絶対パス指定にします。MacでHCP pipelineを使う場合は、ターミナルにフルディスクアクセスを割り当てて、システム整合性保護(System Integrity Protection, SIP)を無効化します。
HCP pipelineを動かすには、バックグラウンドで動作してくれる下記のソフトをインストールしておく必要があります。
FSL本体 バージョン6.0.2以降
Freesufer バージョン6.0 (7.0ではまだ動作しません)
Connectome workbench バージョン1.4.2
MSM バージョン3.0
FSL FIX バージョン1.06.12以降
Matlab 2017b以降とsignal processing toolbox(Matlabがない場合はruntime 2017bも利用可)
R バージョン3.6.3
※FSLを最新版にアップデートするには、以下のコマンドをターミナルで実行します。
> python2.7 /usr/local/fsl/etc/fslinstaller.py -c
※HCP pipelineの一部(tkregister2を使う、FreeSurferPipelineとGenericfMRIVolumeProcessingPipeline)はそのままではCatalinaで動作しません。Mojaveなどの少し前のOSで実行するか、tkregister2をtkregister2_cmdlに書き換えます。
撮像について
HCP pipelineを使ってデータ解析をしようと決めて、課題の準備もでき、いざ撮像プロトコルを決めようとしたとき、どのような撮像をすればよいのか分からないということが考えられます。そのような場合のために、ここではどのような撮像が最低限必要なのかを説明していきます。MRIはSiemens Prismaで、Multibandを用いて撮像することを前提としています。
まず必要になるのが、T1強調画像(T1w)、T2強調画像(T2w)、Fieldmapです。T1wとT2wの解像度をともに0.8mm isovoxelとして、同じ撮像範囲を設定します。FieldmapはGradient EchoとSpin Echoの2種類から選べますが、時間が短くて済むのでSpin Echo (SE)がおすすめです(GREは約2分で、SEは約20秒)。Fieldmapの内容については下記を参照してください。
課題中の脳活動を計測するBOLD EPI画像の取得では、Multibandによる撮像を行います。解像度を2 mm x 2 mm x 2 mmとして、Multiband factor = 6に設定し、撮像枚数は72 slicesとします。in-plane FOVは96 x 96です。撮像範囲は、192 mm x 192 mm x 144 mmとなります。TRをできるだけ短く、1秒以下にすることをおすすめします。Flip angleはTRに対応するErnst角に設定します。例えば、TR = 1 sならFA = 59°、TR = 0.9 sならFA = 56°くらいです。パラレルイメージング(GRAPPA)は設定しません。TEは30 ms(30を超えてしまう場合は最小値)、Partial Fourierはオフにします(ただし、OFCをターゲットにしたい場合など、TEを短くする実験の場合は入れたほうが良い)。MultibandのSpecialカードで、Single-band imagesとMB LeakBlock kernelのオプションをonにします。エンコーディング方向をPAで撮像するのですが、Routineの撮像方向はAPにしておきます。SpecialカードでInvert RO/PE polarityをonにして、エンコーディング方向をPAに設定します。
課題の各runとT1wの直前にSE Fieldmapを取得します。課題の脳活動を計測するEPI画像と同じ解像度、同じEcho spacing(Bandwidth)、同じ撮像範囲で行います。エンコーディング方向がAPの画像を1 volume、PAの画像を1 volume撮像します。ここでも、MultibandのSpecialカードで、MB LeakBlock kernelのオプションをonにします。エンコーディング方向をPAで撮像する場合は、Routineの撮像方向をAPにしておき、SpecialカードでInvert RO/PE polarityをonにします。
課題遂行時のEPI撮像の他に、Resting-state(安静時)のEPI撮像を行うこともできます。データ解析時の位置合わせに使うのと、Resting-stateでの機能的結合の解析に使えます。
データ解析
以下では、小川の行っているデータ解析の処理内容に沿って、どのようにしてHCP pipelineを使うかを説明していきます。
準備
取得した画像のDICOMデータを、dcm2niixを用いてnii.gz (compressed nifti形式)に変換します。
SetUpHCPPipeline.shの内容を、自分の環境に合わせて書き換えます。このスクリプトで行うのは主にパスの設定です。
構造画像の処理
HCP pipelineでは、脳構造を調べるにも脳活動を調べるにも、まずは構造画像の処理が必要になります。
必要なファイル:
・T1w画像
・T2w画像
・Spin Echo FieldMap (エンコーディング方向がPAとAP) ※GRE FieldMap (MagnitudeとPhase)も可
以上の計4つのファイルが必要になります。*T1w画像のみでHCP pipelineを使いたい場合は下記を参照してください。
T1w画像とT2w画像は解像度・FOV・撮像位置が一致していることが前提です。
FieldMapの画像はT1w 画像の撮像範囲と一致していなくても良いです。(というよりも、一致させられないでしょう)
実行するシェルスクリプト
1.PreFreeSurferPipelineBatch
2.FreeSurferPipelineBatch
3.PostFreeSurferPipelineBatch
機能画像の処理
Preprocessingに当たる処理を行います。
必要なファイル
・GRE EPI画像
・上記のSigleband reference画像(SBRef)
・Spin Echo FieldMap (エンコーディング方向がPAとAP)
各Runで以上の4つのファイルが必要になりますので、SpecialカードでSBRefを取得するように設定し、すべてのRunにおいて脳活動のEPI撮像の直前にSE Fieldmapを撮像しておくのが良いです。ただ、SE Fieldmapは必ずしも毎run行わなくても良いでしょう。
・Topupのためのパラメタファイル(acqparams.txt)とコンフィグファイル(b02b0.cnf)
歪み補正を行うTopup処理のパラメタファイルの第4パラメタは、取得する画像に応じて計算する必要があります。SE FieldmapをNIFTIに変換するときにできるjsonファイルを見ると、EffectiveEchoSpacingとEchoTrainLengthという項目があります。この2つの数値の積が第4パラメタになります。
スライス枚数が偶数の場合はデフォルトのコンフィグファイルで良いのですが、スライス枚数が奇数の場合はsubsampオプションを以下のように書き換えます。
--subsamp=1,1,1,1,1,1,1,1,1
実行するシェルスクリプトは以下のとおりです。
1.GenericfMRIVolumeProcessingPipelineBatch
2.GenericfMRISurfaceProcessingPipelineBatch(ここではResting-stateのみ)
脳活動の計算(Resting-state)
頭部の動き、白質、脳脊髄液などに由来するノイズ等を取り除くため、ICA+FIXの処理を行います。
続いてMSMAllの処理を行い、MSMAllの後処理に当たるDeDriftの処理を行います。
実行するシェルスクリプトは以下のとおりです。
1.IcaFixProcessingBatch
2.MSMAllPipelineBatch
3.DeDriftAndResamplePipelineBatch
脳活動の計算(Task)
HCP pipelineでは、FSLでTaskの脳活動解析を行います。SPMで解析をしたい方は、GenericfMRIVolumeProcssingPipelineBatchの出力を解析すると良いでしょう。そのときは、SPMとFSLで標準脳の位置表示がずれることがあることに注意します。
残りはTBA
実行するシェルスクリプト
1.GenericfMRISurfaceProcessingPipelineBatch(Registrationの対象をMSMAllにする)
2.PrepareForTaskAnalysis (なぜか、3T_Size_Checkingフォルダに入っている)
2-1.copy_evs_into_results
2-2.generate_level1_fsf
3.TaskfMRIAnalysisBatch
T1w画像のみでHCP pipelineを使う方法
Freesurferに関係する3つのスクリプトに、以下のオプションを追加します。また、T2w画像やFieldmap画像などの設定をNONEにします。
--processing-mode="LegacyStyleData"
機能画像の処理においてSlice timing correctionを行う方法
TRが長い場合はSlice timing correction(スライスタイミング補正)を行う方が良い場合があります。また、Dynamic Causal Modeling (DCM)を行うにはSlice timing correctionが必須となります。Slice timing correctionを行うには、機能画像のスクリプトに以下のオプションを追加します。
--processing-mode="LegacyStyleData"
--doslicetime="TRUE"
--slicetimerparams="--tcustom=$filename"
$filenameには、スライスタイミングを記録したファイル名(テキストファイル)を入れます。そのファイルには、1行に1数値を入れておきます(N行目がNスライス目に相当)。単位はTRなので、0から1までの数値となります。
スクリプトを実行しようとすると、/bin/bash: bad interpreter: Operation not permittedというエラーが出る場合の対処方法
拡張ファイル属性にcom.apple.quarantineが存在するため実行できないようになっています。以下のようにして、この属性を削除すると実行できるようになります。
xattr -d com.apple.quarantine filename