'뇌영상 분석'에 해당되는 글 6건

  1. Contrasts in Neuroimaging Data Anlaysis (1)
  2. Slice Timing Correction
  3. smoothness estimation in SPM and AFNI (2)
  4. Voxel size determination through the MonteCarlo Simulation in AFNI
  5. AFNI에서 ROI 만들기.
  6. Data Analysis (1): Neuroimaging Data loading using SPM8 toolbox (11)

SPM 등의 뇌영상 데이터 분석 툴을 이용한 뇌영상 데이터의 분석은 기본적으로 각 복셀의 영상에 할당된 데이터 값을 일반 선형 모델 (General Linear Model, GLM)을 이용하여 모델링하고, 실제 데이터와 모델이 얼마나 잘 맞는지 통계적으로 테스트 하는 것이다. 특정 복셀 $i$에 대해서 $Y_{i} = XB_i + E_i$로 모델링 했을때 $X$는 디자인 행렬이되고, 벡터 $B_i$는 분석을 통해서 추정되는 파라미터이며, $E_i$는 에러를 의미한다. 이때 contrast는 $c'B$를 통해서 계산된다. 뇌영상 데이터에서 $c$는 보통 행벡터(column vector)를 의미하고, $c$를 통해서 다양한 contrasts로 결과를 확인할 수 있다. 

벡터 $c$는 contrasts의 가중치를 의미하고, 확인하고자 하는 대부분의 contrast의 Null hypothesis는 $c'B=0$이다. 구체적으로 contrast 가중치를 구성하는 방법에 대해서 다음의 예를 통해서 설명하고자 한다.

위 그림과 같이 세 그룹에 대한 데이터가 있을때, 가장 먼저 생각해 볼 수 있는 통계 분석은 일원분산분석(Analysis of variance, ANOVA)이고 이때의 contrast는 F-contrast로 세 그룹 중 발생할 수 모든 차이를 보고자 할때는 다음과 같은 방법으로 만들 수 있다.

또는 결과는 동일하지만 다른 방법으로 표현하고 싶다면, 다음과 같은 방법도 가능하다.

이제 각 그룹간 차이를 보고자 한다면, 다음과 같은 방법으로 분석이 가능하다. 하지만 분석의 초점이 ANOVA에서 그룹간 차이가 나는 영역을 찾고 해당 영역 내에서 세 그룹간 값이 어떻게 되는지 통계적으로 분석하는 것이 목표라면, F-contrast의 결과로 나온 영역에 대해서 평균 Beta 값을 추출해서 SPSS 등에서 post-hoc 분석을 진행하는 것이 옳다. MATLAB에서 post-hoc 분석이 가능하다면 말리지는 않겠다. SPM에서 T-contrast는 다음과 같은 방법으로 만들 수 있다.

  • -1    1    0:  spm{T} contrast for sample2 > sample1

  •  1     0   -1:  spm{T} contrast for sample1 > sample3

  •  2   -1   -1:  spm{T} contrast for sample1 > (sample2 + sample3)

이제 confounding effect를 제어한 후에 그룹간 차이를 비교하는 분석을 생각해 보자. 이러한 분석을 공변량분석 (Analysis of Covariance, ANCOVA) 이라 불린다. 다음과 같은 데이터가 있다고 가정해 보자.


이제 공변량(covariate1)의 효과를 제거한 후에 두 그룹간의 차이를 비교하는 F-contrast를 생각해 보다. 다음과 같은 spm{F} contrast를 통해서 sample1과 sample2 간의 값의 차이를 비교할 수 있다.



신고

Slice Timing Correction

하나의 3차원 뇌영상 데이터는 여러개의 단면영상(Slice Image)으로 구성되어 있다. 기능자기공명영상(functional magnetic resonance imaging, fMRI)의 경우에는 보통 매2초마다 하나의 3차원 볼륨 영상을 획득하게 된다. 다시 말해서 2초동안 여러개의 단면영상을 획득하게 되는데, 그중에서 제일 처음에 획득한 단면영상과 맨 마지막에 획득한 단면영상 간에는 최고 2초 정도의 시간 차이가 발생하게 된다. 이러한 시간 차이를 보정해 주는 것이 Slice timing correction이라 물리는 전처리 과정이다. SPM의 경우에 slice timing correction 을 적용하게 되면 'a'를 어두로 하는 새로운 뇌영상 파일이 생성된다.

