fMRI預處理實戰(zhàn):DPABI+SPM12參數(shù)設置與避坑指南)
先放個結論任務態(tài)fMRI預處理沒有必要去專門找一套“任務態(tài)專用工具包”DPABI SPM12 這套組合完全夠用而且大部分核心算法本來就是 SPM12 提供的DPABI 更像是一個把參數(shù)串起來、批量跑并幫你省掉大量手寫 batch 的調度界面。真正容易出問題的從來不是軟件本身而是任務態(tài)數(shù)據(jù)在參數(shù)選擇上跟靜息態(tài)的那幾個微妙差異。我最早自己跑任務態(tài)數(shù)據(jù)時也踩過坑按照靜息態(tài)的預處理習慣順手勾了帶通濾波、順手把全局信號回歸也打開了結果一級分析里任務條件的主效應怎么都不顯著后來排查半天發(fā)現(xiàn)是預處理把任務相關的低頻成分濾掉了。那次之后我把任務態(tài)預處理的每個環(huán)節(jié)重新捋了一遍這篇筆記就是把“為什么這么填”和“實操中怎么避坑”一起整理出來主要面向第一次用 DPABI 處理任務態(tài)數(shù)據(jù)、又不想在預處理階段埋雷的同行。1. 為什么任務態(tài)數(shù)據(jù)不能直接照搬靜息態(tài)流程1.1 任務態(tài)與靜息態(tài)在信號假設上的根本差異靜息態(tài)數(shù)據(jù)關注的是被試在無任務狀態(tài)下自發(fā)的 BOLD 低頻漲落默認的分析窗口是 0.01~0.1 Hz對應的是幾十秒到上百秒的慢波成分。所以靜息態(tài)預處理里常常要做帶通濾波把這個頻段之外的心跳、呼吸漂移、掃描儀慢漂移盡量過濾掉方便后面算 ALFF、ReHo、功能連接這些指標。任務態(tài)的邏輯完全不同。任務態(tài)關心的是 BOLD 信號是否跟著刺激或任務條件同步變化事件相關設計里甚至要求分辨出一兩秒內出現(xiàn)的峰值變化這是相對高頻的信息。如果你把靜息態(tài)那套 0.01~0.1 Hz 帶通直接套過來等于把低頻端砍掉——當一個任務的組塊節(jié)律較長時任務相關成分本身就在 0.01 Hz 以下濾掉之后就什么都看不出來了。這也是我在開頭說的那次翻車經(jīng)歷。所以任務態(tài)預處理里通常不主動做帶通濾波最多做一下高通去漂移低頻慢漂移部分交給后面一級分析SPM 里的 high-pass filter去處理。這個思維轉變是所有參數(shù)調整的前提。1.2 DPABI與SPM12的分工誰負責調度誰負責算賬很多初學者打開 DPABI 之后會有點困惑為什么這個界面里也能看到 SPM12 的東西因為 DPABI 本質上是構建在 SPM12 之上的批處理工具SPM12 負責具體算法比如 DICOM 導入、Slice Timing、Realign、Coregister、Segment、Normalize、Smooth 這些底層操作DPABI 負責把操作串成管道、給每個被試批量執(zhí)行、輸出便于核查的目錄結構。環(huán)節(jié)主要負責者圖形界面與流程編排DPABI批量循環(huán)、數(shù)據(jù)目錄組織DPABIDICOM 轉 NIFTISPM12 / dcm2nii 接口時間層校正、頭動校正、配準、分割、標準化、平滑SPM12預處理質量指標生成DPABI一句話總結你把樣本信息、TR、層數(shù)、層順序、平滑核大小這些參數(shù)填在 DPABI 面板上DPABI 生成對應的 SPM12 batch再用 MATLAB 指令逐被試執(zhí)行。所以如果 SPM12 路徑?jīng)]配好或者版本不匹配DPABI 界面再正常也沒用。1.3 任務態(tài)預處理流程總覽與文件命名約定任務態(tài)預處理的基本鏈路如下DICOM 轉 NIFTI如果數(shù)據(jù)采集端已經(jīng)是 NIFTI 則跳過去除前幾個時間點可選推薦時間層校正 Slice Timing頭動校正 Realign功能像與結構像配準 Coregister結構像分割 Segmentation標準化 Normalize到 MNI 空間平滑 SmoothSPM12 生成文件時有非常固定的前綴規(guī)則這個一定要記住否則你在結果目錄里看到一堆亂碼一樣的文件名會直接懵前綴含義a時間層校正后r頭動校正realigned后w標準化normalized后s平滑smoothed后c1/c2/c3分割出的灰質/白質/腦脊液rc1/rc2/rc3重采樣到 MNI 空間的灰質/白質/腦脊液rp_頭動參數(shù)文件每幀 6 列X/Y/Z平移 roll/pitch/yawmean*頭動校正后的平均功能像如果整個鏈路完整跑完最后用于統(tǒng)計的文件通常形如swar*.nii先做了 slice timinga再做 realignr再 normalizew最后 smooths。通過文件前綴就能快速反推這個被試跑到哪一步排查問題會快很多。2. 數(shù)據(jù)整理與環(huán)境配置先把地基打牢2.1 目錄結構與命名規(guī)范DPABI 對數(shù)據(jù)目錄有一個隱含要求每個被試一個獨立文件夾T1 和功能像分別放好。我推薦在開始之前就固定好目錄結構不要邊跑邊改。我自己的習慣是這樣的Project/ ├── Participants/ │ ├── 001/ │ │ ├── T1/ │ │ │ └── *.dcm │ │ └── Func/ │ │ └── *.dcm │ ├── 002/ │ │ ├── T1/ │ │ │ └── *.dcm │ │ └── Func/ │ │ └── *.dcm │ └── ...幾個容易踩的細節(jié)被試編號不要用純中文、不要帶空格用001、002這種定寬數(shù)字避免后面排序變成 1、10、2 的詭異順序。如果每個被試有多段任務多個 run建議在 Func 目錄下再拆分子目錄比如Func_Run1、Func_Run2或者按 DPABI 支持的多功能層目錄去組織。DICOM 文件來源檢查一下同一個 DICOM 文件夾里如果混入了定位像或場圖序列DPABI 的 DICOM 導入階段會分錯文件預處理出來的時間點數(shù)對不上。2.2 MATLAB、SPM12、DPABI的版本與路徑配置這三個軟件的版本搭配很關鍵。SPM12 官方支持 MATLAB 的時間范圍比較寬但太新的 MATLAB 偶爾會出現(xiàn)圖形句柄或兼容性問題我自己的經(jīng)驗是 MATLAB R2020b 到 R2022b 這段區(qū)間比較穩(wěn)妥。DPABI 建議用最新版 V6.x老版本在任務態(tài)數(shù)據(jù)處理上的選項不如新版本完整。在 MATLAB 里配置路徑的方式addpath(genpath(/path/to/spm12)); addpath(genpath(/path/to/DPABI_V6.x)); savepath;注意兩點先加 SPM12再加 DPABI。因為 DPABI 會調用 SPM12 的函數(shù)如果路徑順序反了某些同名函數(shù)可能被覆蓋跑出來的行為會非常詭異。配置完路徑后最好重啟一次 MATLAB再在命令行里輸入spm和dpabi驗證能否正常啟動。如果spm打開窗口正常、dpabi也能打開主界面說明環(huán)境基本沒問題。2.3 去除前幾個時間點任務前靜息與dummy scan怎么處理fMRI 掃描開始的幾個 volume 往往信號極不穩(wěn)定因為磁化狀態(tài)還沒有達到穩(wěn)態(tài)T1 對比度會異常偏高。所有預處理流程里我都會先切掉頭幾個時間點再進入正式步驟DPABI 里也有對應的選項。切幾個常見做法是 4~10 個。這個取決于掃描儀和序列穩(wěn)妥一點的做法是先看原始數(shù)據(jù)的時序均值圖像如果前幾幀明顯比后面亮或出現(xiàn)整體強度漂移就切掉 10 個如果本身穩(wěn)定切 4 個也可以。切太多會損失數(shù)據(jù)尤其是實驗本身時長較短時。有一個“雙重切除”的坑如果采集時掃描序列已經(jīng)設置了 dummy scan即前 n 幀采集但不保存那接收到的 NIFTI 數(shù)據(jù)可能已經(jīng)自動少了那幾個 volume你又在 DPABI 里切 10 幀等于白白損失數(shù)據(jù)。所以開始前一定問清楚掃描參數(shù)里有沒有 dummy scan再決定 DPABI 里要不要切、切多少。3. 預處理核心參數(shù)逐項拆解每一站該怎么填3.1 時間層校正參考層和slice order的確定方法功能像是一個 volume 一個 volume 拍的但每個 volume 內部并不是同一瞬間完成的掃描儀在 TR 時間內逐層采集所以每一層對應的真實采集時間有一個偏移。時間層校正的作用就是通過插值把同一 volume 所有層的信號對齊到某個參考時間點上。對任務態(tài)來說這一步比靜息態(tài)更敏感。事件相關設計中 HRF 的峰值本來就短幾百毫秒的層間偏移如果不校正不同腦區(qū)的 BOLD 響應看起來會像在不同時間點到達時間精度不夠的后果就是激活檢測效率下降。需要填寫兩個關鍵參數(shù)slice number 和 slice order。slice number直接等于你掃描序列的層數(shù)比如 32 層就填 32。slice order指的是每一層實際被采集的先后順序這個必須從掃描序列參數(shù)里確認不要猜。常見的三種輸入形式掃描方式Slice Order 填寫示例32層說明連續(xù)升序從上到下或從下到上1:32或32:-1:1飛利浦、GE 常見隔層先奇數(shù)后偶數(shù)1:2:31, 2:2:32西門子常見隔層先偶數(shù)后奇數(shù)2:2:32, 1:2:31需從掃描參數(shù)確認參考層reference slice一般取全部層的中位時間點那一層比如 32 層通常選第 16 層。如果某個興趣區(qū)在腦頂部也可以考慮把參考層往頂部區(qū)域靠一下盡量減少該區(qū)域相對參考層的時間偏移。3.2 頭動校正任務誘發(fā)的微動作如何記錄和評估頭動校正是用剛體變換把每一幀對齊到參考幀通常是第一幀或平均幀估計 6 個參數(shù)X/Y/Z 平移以及繞三軸的旋轉。SPM12 的 realign 會輸出rp_*.txt文件里面每一行對應一個 frame 的 6 個參數(shù)后續(xù)分析里經(jīng)常被當作協(xié)變量回歸進模型。判斷頭動是否嚴重不能只看平移還要看旋轉。經(jīng)驗做法是計算 frame-wise displacementFDPower 2012 年提出的公式大致是FD(i) |dX(i)| |dY(i)| |dZ(i)| |alpha(i)| * 50 |beta(i)| * 50 |gamma(i)| * 50旋轉角度單位是弧度乘以 50 近似等于頭部半徑約 50 mm處的弧長位移。我的實操閾值參考指標輕度需要警惕通常建議剔除最大平移 1.5 mm1.5~3 mm 3 mm最大旋轉 1.5°1.5°~3° 3°平均 FD 0.2 mm0.2~0.5 mm 0.5 mm壞幀比例FD 0.5 mm 10%10%~20% 20%任務態(tài)數(shù)據(jù)頭動問題比靜息態(tài)更突出被試在任務中要做按鍵、說話、咀嚼等動作頭動不僅幅度大而且往往和任務條件時間鎖相關。參數(shù)填得好不好直接影響后面一級分析是否出現(xiàn)任務相關偽影。3.3 配準與分割為什么先配準后標準化順序不能反預處理的后半段思路是先把功能像和結構像T1對到同一個空間再用 T1 分割得到從個體空間到 MNI 標準空間的變形場最后把這個變形場施加到功能像上。注意順序不能反原因很簡單如果你先單獨分割 T1 并把它標準化到 MNI再反過來把功能像配準到已經(jīng)標準化的 T1就會有兩次獨立的插值誤差功能像會變得很“糊”組分析的時候敏感度下降。DPABI 里配準一般選擇 T1 到功能像的平均像做 coregistration然后用平均像目檢。分割用的是 SPM12 的 New Segment會生成 c1、c2、c3 三張圖對應灰質、白質、腦脊液。配準做完一定要看質量。SPM12 里可以用 Check Reg 同時顯示 T1 和功能平均像肉眼檢查邊緣是否對齊如果發(fā)現(xiàn)明顯錯位不要急著繼續(xù)標準化先回到原始影像看看是不是方向標簽有問題必要時手動 reorient。3.4 標準化到MNI空間普通標準化與DARTEL怎么選標準化是把不同被試的大腦形態(tài)統(tǒng)一到同一模板空間這樣后續(xù)組分析才能逐體素比較。DPABI 里提供普通標準化和基于 DARTEL 的標準化兩個方向。普通標準化就是直接用 SPM12 的 unified segmentation 得到變形場再寫到功能像上。速度快、默認參數(shù)成熟絕大多數(shù)任務態(tài)研究完全夠用。DARTEL 的做法是先在樣本內部生成一個平均模板再用這個模板對每個被試做更精細的配準精度更高但計算量大很多而且對單個被試的數(shù)據(jù)質量更敏感。小規(guī)模精細空間歸一或者樣本腦形態(tài)差異偏大時可以選 DARTEL常規(guī)任務態(tài)全腦分析直接普通標準化即可。方式速度精度適用場景普通標準化快足夠常規(guī)任務態(tài)組分析DARTEL慢更優(yōu)VBM、樣本間結構差異大、需要精細配準時標準化時還要確定體素大小。我一般設置 2×2×2 mm 或保持原始分辨率也有很多人習慣 3×3×3 mm計算更快、后續(xù)平滑后統(tǒng)計更容易滿足隨機場假設。如果被試樣本量不大我建議 2 mm保留更多空間細節(jié)。3.5 平滑核大小與任務激活檢測的關系平滑是一個容易被低估的步驟。高斯平滑核FWHM的作用不只是把圖變糊它有三個實際收益提高信噪比減少高頻噪聲對統(tǒng)計量的影響讓不同被試之間功能激活位置稍微偏移時不至于完全配不上使誤差場更接近高斯隨機場模型滿足基于 GRF 的體素水平校正假設。FWHM 的選擇核心看數(shù)據(jù)分辨率和興趣區(qū)大小。常見原則是取原始體素大小的 2 倍左右例如 2 mm 體素用 6 mm FWHM3 mm 體素用 6 mm 也常見。如果分析目標是非常小的核團比如杏仁核、導水管周圍灰質FWHM 太大很容易把小激活抹平建議用 4 mm 甚至更小。反之如果只做全腦大規(guī)模激活區(qū)8 mm 也不夸張。DPABI 里平滑核的填寫格式是三個數(shù)字比如[6 6 6]分別代表 X/Y/Z 方向。4. 任務態(tài)預處理中必須避開的幾個坑4.1 頭動與任務條件相混淆時該怎么辦這是一個非常隱蔽但傷害極大的坑。假設你的實驗組塊是“任務 30 秒 休息 30 秒”被試在任務期做高頻按鍵頭動也跟著任務期同步增加。這種頭動和任務條件高度相關會產(chǎn)生兩個后果一是頭動參數(shù)回歸進 GLM 后任務條件本身的一部分真實差異也被吸走激活檢測效力下降二是如果頭動與任務完全同步SPM 的模型有可能分不清到底是 BOLD 激活還是運動偽影。判斷方法很簡單畫一條頭動曲線再畫一條任務開/關方塊圖肉眼對比兩者是否同步更嚴謹一點可以算一下每個幀的頭動 FD 和任務條件指示變量之間的相關系數(shù)。如果相關系數(shù)很高這個被試的數(shù)據(jù)就要特別小心。處理手段按嚴重程度排序輕度同步頭動可以在 GLM 中添加額外的微運動回歸器或者用 F 檢驗對比有無頭動協(xié)變量時的任務效應是否穩(wěn)定中度同步頭動可以考慮按壞幀剔除scrubbing或僅剔除某幾個幀嚴重同步頭動直接剔除被試可能是更理智的選擇因為模型已經(jīng)很難分離真實激活和運動偽影。4.2 帶通濾波在任務態(tài)中的風險前面提過帶通濾波是靜息態(tài)預處理的常規(guī)操作但任務態(tài)不建議加。舉個具體例子一個長組塊設計任務塊 120 秒、控制塊 120 秒整個周期是 240 秒對應的頻率是 1/240 ≈ 0.0042 Hz。如果你用 0.01~0.1 Hz 的帶通這個任務相關頻率就被徹底濾掉了。即使任務塊沒那么長只要任務節(jié)律偏低頻帶通都會造成不同程度的信息損傷。正確做法是預處理階段不濾波最多做線性去漂移或高通低頻漂移的抑制交給 SPM12 的一級分析模塊SPM 的 high-pass filter 默認 128 秒截止你也可以根據(jù)任務設計改成 256 秒甚至更長。如果在 DPABI 里看到濾波選項確認它是“僅去線性漂移”還是“帶通濾波”任務態(tài)數(shù)據(jù)一律選不濾波或只做去線性漂移。4.3 多run數(shù)據(jù)如何整理與預處理一個被試跑了兩個 run預處理時不要圖省事把兩個 run 拼成一個四維文件。正確做法是讓 DPABI 分別對每個 run 做時間層校正和頭動校正之后標準化時共用同一個 T1 變形場。目錄組織上我建議001/ ├── T1/ │ └── *.nii ├── Func_Run1/ │ └── *.nii └── Func_Run2/ └── *.nii在 DPABI 里可以添加多個功能像目錄DPABI 會逐個處理。之后在做 SPM12 一級分析時兩個 run 作為兩個 session 輸入PD 模型里設置相同的回歸量即可。如果兩個 run 之間的 TR、層數(shù)、層順序不一樣不常見但換序列時會出現(xiàn)就必須當作兩批數(shù)據(jù)分別設置預處理參數(shù)不能一個 settings 全跑。4.4 預處理結果質檢清單預處理跑完不是結束質檢必須做。我的習慣是建立一個簡單的檢查清單每個被試過一遍時間點數(shù)是否符合預期DICOM 轉換后是否有缺失頭動曲線是否有突然的大尖峰T1 與功能平均像的配準邊緣是否吻合標準化后的功能像在 MNI 模板上看起來是否“正?!庇袥]有明顯扭曲平滑后的圖像有沒有偽影或強度異常rp_*.txt 文件是否完整幀數(shù)和原始 volume 數(shù)一致。每次跑完預處理后我都會把每個被試的質檢結果記錄在一個表格里標記 PASS / Warning / Fail再決定哪些數(shù)據(jù)進入后續(xù)統(tǒng)計分析。這一步花不了多少時間但能省掉后面大量返工。5. 我在實操中碰到的具體報錯與排查過程5.1 路徑和版本導致的MATLAB報錯最常見的報錯是打開 DPABI 時提示找不到 SPM 函數(shù)或者運行時報Undefined function or variable spm。這個基本就是 SPM12 路徑?jīng)]加上去重新確認addpath之后重啟 MATLAB 就好了。另一個比較隱蔽的是 MATLAB 版本過新導致的圖形窗口報錯例如Error using matlab.ui.Figure或者某些 SPM 的 GUI 函數(shù)異常。SPM12 本身一直有更新但如果你手頭的 SPM12 是幾年前的舊版配 MATLAB R2023b 以上就可能出問題。排查方向很簡單先看 SPM12 的版本日期再看 MATLAB 版本兩個跨度太大就換掉其中一個。我個人的經(jīng)驗值是 MATLAB R2020b 最新 SPM12 最新 DPABI 最穩(wěn)。5.2 配準結果錯位的處理有次我給一個數(shù)據(jù)集跑完預處理QC 階段發(fā)現(xiàn)某個被試的功能像和 T1 明顯錯位腦組織邊緣是錯開的功能像整體相對 T1 旋轉了一點。查來查去問題出在這個被試采集時的頭位明顯偏轉DICOM 轉 NIFTI 后圖像的齊次坐標和 T1 不一致。這種情況下強制做 coregistration 很容易配歪。解決思路是在正式預處理前加一步手動 reorient在 SPM12 的 Display 里打開功能像和 T1用“Reorient matrix”把兩者大致擺到同一朝向再重新跑預處理。DPABI 里也提供 Reorient 相關選項可以在預處理前對所有被試統(tǒng)一設置 AC-PC 對齊。這個坑很難通過調參數(shù)解決因為根源在原始數(shù)據(jù)方向信息。所以我的習慣是預處理前先隨機抽 2~3 個被試用 SPM12 打開圖像看看方向是否正確再決定跑全量。5.3 頭動過大的被試是剔除還是保留遇到頭動大的被試很多人會糾結到底剔除還是保留。我的決策流程是先看這個被試是否全 run 頭動都大還是只有某個 run 頭動大看頭動是否和任務條件同步看核心 ROI 區(qū)域是否出現(xiàn)明顯偽影或信號異常用含頭動協(xié)變量的模型跑一次對比結果穩(wěn)定性。如果只是某個 run 頭動大我會選擇只剔除該 run 而不是整個被試保留另一個合格 run 的影像數(shù)據(jù)。如果頭動普遍超標且與任務同步我會直接剔除哪怕樣本量變小也不硬留。因為任務態(tài) GLM 對運動偽影的容忍度比靜息態(tài)更低留一個嚴重頭動被試進組分析可能產(chǎn)生假激活也會掩蓋真實效應得不償失。5.4 預處理完成后給下一步留好哪些文件預處理完成并不代表萬事大吉最終進入一級分析前要確保以下文件都完整每個被試標準化并平滑后的功能像文件swar*.nii或類似前綴每個 run 的頭動參數(shù)文件rp_*.txt后續(xù) GLM 要作為協(xié)變量每個 run 對應的幀數(shù)信息便于檢查模型自由度如果做了壞幀剔除保留剔除幀編號列表方便后續(xù)寫模型時做 missing frame 處理。這些文件我會統(tǒng)一在Preprocessed目錄下按被試歸檔和原始數(shù)據(jù)分開。這樣一來下一步在 SPM12 里建一級模型時只需要按被試填功能像路徑和rp_文件路徑就能順暢地接上。我自己跑任務態(tài)數(shù)據(jù)到現(xiàn)在最深的感受是預處理這部分真正決定分析質量的往往不是軟件多復雜而是你有沒有理解每一步的參數(shù)從哪來、為什么要這么設。頭動、時間層校正、濾波設置這幾個環(huán)節(jié)尤其值得多花時間核對。如果你也準備用 DPABI 跑任務態(tài)預處理我建議第一次先抽 3 個被試把整條鏈路跑通、把所有中間結果都檢查一遍再大批量跑這樣會比直接懟全樣本數(shù)據(jù)穩(wěn)妥得多。