-
Notifications
You must be signed in to change notification settings - Fork 6
Mapping Anatomical Correlations Across Cerebral Cortex (MACACC Method)
This tutorial is for applying Jason Lerch's method "MACACC" on any cortical thickness data you have. If you use it or want more details on the technique, cite:
Lerch, J. P., Worsley, K., Shaw, W. P., Greenstein, D. K., Lenroot, R. K., Giedd, J., & Evans, A. C. (2006). Mapping anatomical correlations across cerebral cortex (MACACC) using cortical thickness from MRI. Neuroimage, 31(3), 993-1003. Link
In short, MACACC is used to find cortical thickness correlations of a seed vertex with whole brain vertices. This means that we will correlate the cortical thickness (in mm) from a vertex against all vertices in the brain and find if there are significant similarities with local and non-local cortical areas.
The technique is simple in the application, the problem is finding a right seed to study and the conclusions you draw from it, and this will depend on your scientific questions and hypotheses.
First you need to run the Cortical Thickness Pipeline using CIVET and have the database in R. This is explained in another subject in this wiki.
This part is the most important. You shouldn't do this at random or correlate everything and find some map that fits your story. I recommend having a seed of interest based on your CT results (i.e. significant group differences), based on the bibliography or based on meta-analysis websites such as Neurosynth. There is, however, a way to do this without a prior hypothesis (see below MACACC-Strength).
Once you have a region of interest, localize the vertex point using Display (minc-tookit) and write down the specific number (i.e. 13088). You can write it in R as a variable seed_point=13089. Notice it's the same number + 1 because R starts at 1 while Display starts at 0. Remember you can have more than one seed per brain, sometimes you want the same point for both hemispheres, or sometimes you are interested in more than one seed.
'Insert Figure 1'
Create a vertexTable for each the right and left hemispheres cortical thickness if it doesn't exist already.
R > lh_ct_table <- vertexTable(left_ct)
Then, make a variable (column) showing the CT value at the seed vertex in each subject (row).
R > dataset$seed_vertex <- lh_ct_table[seed_point,]
Make sure you have the RMINC package loaded. In R, the equivalent for Pearson's Correlation is the use of a Linear Model. The formula looks like this:
R > macacc_corr <- vertexLm(left_ct ~ seed_vertex, dataset)
Where left_ct is the left hemisphere cortical thickness variable and seed_vertex is the seed variable.
It is recommended to use the FDR-corrected values, therefore:
R > vertexFDR(macacc_corr)
Example of FDR output
The important value here is the t-value at FDR 0.05. You will then have t-values but not correlation coefficients. If you want that, look ahead.
As with the CIVET CT, you need to first export the vertex values and you can plot them in Display or brain-view2.
write.table(macacc_corr ['t-value'], file="macacc_FDR.txt")
Threshold the t-value and look at your correlated areas.
'Insert Figure 2'
You may want to show correlation coefficient maps apart from the t-values. To do this, you need to extract convert the r^2 values to r using this forumla:
macacc_r_corr = sqrt(macacc_corr[,"R-squared"])*sign(macacc_corr[,"tvalue-lv1"])
Then, write the table
write.table(macacc_r_corr, col.names=FALSE, row.names=FALSE, file="macacc_r_corr.txt")
The threshold here will be 1 and -1.
'Insert Figure 3'
First, when you write.tableit will create the text files with t-values that will go beyond > 100 sometimes due to the autocorrelation in the seed. You need to change the text files by hand just by finding the seed row and using the maximum t-value in the data (i.e. 70) or changing it to 0 so it appears empty.
Second, once you've changed the files you need to change the civet_create.sh script to create MACACC maps. Gabriel has created a script for creating nice maps using CIVET, but it cannot be used directly to create maps with MACACC. For that you need to use the same script with a little modification. The script will create a gradient depending on the min and max values, which could make your maps look faded if the t-value is 80 for example. Therefore, we need to rescale the maximum to value in the new script.
Remember for this script you need to run either on Scinet, the CIC or the Virtual Machine. Also, in this script you need to modify the location of the file CIVET-CC-mask.txt (can be downloaded from this Cobralab github, in the section minc-toolkit-extras).
This create_civet.sh file looks for the maximum t-value but it will always be 10 for both hemispheres. You can change the maximum value to whatever you need.
#!/bin/bash
#USAGE
#create_civet_image.sh left.obj left_statmap left_column left_FDR right.obj right_statmap right_column right_FDR output.png
#Automatically generates a standardized view of a brain, using stats files for left and right hemispheres (from RMINC)
#Colourizes objects using stats files, ray traces standard views and them merges them together
# To change the maximum (for example in MACACC), change this line $(sort -g $leftstatmap | tail -1) to the MAX t-value.
set -euo pipefail
IFS=$'\n\t'
if [[ $# -ne 9 ]]
then
echo Usage: create_ciet_image.sh left_average.obj left_statmap left_column left_FDR right_average.obj right_statmap right_column right_FDR output.png
exit
fi
TMPDIR=$(mktemp -d)
leftbrain=$1
leftstatmap=$2
leftcolumn=$3
leftFDR=$4
rightbrain=$5
rightstatmap=$6
rightcolumn=$7
rightFDR=$8
output=$9
#Pick out column, multiply by mask and store in new tempfile
cut -d " " -f $leftcolumn $leftstatmap | paste -d\* /media/sf_INP_MRI_Backup/egarza_stuff/CIVET-CC-mask.txt - | sed 's/[eE]+\{0,1\}/*10^/g' | bc > $TMPDIR/$(basename $leftstatmap)
cut -d " " -f $rightcolumn $rightstatmap | paste -d\* /media/sf_INP_MRI_Backup/egarza_stuff/CIVET-CC-mask.txt - | sed 's/[eE]+\{0,1\}/*10^/g' | bc > $TMPDIR/$(basename $rightstatmap)
#Redefine statmap to masked tempfile
leftstatmap=$TMPDIR/$(basename $leftstatmap)
rightstatmap=$TMPDIR/$(basename $rightstatmap)
#Colourize left
if (( $(echo "$leftFDR <= $(sort -g $leftstatmap | tail -1)" | bc) ))
then
colour_object $leftbrain $leftstatmap $TMPDIR/left_pos.obj red_metal_inv $leftFDR 8 transparent black composite
else
cp $leftbrain $TMPDIR/left_pos.obj
fi
if (( $(echo "-$leftFDR >= $(sort -g $leftstatmap | head -1)" | bc) ))
then
colour_object $TMPDIR/left_pos.obj $leftstatmap $TMPDIR/left_comb.obj cold_metal $(sort -g $leftstatmap | head -1) -$leftFDR black transparent composite
else
cp $TMPDIR/left_pos.obj $TMPDIR/left_comb.obj
fi
#Colourize right
if (( $(echo "$rightFDR <= $(sort -g $rightstatmap | tail -1)" | bc) ))
then
colour_object $rightbrain $rightstatmap $TMPDIR/right_pos.obj red_metal_inv $rightFDR 8 transparent black composite
else
cp $rightbrain $TMPDIR/right_pos.obj
fi
if (( $(echo "-$rightFDR >= $(sort -g $rightstatmap | head -1)" | bc) ))
then
colour_object $TMPDIR/right_pos.obj $rightstatmap $TMPDIR/right_comb.obj cold_metal $(sort -g $rightstatmap | head -1) -$rightFDR black transparent composite
else
cp $TMPDIR/right_pos.obj $TMPDIR/right_comb.obj
fi
#Generate views
echo """ray_trace -output $TMPDIR/bottom.rgb $TMPDIR/left_comb.obj $TMPDIR/right_comb.obj -size 1000 1000 -crop -sup 3 -bg black -bottom
ray_trace -output $TMPDIR/top.rgb $TMPDIR/left_comb.obj $TMPDIR/right_comb.obj -size 1000 1000 -crop -sup 3 -bg black -top
ray_trace -output $TMPDIR/front.rgb $TMPDIR/left_comb.obj $TMPDIR/right_comb.obj -size 1000 1000 -crop -sup 3 -bg black -front
ray_trace -output $TMPDIR/back.rgb $TMPDIR/left_comb.obj $TMPDIR/right_comb.obj -size 1000 1000 -crop -sup 3 -bg black -back
ray_trace -output $TMPDIR/left_medial.rgb $TMPDIR/left_comb.obj -size 1000 1000 -crop -sup 3 -bg black -right
ray_trace -output $TMPDIR/left_lateral.rgb $TMPDIR/left_comb.obj -size 1000 1000 -crop -sup 3 -bg black -left
ray_trace -output $TMPDIR/right_medial.rgb $TMPDIR/right_comb.obj -size 1000 1000 -crop -sup 3 -bg black -left
ray_trace -output $TMPDIR/right_lateral.rgb $TMPDIR/right_comb.obj -size 1000 1000 -crop -sup 3 -bg black -right""" | parallel -j8
#Create pngs
for file in $TMPDIR/*.rgb
do
echo convert $file $TMPDIR/$(basename $file rgb)png
done | parallel -j8
#Create colourbars
convert xc:black xc:rgb\(0,0,128\) xc:rgb\(0,128,255\) xc:rgb\(128,255,255\) xc:white -append -filter Cubic -resize 30x600\! -bordercolor black -border 2x2 -bordercolor white -border 2x2 $TMPDIR/coldmetalinv.png
convert xc:black xc:rgb\(128,0,0\) xc:rgb\(255,0,128\) xc:rgb\(255,128,255\) xc:white -append -filter Cubic -resize 30x600\! -bordercolor black -border 2x2 -bordercolor white -border 2x2 $TMPDIR/redmetalinv.png
#Create combined figure
convert -gravity Center -background black +append $TMPDIR/left_medial.png $TMPDIR/right_medial.png $TMPDIR/1.png
convert -gravity Center -background black +append $TMPDIR/left_lateral.png $TMPDIR/right_lateral.png $TMPDIR/2.png
convert -gravity Center -background black +append $TMPDIR/top.png $TMPDIR/bottom.png $TMPDIR/3.png
convert -gravity Center -background black +append $TMPDIR/front.png $TMPDIR/back.png $TMPDIR/4.png
convert $TMPDIR/1.png $TMPDIR/2.png $TMPDIR/3.png $TMPDIR/4.png -gravity Center -background black -append $TMPDIR/unlabelled.png
convert $TMPDIR/unlabelled.png \
-page +1364+1685 $TMPDIR/coldmetalinv.png \
-page +199+1685 $TMPDIR/redmetalinv.png \
-stroke none -fill white -pointsize 45 -annotate +1410+2285 'FDR 10%' \
-stroke none -fill white -pointsize 45 -annotate +1410+1725 'FDR<1%' \
-stroke none -fill white -pointsize 45 -annotate +5+2285 'FDR 10%' \
-stroke none -fill white -pointsize 45 -annotate +5+1725 'FDR<1%' \
-flatten $output
The command in this case includes the stats maps and the correct FDR values for your threshold (i.e. 10%).
./civet_create.sh lh_average.obj left_stat_map.txt 1 1.92 rh_average.obj right_stat_map.txt 1 2.95 imagefile.ppm
Now, if you want to create correlation maps with the Pearson's R text files you output, then you need this other script which will constrain the min and max to -1 1 respectively, and it will have different colors.
#!/bin/bash
#USAGE
#create_civet_image.sh left.obj left_statmap left_column left_FDR right.obj right_statmap right_column right_FDR output.png
#Automatically generates a standardized view of a brain, using stats files for left and right hemispheres (from RMINC)
#Colourizes objects using stats files, ray traces standard views and them merges them together
# To change the maximum (for example in MACACC), change this line $(sort -g $leftstatmap | tail -1) to the MAX t-value.
set -euo pipefail
IFS=$'\n\t'
if [[ $# -ne 9 ]]
then
echo Usage: create_ciet_image.sh left_average.obj left_statmap left_column left_FDR right_average.obj right_statmap right_column right_FDR output.png
exit
fi
TMPDIR=$(mktemp -d)
leftbrain=$1
leftstatmap=$2
leftcolumn=$3
leftFDR=$4
rightbrain=$5
rightstatmap=$6
rightcolumn=$7
rightFDR=$8
output=$9
#Pick out column, multiply by mask and store in new tempfile
cut -d " " -f $leftcolumn $leftstatmap | paste -d\* /media/sf_INP_MRI_Backup/egarza_stuff/CIVET-CC-mask.txt - | sed 's/[eE]+\{0,1\}/*10^/g' | bc > $TMPDIR/$(basename $leftstatmap)
cut -d " " -f $rightcolumn $rightstatmap | paste -d\* /media/sf_INP_MRI_Backup/egarza_stuff/CIVET-CC-mask.txt - | sed 's/[eE]+\{0,1\}/*10^/g' | bc > $TMPDIR/$(basename $rightstatmap)
#Redefine statmap to masked tempfile
leftstatmap=$TMPDIR/$(basename $leftstatmap)
rightstatmap=$TMPDIR/$(basename $rightstatmap)
#Colourize left
colour_object $leftbrain $leftstatmap $TMPDIR/left_pos.obj spectral -1 1 transparent black composite
cp $TMPDIR/left_pos.obj $TMPDIR/left_comb.obj
#Colourize right
colour_object $rightbrain $rightstatmap $TMPDIR/right_pos.obj spectral -1 1 transparent black composite
cp $TMPDIR/right_pos.obj $TMPDIR/right_comb.obj
#Generate views
echo """ray_trace -output $TMPDIR/bottom.rgb $TMPDIR/left_comb.obj $TMPDIR/right_comb.obj -size 1000 1000 -crop -sup 3 -bg black -bottom
ray_trace -output $TMPDIR/top.rgb $TMPDIR/left_comb.obj $TMPDIR/right_comb.obj -size 1000 1000 -crop -sup 3 -bg black -top
ray_trace -output $TMPDIR/front.rgb $TMPDIR/left_comb.obj $TMPDIR/right_comb.obj -size 1000 1000 -crop -sup 3 -bg black -front
ray_trace -output $TMPDIR/back.rgb $TMPDIR/left_comb.obj $TMPDIR/right_comb.obj -size 1000 1000 -crop -sup 3 -bg black -back
ray_trace -output $TMPDIR/left_medial.rgb $TMPDIR/left_comb.obj -size 1000 1000 -crop -sup 3 -bg black -right
ray_trace -output $TMPDIR/left_lateral.rgb $TMPDIR/left_comb.obj -size 1000 1000 -crop -sup 3 -bg black -left
ray_trace -output $TMPDIR/right_medial.rgb $TMPDIR/right_comb.obj -size 1000 1000 -crop -sup 3 -bg black -left
ray_trace -output $TMPDIR/right_lateral.rgb $TMPDIR/right_comb.obj -size 1000 1000 -crop -sup 3 -bg black -right""" | parallel -j8
#Create pngs
for file in $TMPDIR/*.rgb
do
echo convert $file $TMPDIR/$(basename $file rgb)png
done | parallel -j8
#Create colourbars
convert xc:black xc:rgb\(0,0,128\) xc:rgb\(0,128,255\) xc:rgb\(128,255,255\) xc:white -append -filter Cubic -resize 30x600\! -bordercolor black -border 2x2 -bordercolor white -border 2x2 $TMPDIR/spectral.png
convert xc:black xc:rgb\(128,0,0\) xc:rgb\(255,0,128\) xc:rgb\(255,128,255\) xc:white -append -filter Cubic -resize 30x600\! -bordercolor black -border 2x2 -bordercolor white -border 2x2 $TMPDIR/spectral.png
#Create combined figure
convert -gravity Center -background black +append $TMPDIR/left_medial.png $TMPDIR/right_medial.png $TMPDIR/1.png
convert -gravity Center -background black +append $TMPDIR/left_lateral.png $TMPDIR/right_lateral.png $TMPDIR/2.png
convert -gravity Center -background black +append $TMPDIR/top.png $TMPDIR/bottom.png $TMPDIR/3.png
convert -gravity Center -background black +append $TMPDIR/front.png $TMPDIR/back.png $TMPDIR/4.png
convert $TMPDIR/1.png $TMPDIR/2.png $TMPDIR/3.png $TMPDIR/4.png -gravity Center -background black -append $TMPDIR/unlabelled.png
convert $TMPDIR/unlabelled.png \
-page +1364+1685 $TMPDIR/spectral.png \
-page +199+1685 $TMPDIR/spectral.png \
-stroke none -fill white -pointsize 45 -annotate +1410+2285 'r +1' \
-stroke none -fill white -pointsize 45 -annotate +1410+1725 'r -1' \
-stroke none -fill white -pointsize 45 -annotate +5+2285 'r +1' \
-stroke none -fill white -pointsize 45 -annotate +5+1725 'r -1' \
-flatten $output
In this case, the command should be:
./civet_create.sh lh_average.obj left_corr_map.txt 1 0 rh_average.obj right_corr_map.txt 1 0 imagefile.ppm
If you have 2 groups, you can get MACACC maps for each one separately. There is a formula to compare them using interaction (MACACC-Slope), but I haven't seen used in papers yet aside from the core paper.
# Original MACACC model
macacc_lnacc_l <- vertexLm(ct_left ~ l_nacc_vs, data = ct_data_macacc)
# Model with interaction.
macacc_lnacc_l_slope <- vertexLm(ct_left ~ l_nacc_vs + group + l_nacc_vs*group, data = ct_data_macacc)
# vertexFDR(macacc_lnacc_l_slope)
If the interaction term is significant, you have a difference in slopes.