단면영상의 순서(slice order)는 각 단면영상이 획득된 순서를 의미한다. 단면영상에서 순서로 정의해준 순서는 각 단면영상이 획득된 시간 순서를 의미한다. SPM에서 order of slices를 확인하기 위해서는 SPM Display에서 보고자 하는 뇌영상을 선택하고, cross-hairs를 z=1인 위치로 이동시키면 볼륨 영상의 첫번째 단면 영상을 가르키게 된다.

Slice timing 을 통해서 동일한 볼륨 영상 내에 서로 다른 시간에 획득된 단면 영상들이 마치 동일한 시간에 획득된 것 처럼 보정을 하게 된다. 보정 과정은 sinc-interpolation을 이용하고, 이것은 시계열 데이터를 lagging하는 효과를 낸다.

개념을 명확히 하기 위해서 예를 들어 설명해 보고자 한다. 가령 연속되는 단면 영상에 어떤 신경활동이 있다고 가정해 보자. A 단면영상의 값들은 t=0s에 획득되었고, B 단면영상의 값들은 t=1s에 획들되었다. 보정이 없다면 B 단면영상에 있는 혈류동역학 반응(hemodynamic response)은 사실상 B 영상이 획득되기 1초 전에 발생한 신호를의미한다. 이것을 보정하기 위해서는 B 단면영상에 있는 값들을 1초 전에 발생한 시점의 값들로 shifting 해야 한다.

Slice timing correction에서 중요한 점은 뇌영상 데이터를 획득한 순서를 명확히 알아야 한다는 것이다. 만약 잘 모른다면 뇌영상 획득에 도움을 준 방사선사 선생님들께 문의해서 영상 획득 순서를 숙지하고 있는 것이 좋다. 보통은 ascending (from bottom to top), descending (from top to bottom), interleaved (odd number slices and then even number slices)등의 순서로 획득된다.

지멘스Siemens MRI로부터 .ima 파일로 데이터를 저장했다면, MATLAB에서 다음과 같은 명령어로도 order of slice를 확인하는 것이 가능하다.

>>hdr = spm_dicom_headers(’dicom.ima’); >>slice_times = hdr.Private_0019_1029


이 경우 단면영상의 순서는 아래쪽으로부터 위쪽으로 증가하는 방향이다.


신고

Smoothness estimation은 MonteCarlo simulation을 위해서 반드시 필요한 과정입니다.

SPM으로 영상 데이터를 분석했다면, SPM.mat 파일의 Field 값을 확인함으로써 smoothness를 확인할 수 있습니다.

>>load SPM; % SPM 결과 파일이 저장된 폴더에서 실행 >>M = SPM.xVol.M; % 변환행렬 정보를 가져옴 >>VOX = abs(diag(M)); % 대각행렬 정보가 볼셀 사이즈 >>FWHM = SPM.xVol.FWHM; % FWHM in voxel unit >>FWHMmm = FWHM.*VOX(1:3)'; % FWHM in mm unit >>disp(FWHMmm); 

SPM에서 Gaussian random field theory를 기반으로 smoothness를 estimation하고, 이것은 spm_est_smoothness을 통해서 계산됩니다. 

위의 과정을 통해서 확인한 FWHMmm 정보는 AFNI의 3dFWHMx의 명령어를 통해서 estimation한 smoothness와 완전히 동일하지는 않지만, 비슷한 수준의 값이 나와야 합니다. AFNI가 설치되어 있다면 다음과 같이 실행해 보세요. 

~>3dFWHMx -mask mask.hdr -input ResMS.hdr % SPM 결과 파일이 저장된 폴더에서 실행

여기서 얻은 smoothness를 기반으로 MonteCarlo simulation을 시행하면, 이것이 AlphaSim-corrected p-value가 됩니다.

관련 글은 다음의 링크에서도 확인이 가능합니다.
https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=spm;bb3ad6d.1308

신고

뇌영상을 이용한 뇌과학 연구는 대부분 다음과 같은 질문에 답하는 것이다.

  • 3차원의 뇌에 어떤 영역이 활성화 되었는지?  
  • 특정 영역의 시계열 데이터와 상관성이 높게 나오는 뇌영역은 어디인지? 

