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

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

【脳MRI画像解析】ComBatをもちいたデータの調和(Harmonization)

f:id:AIProgrammer:20210318172303p:plain

動かしながら学ぶ PyTorchプログラミング入門

動かしながら学ぶ PyTorchプログラミング入門

  • 作者:斎藤勇哉
  • 発売日: 2020/11/30
  • メディア: Kindle版

目的

  • 多施設MRIデータに含まれる施設間差(MRI装置・撮像条件)の評価
  • ComBatを用いて施設間差を除去(データの調和)
  • 施設間差補正後の評価
  • 拡散定量値を用いて、施設間差を評価

背景

多施設で取得されたMRIデータは、施設ごとに異なる撮像装置および撮像条件で取得されたものであり、施設間差(site effect)を含む。そのため、この施設間差と本来観測したい個人間の差(ex. 疾患・成長・老化による影響)を分離することができず、単に多施設データをするだけでは、感度・特異度に劣る。特に、精神疾患のような脳神経変性がわずかな場合、施設間差による変動に埋もれて、かすかな病変を検出できない。そのため、多施設共同研究において施設間差を除去することは、重要な課題である。

参考資料: mriportal.umin.jp

準備

施設間差を除去するために、ComBatを用いる。ComBatは、MATLAB・Python・Rで動作するが、ここではPython版の使い方を解説する。

Python版のComBat (neuroCombat)をインストールするには、以下のコマンドを実行する。

pip3 install neuroCombat

neuroCombatの使い方は、次の通り。

neuroCombat(dat,
           covars,
           batch_col,
           categorical_cols=None,
           continuous_cols=None,
           eb=True,
           parametric=True,
           mean_only=False,
           ref_batch=None):
    """
    Run ComBat to remove scanner effects in multi-site imaging data
    Arguments
    ---------
    dat : a pandas data frame or numpy array
        - neuroimaging data to correct with shape = (features, samples) e.g. cortical thickness measurements, image voxels, etc
    covars : a pandas data frame w/ shape = (samples, covariates)
        - contains the batch/scanner covariate as well as additional covariates (optional) that should be preserved during harmonization.
        
    batch_col : string
        - indicates batch (scanner) column name in covars (e.g. "scanner")
    categorical_cols : list of strings
        - specifies column names in covars data frame of categorical variables to be preserved during harmonization (e.g. ["sex", "disease"])
    continuous_cols : list of strings
        - indicates column names in covars data frame of continuous variables to be preserved during harmonization (e.g. ["age"])
    eb : should Empirical Bayes be performed?
        - True by default
    parametric : should parametric adjustements be performed?
        - True by default
    mean_only : should only be the mean adjusted (no scaling)?
        - False by default
    ref_batch : batch (site or scanner) to be used as reference for batch adjustment.
        - False by default
        
    Returns
    -------
    - A numpy array with the same shape as `dat` which has now been ComBat-harmonized
    """

使用するデータ

各被験者ごとに計測した定量値を、次のようにCSVに保存する (data.csv)。

以下の例では、拡散MRIから得られた拡散定量値を、白質路(JHU DTI-based white-matter atlases)ごとにサンプルした結果をまとめている。

各変数の詳細は、次の通り。

  • Region:白質路(48領域)
  • Map:拡散定量値(FA, MD, ICVF, ISOVF, OD)
  • Site:MRI撮像施設(SiteA, SiteB, SiteC)
  • Scanner:MRI装置(ScannerA)
  • Protocol:撮像プロトコル(ProtocolA)
  • Subj???:被験者(Subj001, Subj002, Subj003)

    data.csv

    Region Map Site Scanner Protocol Subj001 Subj002 Subj003
    MCP FA SiteA ScannerA ProtocolA 0.598509 0.607412 0.590479
    PCT FA SiteA ScannerA ProtocolA 0.483564 0.524959 0.447631
    GCC FA SiteA ScannerA ProtocolA 0.722272 0.771941 0.75223
    UF_L ISOVF SiteC ScannerA ProtocolA 0.117252 0.126971 0.147088

ソースコード

データ(data.csv)の準備ができたら、以下のスクリプトを実行する。

やっている内容としては、

  • 補正前のデータに対して、相関係数(spearmanの順位相関係数)の混合行列を生成
  • ComBatでデータを補正し、相関係数(spearmanの順位相関係数)の混合行列を生成
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import os
from neuroCombat import neuroCombat


