【初心者向け】基礎&実践プログラミング

初心者がつまづきやすいところ、最短で実力が身につく方法をお伝えします。

【主成分分析】FLICAを用いたマルチモーダル主成分分析

目的

  • FLICAを用いた画像のマルチモーダル主成分分析(ICA)

リポジトリ

FLICAに必要なスクリプトを次のコマンドで取得。今のところプライベートリポジトリにしています。ご興味のある方はコメントしてください。

git clone https://github.com/Hexans/FLICA_corrected.git

必要なデータ

必要なデータは、モダリティごとに全ての被験者の画像がまとめられた4D画像である。例えば、次のようなファイルを用いる。

  • VBM: GM_mod_merg_s?.nii.gz
  • TBSS: all??skeletonised.nii.gz
  • FreeSurfer: ?h.thick.fsaverage.mgh and ?h.pial.area.fsaverage.mgh
  • etc.

さらに、共変量(年齢・疾患グループ・性別)との相関もみることができる。共変量ごとに1列に数値をまとめテキストファイル(.txt)として保存する。

  • 年齢: age.txt
  • 疾患グループ: group.txt
  • 性別: sex.txt

実行

FLICAを実行するには、次のコマンドを実行。

ICAのComponent数を被験者数の4分の1以下(Components < subjects/4)に設定し、Overfittingを避ける。また、被験者ごとのノイズ推定を省く(opts.lambda_dims = '';)ことでもOverfittingを避けることができる。

run_FLICA.sh

%% Set up:
addpath([getenv('FSLDIR') '/etc/matlab/'])
addpath('~/Documents/software/flica')

%% Load data
Yfiles = {'all_FA_skeletonised.nii.gz',
    'all_MD_skeletonised.nii.gz',
    'all_AD_skeletonised.nii.gz',
    'all_RD_skeletonised.nii.gz'};
% NOTE that these should be downsampled to around 20k voxels in the mask,
% per modality... hoping to increase this memory/cpu-related limitation.
[Y, fileinfo] = flica_load(Yfiles);
fileinfo.shortNames = {'FA', 'MD', 'AD', 'RD'};

%% Run FLICA
clear opts
opts.num_components = 5;
opts.maxits = 100;
% opts.lambda_dims = '';  # Reduce subject-noise estimation to avoid overfitting

Morig = flica(Y, opts);
[M, weights] = flica_reorder(Morig);

%% Save results
outdir = [pwd '/result/'];
mkdir(outdir)
flica_save_everything(outdir, M, fileinfo);

%% Produce correlation plots
clear des
des.Group = load('group.txt');
des.Age = load('age.txt');
des.Gender = load('gender.txt');
flica_posthoc_correlations(outdir, des);

実行が成功すると次のように画像を読み込み、計算が始まる。

Loading "all_FA_skeletonised.nii.gz"... 
Converting back to single... Generating mask... Flattening data matrix... Freeing memory... 
Loading "all_MD_skeletonised.nii.gz"... 
Converting back to single... Generating mask... Flattening data matrix... Freeing memory... 
Loading "all_AD_skeletonised.nii.gz"... 
Converting back to single... Generating mask... Flattening data matrix... Freeing memory... 
Loading "all_RD_skeletonised.nii.gz"... 
Converting back to single... Generating mask... Flattening data matrix... Freeing memory... 
Setup... doing PCA... more setup... done (0.72 sec)

Iteration 1... eta H W X.... lambda.... done updates (0.2 sec) F (0.061 sec)No previous lastF
Iteration 2... eta H W X.... lambda.... done updates (0.17 sec) F (0.029 sec), dF = 6.38536e+07
Iteration 3... eta H W X.... lambda.... done updates (0.16 sec) F (0.027 sec), dF = 21559.9
Iteration 4... eta H W X.... lambda.... done updates (0.13 sec) F (0.019 sec), dF = 253.743
...
Iteration 100... eta H W X.... lambda.... done updates (0.13 sec) F (0.019 sec), dF = 253.743

計算が終わったら、flica_report.shを実行することで、HTML形式で結果をまとめることができる。まとめられた結果(HTML)を見るには「index.html」を開く。

flica_report.sh result  # 引数としてresultフォルダを指定

結果

index.html

index.htmlの開くと次のようにFLICAの、まとめの結果を見ることができる。

f:id:AIProgrammer:20201118222137p:plain

コンポーネントごとの結果(All_in_one_page)

ICAの各Componentごとの結果は、index.htmlの「All_in_one_page」から確認することができる。

f:id:AIProgrammer:20201118222606p:plain

モダリティごとの結果(FA, MD, AD, RD, etc.)

モダリティごとの結果は、index.htmlの「モダリティ名」をクリックすることで確認できる。以下の図は、FAの結果である。

f:id:AIProgrammer:20201118224553p:plain

GBSSへの適応

GBSSの結果に対してもFLICAが使えるか、不明であったためできるか試してみた。

特に問題なく、使えそうな感じ。

index.html

f:id:AIProgrammer:20201118223204p:plain

コンポーネントごとの結果(All_in_one_page)

f:id:AIProgrammer:20201118223409p:plain

モダリティごとの結果(FA, MD, AD, RD, etc.)

f:id:AIProgrammer:20201118224503p:plain


頑張れ!喝!!の代わりにB!ブックマークを押していただけるとただただうれしいです(^^)! ↓