생물정보학에서 다루는 Microarray 데이터도 마찬가지 이지만, 뇌영상 데이터도 multiple comparison에서 발생하는 false positive를 조절하는 방법에 대해서다양한 해법들이 있다. 가령, false discovery rate (FDR) 또는 family-wire error rate (FWE) 등이 전통적으로 가장 많이 이용되어 왔던 multiple comparison correction 방법들이다. 

뇌영상 데이터는 특정 복셀에서 통계적으로 유의미한 차이를 보인다고 했을때 "아, 이 볼셀에서 통계적으로 유의미 하니까 해당 영역이 중요하겠네!" 라고 생각하면 안된다. 뇌는 인접한 뇌영역들 같이 비슷한 역할을 하는 뉴런들이 있어서 특정 복셀 1개만 통계적으로 유의미한 차이를 보이고, 인접한 영역에 있는 복셀들에서도 비슷한 양상이 보여야 한다. 전문 용어로 차이가 나는 영역 덩어리 (continuous voxel size)가 얼만큼 큰지?의 여부에 따라서 해당 영역이 false positive인지를 판별해야 한다. 최근 뇌영상 연구에서는 multiple comparison correction을 하지 않은 분석 결과는 좋은 저널에 출판하기 어려울 것이다. 

뇌영상 연구의 Multiple comparison에서 중요한 것은 다음의 두가지 질문일 것이다.

  • 각 복셀에서의 통계적 유의미성이 어느정도로 해야 하는지? (보통 뇌영상(fMRI or MRI) 분석에서는 p<0.005 or p<0.001을 많이 사용) 
  • 해당 threshold에서 덩어리져 있는 영역-덩어리를 찾아야 할텐데, 덩어리(or continuous voxel size)의 크기는 얼만큼이 적당한지?  

위의 질문은 MonteCarlo (MC) simulation을 통해서 할 수 있는데, SPM with REST toolbox와 AFNI의 AlphaSim을 통해서 시뮬레이션이 가능하다. MC simulation을 위해서는 여러가지 파라미터가 필요한데, 각각의 파라미터의 이름은 프로그램마다 조금씩 다르지만, 대략적으로 다음과 같은 의미를 가진다고 할 수 있다.

  • mask: MC simulation을 하는 voxel size, image dimension, and brain boundary 등의 정보를 담고 있음.
  • rmm: 얼만큼 떨어져 있는 복셀을 인접한 복셀로 볼 것인지? 볼셀의 중심과 중심의 거리
  • NN: nearest neighbor를 어떻게 정의할 것인지?  가령, 각각의 볼셀이 6면체로 되어 있는데 복셀의 face가 붙어 있는 것을 연결된 복셀로 볼 것인지? edge가 붙어 있는 것을 연속적으로 연결된 복셀로 볼것인지? 아니면 corner가 붙어 있는 것을 연속적으로 연결된 볼셀로 볼 것인지? 의 여부를 결정할 수 있다.
  • FWHM: 이것은 결과 데이터의 smoothness를 의미함. 단순히 smoothing kernel의 FWHM 값을 입력하면 안되고, AFNI나 SPM의 smoothness estimation 방법을 통해서 계산한 결과 값을 입력해야 한다. smoothness estimation을 다룬 포스팅을 확인하기.



신고

AFNI에서 ROI 만들기.

기능 뇌영상 연구에서는 특정한 관심 영역(Region of Interest, ROI)를 설정하고, 해당 ROI에 대한 시계열 데이터를 추출해서 분석을 해야 하는 경우가 많이 있습니다. AFNI를 통해서 ROI를 그리는 방법에 대해서 잘 설명되어 있는 영문 포스트가 있어서 이곳에 한글로 옮깁니다.

Afni 프로그램은 리눅스나 매킨토시 계열의 컴퓨터에서 작동하는 프로그램으로, 프로그램을 설치하고 쉘Shell에서 PATH 설정을 제대로 했다면, @GetAfniBin 명령어를 통해서 Afni가 설치되어 있는 경로를 return 받을 수 있습니다. 

이제 중심 좌표의 위치가 (x=23,y=21,z=-6.5 in talairach coordinate)이고 반지름이 3mm인 왼쪽 hippocampus 영역의 일부 영역을 구Shpere 형태의 ROI로 만들고 싶다면, 다음의 명령어를 통해서 가능합니다.