# Define
filename = 'data.csv'
cmap = 'Spectral_r'
hm_vmax = 1
hm_vmin = 0.8
hm_center = 0.9


def arrange_data(df, map):
    columns = ["{}_{}_{}_{}".format(id, site, scanner, protocol)
               for id in df.columns[5:]
               for site in df['Site'].unique()
               for scanner in df['Scanner'].unique()
               for protocol in df['Protocol'].unique()]
    df_arrange = pd.DataFrame(index=df['Region'].unique(), columns=columns)
    df_arrange.index.name = 'Region'
    for id in df.columns[5:]:
        for i in range(df.shape[0]):
            col = "{}_{}_{}_{}".format(id,
                                       df.loc[i, 'Site'],
                                       df.loc[i, 'Scanner'],
                                       df.loc[i, 'Protocol'])
            df_arrange.loc[df.loc[i, 'Region'], col] = df.loc[i, id]
    df_arrange = df_arrange.astype(float)
    return df_arrange


def corr_heatmap(df, output_base, map, hm_vmax, hm_vmin, hm_center):
    # Arrange data array
    df_arrange = arrange_data(df, map)
    df_arrange.to_csv(output_base+'.csv')
    # Correlation analysis
    df_corr = df_arrange.corr(method='spearman')
    # print(map, df_corr.min().sort_values())
    # Heatmap
    sns.heatmap(df_corr,
                vmax=hm_vmax, vmin=hm_vmin, center=hm_center,
                cmap=cmap, annot=True, square=True, fmt='.3f', annot_kws={'size': 8})
    plt.title(output_base.split('/')[1])
    plt.tick_params(labelsize=8)
    plt.savefig("{}.png".format(output_base),  dpi=300, bbox_inches="tight")
    plt.close()


def ComBat(df, map, categorical_cols=None, continuous_cols=None, parametric=True):
    dat = df.drop(['Region', 'Map', 'Site', 'Scanner', 'Protocol'], axis=1).T
    covars = pd.DataFrame(df['Site'].replace({'SiteA': 1, 'SiteB': 2, 'SiteC': 3}))
    batch_col = 'Site'
    # print('dat:', dat)
    # print('covars:', covars)
    # print('batch_col:', batch_col)
    dat_combat = neuroCombat(dat,
                             covars,
                             batch_col,
                             categorical_cols=categorical_cols,
                             continuous_cols=continuous_cols,
                             eb=True,
                             parametric=parametric,
                             mean_only=False,
                             ref_batch=None)
    df_combat = df.copy()
    df_combat.iloc[:, 5:] = dat_combat.T
    # print('df_combat:', df_combat)
    return df_combat


def main():
    # Make directory
    if not os.path.isdir('result'):
        os.makedirs('result')
    # Read csv
    df = pd.read_csv(filename)
    for map in df['Map'].unique():
        df_map = df[df['Map'] == map].reset_index(drop=True)
        # Make correlation heatmap
        output_base = "result/{}_before".format(map)
        corr_heatmap(df_map, output_base, map, hm_vmax, hm_vmin, hm_center)
        # ComBat
        df_combat = ComBat(df_map, map,
                           categorical_cols=None,
                           continuous_cols=None,
                           parametric=False)
        # Make correlation heatmap
        output_base = "result/{}_after".format(map)
        corr_heatmap(df_combat, output_base, map, hm_vmax, hm_vmin, hm_center)


if __name__ == '__main__':
    main()

Result

処理が完了すると、Resultフォルダが生成され、そこに結果が保存される。

「before」は補正前、「after」は補正後の結果である。

result/
├── FA_after.csv
├── FA_after.png
├── FA_before.csv
├── FA_before.png
├── ICVF_after.csv
├── ICVF_after.png
├── ICVF_before.csv
├── ICVF_before.png
├── ISOVF_after.csv
├── ISOVF_after.png
├── ISOVF_before.csv
├── ISOVF_before.png
├── MD_after.csv
├── MD_after.png
├── MD_before.csv
├── MD_before.png
├── OD_after.csv
├── OD_after.png
├── OD_before.csv
└── OD_before.png

例えば、FAの場合、補正前後での相関係数混合行列は、次のようになる。

f:id:AIProgrammer:20210318172829p:plain


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

動かしながら学ぶ PyTorchプログラミング入門

動かしながら学ぶ PyTorchプログラミング入門

  • 作者:斎藤勇哉
  • 発売日: 2020/11/30
  • メディア: Kindle版