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

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

【MATLAB】ComBatを用いた交絡因子の排除 (実践編)

f:id:AIProgrammer:20200507234419p:plain

問題

グループ間の差を解析する際、計測した値は、年齢、性別、計測器等の交絡因子の影響を含む。

目的

年齢、性別、計測器等の交絡因子の影響を排除するように計測値を補正する。 詳細は、こちらのManualをみること。

準備

ComBatのDownload

データの補正には、ComBatというSoftwareを用いる。 ComBatのサイトから、必要なファイルをダウンロードする。 (Qiitaの画像アップロード月制限で画像載せれません。)

パスを通す

MATLABを起動し、Downloadした、"./Matlab/scripts/"にパスを通す。 パスの追加方法は、こちら

または、[パスの設定] ダイアログ ボックスを使用して MATLAB 検索パス全体を表示します。[ホーム] タブの [環境] セクションで、[パスの設定] をクリックします。[パスの設定] ダイアログ ボックスは、検索パス上のすべてのフォルダーを一覧表示します。[パスの設定] ダイアログ ボックスの使用の詳細については、検索パス上のフォルダーの変更を参照してください。

Datasetの準備

以下のように、列がSubject、行が交絡因子や計測値やROIの種類となるようにデータ(CSVファイル)を用意。

  • 1行目は、Cite(計測器, ex.) MRI scaner)。すべての被験者が同じMRIで撮像していればすべての被験者に1をふる。番号は1からふらないとエラーになります。0からふらない。
  • 2行目は、年齢
  • 3行目は、性別を記入。1が男性、2が女性となるように番号をふる。番号は1からふらないとエラーになります。0からふらない。
  • 4行目は、Group。必ず健常者が1になるように番号をふる。また、番号は1からふらないとエラーになります。0からふらない。
Sub001 Sub002 Sub003 ・・・ SubN
cite(1:A univ., 2:B univ, 3:C univ) 1 3 2 ・・・ XXX
Age(y) 56 56 56 ・・・ XXX
Sex(M:0,F:1) 1 1 1 ・・・ XXX
Group 1 1 1 ・・・ XXX
Value_A 2040 2755 2365 ・・・ XXX
Value_B 882 2259 1397 ・・・ XXX
Value_C 4506 6122 6384 ・・・ XXX
・・・ ・・・ ・・・ ・・・ ・・・ XXX
Value_? XXX XXX XXX XXX XXX

Datasetの欠損値(NAN)の確認

Datasetの中には、欠損値(NAN)が含まれているときがある。 ComBatでは、欠損値が含まれているとエラーが起きるため、NANは、0と置き換えて処理している。 ComBatを走らせる時点で、NANが含まれたDataset(CSVファイル)をImportすることができるが、まれに"#NAME?"というセルが含まれていることがあり、これをNANに置き換えてからでないとエラーになる。

ComBatの実行(交絡因子の排除)

以下のコードを作成し、"run_ComBat.m"として保存。 このとき、先ほど作成し保存したDataset(CSVファイル)と同じフォルダに保存する。

clear
clc

%import data
[file,path] = uigetfile('*.csv','ComBatにinputするファイルを選択してください');
all_data = importdata([path file]);
% 欠損地(NAN)を0で埋める。
all_data.data(isnan(all_data.data)) = 0;
% Valueの最初の行数
p = 5;
s = input('Parametric(1) or Non-Parametric(0):');
value = all_data.data(p:end,:);
cite = all_data.data(1,:); % 最小値は1。0から番号を振らない。
age = all_data.data(2,:)';
sex = all_data.data(3,:)';
group = dummyvar(all_data.data(4,:)'); % 最小値は1。0から番号を振らない。
mod = [age sex group];
data_harmonized = combat(value,cite, mod, s);
csvwrite('post_ComBat_data.csv',data_harmonized);

作成した"run_ComBat.m"を実行。スクリプトをドラッグ&ドロップでも可能。

>> run('run_ComBat.m')

コマンドラインで以下のような、質問をされる。 データが正規分布(Parametric)にしたがっていれば1。 非正規分布(Non-Parametric)であれば0。

Parametric(1) or Non-Parametric(0):

以上の質問に答えて数字を打ち込むと、以下のようなログが流れる。

[combat] Found 1 batches
[combat] Adjusting for 2 covariate(s) of covariate level(s)
[combat] Standardizing Data across features
[combat] Fitting L/S model and finding priors
[combat] Finding parametric adjustments
[combat] Adjusting the Data

ComBatの結果を確認。(交絡因子が排除されたデータ)

ComBatを実行すると、交絡因子が排除された(補正された)データがpost_ComBat_data.csvとしてOutputされる。 中身をみるとこんな感じ。 ComBatの出力される結果は、計測値(ROIの値)のみです。 Inputしたdatasetの6行目以降(Groupの行より下)の部分のみが出力されます。

0.10068 0.084548 0.064942 0.079538 0.076502
0.06857 0.093662 0.097684 0.10724 0.079938
0.14628 0.11298 0.084449 0.1087 0.084588
0.066614 0.037983 0.0047462 0.028432 0.031116
0.09353 0.10883 0.056757 0.07756 0.080544

例えば、このようなdatasetをInputした出力されるのは、Value_A, Value_B, Value_C,...Value_?までの範囲です。 元ある、値とComBatの出力結果を入れ替えて終了です。

Sub001 Sub002 Sub003 ・・・ SubN
cite(1:A univ., 2:B univ, 3:C univ) 1 3 2 ・・・ XXX
Age(y) 56 56 56 ・・・ XXX
Sex(M:0,F:1) 1 1 1 ・・・ XXX
Group 1 1 1 ・・・ XXX
Value_A 2040 2755 2365 ・・・ XXX
Value_B 882 2259 1397 ・・・ XXX
Value_C 4506 6122 6384 ・・・ XXX
・・・ ・・・ ・・・ ・・・ ・・・ XXX
Value_? XXX XXX XXX XXX XXX



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