3dcalc -a `@GetAfniBin`/TT_N27+tlrc.HEAD \

-prefix ShpereROI \

-expr 'step(3*3 - (x-23)*(x-23) - (y-21)*(y-21) - (z+6.5)*(z+6.5))'

아주 간단하죠?

TT_N27+tlrc는 Afni나 Caret 등에서 볼 수 있는 Talairach 표준 템플릿 입니다.


신고

일반적인 데이터 분석에 관련된 기술을 포스팅하고 싶었지만, 데이터 분석은 data specific한 부분들이 있기 때문에 일반적인 데이터 분석 및 시각화 기법에 대한 강의는 어려울것 같고, 뇌영상 데이터와 트위터 데이터 수집 및 분석과 관련된 포스팅을 연재하고자 합니다.

일단 뇌영상 데이터 분석 기법으로 강의를 시작하는 이유는 최근 5년 동안 제가 연구해온 분야로 expert까지는 아니어도 intermediate 이상의 실력은 된다고 생각하고 있기에, 지금까지 습득한 노하우를 관련 분야에 계신 분들께 나누고 싶은 생각이 들었습니다.

뇌영상 분석에 사용되는 툴은 여러가지가 있지만, 보통은 MATLAB 기반의 SPM8, linux 기반에서 작동하는 AFNI 또는 FSL 등이 가장 많이 이용되고 있습니다.어떤 소프트웨어가 가장 좋은가? 라는 질문의 대한 답은... 내가 사용하기 편한 소프트웨어. 가 아닐까? 생각됩니다. 우수한 프로그램은 고성능의 알고리즘으로 무장하여 계산속도를 극대화 하는 것도 중요하지만, 사용자가 얼마나 편하게 접근할 수 있는가?에도 많은 노력을 기울여야 한다고 생각합니다.

필자는 MATLAB이 사용하기 편하고, engineering 기반이 아닌 사람들과 소통하기에 (뇌연구는 심리학, 정신의학, 물리학, 통계학, 전자공학, 의공학 등 다양한 전공을 가진 연구자들이 있기 때문에 공동연구와 연구자들간의 소통이 중요함) 가장 적합한 프로그램이라 생각하고, MATLAB을 사용하기 시작했습니다. 앞으로 'Data Analysis (for neuroscientist)'라는 주제로 한글 인터넷 강좌를 연재하고, MATLAB 뿐만이 아니라 python 및 linux script에 대한 내용도 다룰 예정입니다. 뇌영상 데이터를 연구 재료로 사용하고 계신 연구자들 중에서, GUI 환경의 분석에서 벗어나고자 하시는 분들께 조금이라도 도움이 되었으면 좋겠네요^^ 오늘은 그 첫번째 강좌로 "Data loading and writing using SPM8 toolbox"에 대한 강의를 하고자 합니다.

SPM8은 프로그래밍을 잘 몰라도 사용할 수 있는 막강한 GUI환경을 제공합니다. 이 GUI를 이용하면 뇌영상 데이터(fMRI, EEG, PET 등)의 전처리를 수행할 수 있습니다. 물론 GUI를 사용하지 않고, MATLAB command line을 통해서 더욱 편리하게 전처리를 할 수도 있지만, 그러기 위해서는 프로그래밍에 대한 기초 지식이 있어야 해요. 하지만 겁먹을 필요는 없고요, 앞으로 하나하나 배워나가면 됩니다.

일단 MATLAB에서 SPM8의 툴발스를 이용하여 뇌영상 데이터를 불러들이고, 또 어떤 연산 결과를 새로운 이름의 영상으로 저장하고자 한다면,

(1) SPM8 toolbox를 다운로드 받아야 합니다.

(2) 예제로 사용할 뇌영상 데이터가 필요합니다. 데이터가 없다면 SPM 홈페이지에 있는 auditory fMRI data를 사용하실 수 있습니다. (MoAEpilot.zip)


이제 다운로드 받은 영상 데이터를 matrix 형태로 불러들이기 위해서는 MATLAB을 실행시키고, 다음과 같이 명령어를 순서대로 입력하면 됩니다.

