Beauchamp:ROIanalysis

If you would like to compare activity across subjects in threshold-defined areas, an ROI analysis may be right for you!

For example, we may be interested in comparing amplitude of activity in auditory cortex during different types of speech stimuli in people with different behavioral characteristics. While we could use a group whole-brain analysis to see if there are any voxels in a combined map within auditory cortex that show a differential response, it may be better to separately define each auditory ROI in each subject using strict and reproducible criteria, such as T stat > 3.

To facilitate this process, the surfaces in SUMA come already parcellated into useful regions. With the parcellation, each anatomical region is assigned a value (same in left and right hemisphere). So, for each subject, we can find active areas within these anatomical locations to use as ROIs. Depending on when the surface was made, the parcellation is done either with a 2005s or 2009s atlas.

Here are the parcellation values corresponding to 5 regions of interest using the 2005s atlas:

 Auditory cortex                        33, 81
 Extrastriate visual cortex (V5/MT)     60
 Fusiform gyrus                         17, 63
 Inferior frontal gyrus                 6,  72
 Superior temporal sulcus               80

And for the 2009s atlas:

 Auditory cortex                        466, 514
 Extrastriate visual cortex (V5/MT)     493
 Fusiform gyrus                         450
 Inferior frontal gyrus                 285, 505
 Superior temporal sulcus               513

To use these parcellations, first copy over the parcellation files from the SUMA folder:

 cd /Volumes/data9/surfaces/last_name/subjID/SUMA/				
 cp lh.aparc.a2005s.annot.1D.roi /Volumes/data1/UT/{$ec}/afni				
 cp rh.aparc.a2005s.annot.1D.roi /Volumes/data1/UT/{$ec}/afni				
 cd /Volumes/data1/UT/{$ec}/afni				

Then, convert them into BRIK files for AfNI analysis:

 3dSurf2Vol -spec /Volumes/data9/surfaces/netek_anne/HO/SUMA/both.spec -surf_A lh.smoothwm.asc \				
 -grid_parent {$ec}v1mr+orig -sv {$ec}_SurfVol_Alnd_Exp+orig -map_func max -prefix L_parc -sdata_1D lh.aparc.a2005s.annot.1D.roi				

 3dSurf2Vol -spec /Volumes/data9/surfaces/netek_anne/HO/SUMA/both.spec -surf_A rh.smoothwm.asc \
 -grid_parent {$ec}v1mr+orig -sv {$ec}_SurfVol_Alnd_Exp+orig -map_func max -prefix R_parc -sdata_1D rh.aparc.a2005s.annot.1D.roi

The Freesurfer parcellation gives us the entire length of the STS, but what if you are interested only in the posterior portion of this sulcus? It seems fair to consider only the posterior half of the STS for further analysis when studying STSms. To do this, first create an STS ROI for the entire length. In this example, the left STS in a 2009s parcellation:

 3dcalc -prefix {$ec}_L_STS_anat_parc -a0 L_parc+orig -expr "step(equals(a,513))"

Now, let's see how far anterior and posterior this ROI actually goes:

 3dclust 5 100 {$ec}_L_STS_anat_parc+orig

with output of

 #Volume  CM LR  CM PA  CM IS  minLR  maxLR  minPA  maxPA  minIS  maxIS    Mean     SEM    Max Int  MI LR  MI PA  MI IS
 #------  -----  -----  -----  -----  -----  -----  -----  -----  -----  -------  -------  -------  -----  -----  -----
   16516  -52.8  -22.7   -6.9  -66.0  -38.5  -74.2   33.1  -25.9   13.1        1        0        1  -44.0   27.6  -25.9 

We can determine the "cutoff" for y-values in the posterior segment:

 cutoff = ((maxPA-minPA) * 0.5) + minPA

In this example, cutoff = (33.1 + 74.2)*0.5-74.2 = -20.55.

Set this cutoff in the script:

 set cutoff = -20.55

