使用enigma-toolbox可视化分区结果

技术教程
1

一、背景

分区结果指的是基于大脑分区模板获得的不同脑区的特征或统计结果。比如，将大脑皮层划分为100个脑区，获得每个被试不同脑区的皮层厚度特征，然后对每个脑区的皮层厚度进行统计分析，最终得到每个脑区的T值或者P值。本文介绍如何使用enigma-toolbox对分区结果进行可视化。

二、下载和安装enigma-toolbox

  1. 在官方github仓库下载，我这里使用的是2.0.3版本。

  2. enigma-toolbox可以在Python或Matlab环境下使用，我这里使用的是Matlab。将上一步下载的压缩包解压后（假设解压后文件夹名为ENIGMA-v2.0.3），将ENIGMA-v2.0.3文件夹下的matlab子文件夹添加到Matlab搜索路径即可，比如:

addpath(genpath('/home/alex/software/ENIGMA-v2.0.3/matlab'))

三、使用enigma-toolbox自带分区模板

  1. 通常查看相应的函数，似乎enigma-toolbox只提供了DK、HCP-MMP（也称为Glasser360）、Schaefer（100-400脑区）和Economo-Koskinas皮层分区模板，以及Aseg皮层下分区模板。

  2. 以下的代码用于可视化DK分区的结果，DK分区总共有68个脑区，为了方便定位脑区，我这里给左半球的34个脑区赋值为-1到-34，给右半球的34个脑区赋值为1到34（在实际使用中，要可视化的信息可能是每个脑区的T值或者效应量）：

%% Plot cortical data
%% Prepare data
cort_dat = [-(1:34), 1:34];
%% Map data into fsaverage5 space
cort_dat_fsa5 = parcel_to_surface(cort_dat, 'aparc_fsa5');
%% Plot & Save
fig_dat = figure
plot_cortical(cort_dat_fsa5, 'surface_name', 'fsa5', 'color_range', [-34 34], 'cmap', 'RdBu_r')
print(fig_dat, '-dpng', '-r300', 'DK.png');
close(fig_dat);

一个不太确定的问题是，这68个脑区的顺序是什么？也就是作为输入的68个值分别表示什么脑区。因为在enigma-toolbox的文档以及代码里似乎没有明确提到。不过我猜测应该是与annotation文件里的脑区顺序对应的，在下一节中使用另一种方法来绘图时可以印证这一点。

fig01
fig011000×400 149 KB

  1. 以下代码用于可视化Aseg分区的结果，可以显示14个皮下核团结构，为了方便定位，我给左半球的7个核团赋值为-1到-7，给右半球的7个核团赋值为1到7：
%% Plot subcortical data
%% Prepare data
subcort_dat = [-(1:7), 1:7];
%% Plot & Save
fig_dat =figure;
plot_subcortical(subcort_dat, 'color_range', [-7 7], 'cmap', 'RdBu_r', 'ventricles', 'False')
print(fig_dat, '-dpng', '-r300', 'Aseg.png');
close(fig_dat);

作为输入的14个数值，前7个表示左半球，后7个表示右半球。脑区顺序为：accumbens、amygdala、caudate、hippocampus、pallidum、putamen、thalamus。

fig02
fig021000×400 60.1 KB

四、使用其他皮层分区模板

由于enigma-toolbox提供的皮层分区模板有限，如果想可视化其他分区模板，最关键的一步就是如何把分区结果（通常几十到几百个脑区）映射到皮层（通常几万到十几万个顶点）上。在FreeSurfer的annotation文件中提供了脑区与顶点之间的映射关系，因此只要有分区模板的annot文件，就可以完成这一步。为了与前面的绘图结果相互印证，我这里仍然可视化DK模板。以下是具体代码：

%% Plot DK based on annotation file
%% Load annot file
[lh_vertices, lh_label, lh_colortable] = read_annotation('/home/alex/software/ENIGMA-2.0.3/matlab/shared/annot/fsa5_lh_aparc.annot');
[rh_vertices, rh_label, rh_colortable] = read_annotation('/home/alex/software/ENIGMA-2.0.3/matlab/shared/annot/fsa5_rh_aparc.annot');
%% Map data into fsaverage5 space
%% Deal with left and right hemi separately
%% The first and fifth values represent unknown and corpuscallosum regions
lh_cort_dat = [0, -(1:3), 0, -(4:34)];
lh_cort_dat_fsa5 = zeros(1, size(lh_label,1));
for roi_idx = 1:length(lh_colortable.struct_names)
    vtx_idx = lh_label == lh_colortable.table(roi_idx,5);
    lh_cort_dat_fsa5(vtx_idx) = lh_cort_dat(roi_idx);
end
rh_cort_dat = [0, 1:3, 0, 4:34];
rh_cort_dat_fsa5 = zeros(1, size(rh_label,1));
for roi_idx = 1:length(rh_colortable.struct_names)
    vtx_idx = rh_label == rh_colortable.table(roi_idx,5);
    rh_cort_dat_fsa5(vtx_idx) = rh_cort_dat(roi_idx);
end
%% Combine left and right data
cort_dat_fsa5_v2 = [lh_cort_dat_fsa5, rh_cort_dat_fsa5];
%% Check equivalence with parcel_to_surface
all(cort_dat_fsa5 == cort_dat_fsa5_v2)
%% Plot & Save
fig_dat = figure;
plot_cortical(cort_dat_fsa5_v2, 'surface_name', 'fsa5', 'color_range', [-34 34], 'cmap', 'RdBu_r');
print(fig_dat, '-dpng', '-r300', 'DK_v2.png');
close(fig_dat);

这里使用enigma-toolbox自带的annot文件，读取annot文件使用enigma-toolbox自带的read_annotation函数（打开函数会发现该函数是直接从FreeSurfer复制的），annot数据（以左半球为例）中lh_colortable.struct_names表示脑区名字和顺序，lh_colortable.table的第5列表示每个脑区对应的annotation值（也就是每个脑区的数字标签，只是不是1, 2, 3, ...这种编码方式），lh_label表示每个顶点的annotation值。关于annotation文件格式请参考FreeSurfer的文档。通过比较根据enigma-toolbox的parcel_to_surface函数得到的cort_dat_fsa5和根据annotation文件得到的cort_dat_fsa5_v2，会发现完全一致。同样地，使用plot_cortical绘图会得到相同的结果。

（本文部分内容最初于2019和2023年发布在个人公众号上，这里重新进行了整理，并补充了一些内容。）

1 个赞