>addpath /Users/skyeong/spm8; % spm8에서 제공하는 모든 함수를 어떠한 경로에서도 사용할 수 있도록 설정. >data_path = '/Users/skyeong/Downloads/MoAEpilot/fM00223' ; % fMRI 데이터 위치의 절대경로 >fns = spm_select('FPList',fullfile(data_path),'^fM.*\.img$'); % 파일 이름이 fM으로 시작하고, 확장자가 .img인 모든 파일을 선택. >vs = spm_vol(fns); % 뇌영상의 header 정보 불러오기. >IMG = spm_read_vols(vs); % 뇌영상 데이터를 matrix 형태로 불러오기.

vs 변수에는 뇌영상 데이터의 header 정보가 저장되어 있습니다. 다음과 같이 MATLAB command line에서  vs(1) 이라고 입력하면, 첫번째 뇌영상 볼륨의 정보를 볼 수 있어요. 

>vs(1) ans = fname: '/Users/skyeong/Downloads/MoAEpilot/fM00223/fM00223_004.img' dim: [64 64 64] mat: [4x4 double] pinfo: [3x1 double] dt: [4 1] n: [1 1] descrip: 'SPM compatible' private: [1x1 nifti]

여기서 fname은 첫번째 볼륨vs(1)이 위치한 full file path를 의미하고, dim은 뇌영상 데이터의 dimension정보, dt(=data type)는 뇌영상 데이터의 저장 형태입니다. 보통 프로그래밍 언어에서는 boolean, integer, float 에 따라서 데이터를 저장하는 용량이 다른데, [4 1]은 복셀 하나에 저장되는 숫자가 4 bytes로 표현된다는 의미이고, 보통 4 bytes 이면 양의 정수를 표현하기에 적합하다고 생각하면 됩니다. 만약 어떤 연산의 결과가 소수점을 동반한다면 vs(1).dt=[16 1]로 설정해야 하고, 이렇게 하면 한 복셀당 16 bytes의 용량이 필요하기 때문에, 4 bytes로 저장될때보다 뇌영상 데이터의 용량이 4배로 늘어나게 되요. data type이  높다고 좋은것은 아니고, 계산결과가 소수형인지, 양의 정수형인지, 음의 정수를 포함하는지를 잘 판단하여 설정하는 것이 하드디스크의 용량을 적게 사용할 수 있는 노하우입니다. MATLAB의 Numeric Types에 대해서는 MathWorks의 홈페이지를 참고하면 자세히 알아볼 수 있어요.

이제 IMG 변수에 담긴 숫자들은 뭘 의미하는지 알아 볼까요? IMG변수에는 뇌의 활동에 따른 BOLD signal이 숫자의 형태로 각 볼셀에 저장되어 있습니다.

>size(IMG) ans = 64 64 64 96

MATLAB command line에 size(IMG)를 입력하면 IMG 변수가 4차원 형태로 되어 있음을 알 수 있다. 처음 3개 숫자는 x, y, z 축의 복셀 숫자를 의미하고, 네번째 숫자인 96은 시간축 정보를 의미합니다. 즉, 실험을 통해서 96개의 volumes을 획득한 것이고, 만약 repetition time (TR)이 2초 였다면, 96x2 sec = 192 sec 동안 fMRI 영상을 획득한 셈이겠죠? 

이제 첫번째 볼륨vs(1)의 intensity 값의 분표를 MATLAB의 hist() 함수를 이용하여 확인해 볼께요.

>figure; >data = IMG(:,:,:,1); % 첫번째 volume의 intensity 정보를 data라는 변수에 저장. >hist(data(:)); % 모든 복셀의 intensity 값의 분포를 확인. >figure; >hist(data(:),100); % 모든 볼셀의 intensity 값의 100개의 bins으로 나눠서 분포를 확인.


강좌를 올리는게 쉽지는 않네요, 몇줄 안되는 코드이지만 프로그래밍은 직접 실습해보는게 중요합니다. 관심 있으신 분들은, 꼭 집접 확인해 보시길. 다음 강좌에서는 MATLAB을 이용한 regression analysis, bandpass filtering, network analysis 등을 소개할 예정입니다. 질문은 다른 분들도 함께 공유할 수 있도록 블로그나 패이스북의 덧글로 부탁드립니다.

신고