Now, we can find voxels that are within the STS anatomical parcellation (in a 2009 parcellation in this case), active during auditory and visual blocks *and* are posterior to y = -20.55:

 3dcalc -prefix {$ec}_L_pSTS_auto_parc -a0 L_parc+orig -b17 {$ec}v1mr+orig -c20 {$ec}v1mr+orig -expr "step(b-2)*step(c-2)*step(equals(a,513))*step(y+$cutoff)"

 (I realize that the "step(y+$cutoff)" command is kind of backwards. If we want all the values more negative than -20.55, 
 I thought it would've been "isnegative(y-$cutoff)". But, this way works.)

Another example: here we find voxels within the left auditory cortex that are active (T > 2) during auditory blocks in a subject with an old-school 2005s parcellation:

 3dcalc -prefix {$ec}_L_aud_parc -a0 L_parc+orig -b17 {$ec}v1mr+orig -expr "step(b-2)*step(equals(a,33)+equals(a,81))"

---

Perhaps you would like to know the strength of activity averaged across these voxels. It's nice to use percent signal change values, since these are more easily comparable across subjects. First, use tent functions withing 3dDeconvolve to create a nice hemodynamic response for each voxel for each stimulus type:

 set v = 2
 3dDeconvolve -fout -tout -full_first -polort a -concat runs.txt \
 -input {$ec}Albl+orig -num_stimts 12 -nfirst 0 -jobs 2 \
 -mask {$ec}maskAlbl+orig  \
 -stim_times 1 McG.txt 'TENT(0,16,9)' -stim_label 1 McGurk \
 -stim_times 2 NonMcG.txt 'TENT(0,16,9)' -stim_label 2 InC \
 -stim_times 3 Cong.txt 'TENT(0,16,9)' -stim_label 3 Cong \
 -stim_times 4 Target.txt 'TENT(0,16,9)' -stim_label 4 Target \
 -stim_times 5 Ablock.txt 'TENT(0,26,14)' -stim_label 5 Ablock \
 -stim_times 6 Vemotion.txt 'TENT(0,26,14)' -stim_label 6 Vemotion \
 -stim_file 7 {$ec}rall_vr_motion.1D'[0]' -stim_base 7 \
 -stim_file 8 {$ec}rall_vr_motion.1D'[1]' -stim_base 8 \
 -stim_file 9 {$ec}rall_vr_motion.1D'[2]' -stim_base 9 \
 -stim_file 10 {$ec}rall_vr_motion.1D'[3]' -stim_base 10 \
 -stim_file 11 {$ec}rall_vr_motion.1D'[4]' -stim_base 11 \
 -stim_file 12 {$ec}rall_vr_motion.1D'[5]' -stim_base 12 \
 -iresp 1 irf1_tent -iresp 2 irf2_tent -iresp 3 irf3_tent -iresp 4 irf4_tent -iresp 5 irf5_tent -iresp 6 irf6_tent \
 -prefix {$ec}v{$v}mr

 3dTcat -prefix {$ec}v{$v}irf irf1_tent+orig irf2_tent+orig irf3_tent+orig irf4_tent+orig  irf5_tent+orig irf6_tent+orig

Then, convert these impulse response functions into percent signal change:

 3dcalc -a {$ec}v2irf+orig -b {$ec}EPIanatAlbl+orig -expr '100 * a/b * ispositive(10-a/b)' -prefix {$ec}psc_AV

Perhaps we are interested in the amplitude of response during each of four event-related conditions listed (McGurk, incongruent, congruent and target). The hemodynamic response peaks between 4-6 seconds, so it works to average the 3rd and 4th points together to estimate the height of the response in percent signal change:

 # Psc files-- takes average of points 2 and 3 of the HRF (the numbering is 0-8, so points 2 and 3 are the 3rd and 4th points of the series)
 3dMean -prefix {$ec}McGpsc {$ec}psc_AV+orig'[2]' {$ec}psc_AV+orig'[3]' 
 3dMean -prefix {$ec}InCpsc {$ec}psc_AV+orig'[11]' {$ec}psc_AV+orig'[12]' 
 3dMean -prefix {$ec}Congpsc {$ec}psc_AV+orig'[20]' {$ec}psc_AV+orig'[21]' 
 3dMean -prefix {$ec}Targetpsc {$ec}psc_AV+orig'[29]' {$ec}psc_AV+orig'[30]' 

Now, you can find the average percent signal change in the anatomical-functional ROI during each of the four conditions:

 3dROIstats -quiet  -mask  {$ec}_L_aud_parc+orig {$ec}McGpsc+orig
 3dROIstats -quiet  -mask  {$ec}_L_aud_parc+orig {$ec}InCpsc+orig
 3dROIstats -quiet  -mask  {$ec}_L_aud_parc+orig {$ec}Congpsc+orig
 3dROIstats -quiet  -mask  {$ec}_L_aud_parc+orig {$ec}Targetpsc+orig