I would like to use the brainstat library to relate cortical thickness to demographic measures. For that, I need a two-dimensional input matrix that represents thickness values for n_subjects * n_vertices. The columns should contain vertex measures for both hemispheres.
Freesurfer’s recon-all has already been run on the dataset, now my question is which files I should use to create this matrix and if some additional preprocessing steps have to be applied to these files?
The documentation of the output of recon-all does not describe all the files. It looks like the ?.thickness files could be used as input?
Can I use these files “out-of-the-box”? That is, can I simply concatenate the hemispheres for each subject and then concatenate all subjects? Or do I need to apply more preprocessing steps?
E.g., are they already in the same space? What if I would like to work with different smoothing options?
Typically if you are using multiple subjects you need to put them in a standard space, usually fsaverage. Also, typically srface morphometric values are smoothed (surface area with a larger FWHM compared to cortical thickness). Running recon-all with the qcache flag will automatically normalize the morphometrics and smooth them at various smoothing kernel sizes.
1.) I am not the maintainer of the dataset (the data has to be transferred from a collaborating group). That’s why I would be interested in retrieving the minimum set of files that would allow me to create the matrix myself. I saw that the -qcache option returns also non-smoothed but normalized versions (0mm). Do you think it would be smart to ask for these non-smoothed files and do the smoothing myself? This way the amount of data that has to be transferred could be kept small and I wouldn’t have to commit to a specific smoothing kernel. If yes, how would I do the smoothing? Could I use mris_smooth and I apply it to each of the ?h.thickness.fwhm0.fsaverage.mgh files?
2.) The maintainer already provided me the output from aparcstats2table for all atlases, measures, hemispheres & subjects. I was already able to create a 2D matrix from these files by choosing one of the atlases, average thickness as measure and concatenate the hemispheres. Similar question as above / sanity check: Can I use this matrix as an input matrix? Are the data already normalized here (i.e. can I compare values between subjects?). I cannot find any information in the documentation if the data has been normalized. The ?parc.stats files seem to be generated by the function mris_anatomical_stats but it is not described if this also includes normalization:
Parcel-wise data can be directly compared across subjects, as long as the same atlas is being applied to each subject. The spatial normalization is more important for the vertex-wise data, so you can be more sure a given vertex corresponds to the same position across subjects. Just make sure to use proper covariates as needed (e.g., control for total brain size if using parcel volume as a metric of interest).
Found out that the -qcache option calls the function mris_preproc behind the scenes. If you scroll down here to Example 3 this pretty much covers what would need to be done:
So in theory I could also ask the maintainer to give me the ?h.thickness files + a bunch of other files and do the normalization and (optional) smoothing myself. Problem: Just from looking at the generated log-file it’s hard to tell which files exactly are needed:
lh.thickness
lh.sphere.reg
lh.cortex.label
...?
So I think the solution above (asking for the unsmoothed but normalized versions, then manually smooth them) is the most easy solution. Here’s the full log-file, if someone is curious:
Wed Oct 9 17:04:50 CEST 2024
setenv SUBJECTS_DIR .
cd /zi/home/johannes.wiesner/work/testing/data_for_johannes
/zi/apps/opt/freesurfer_7.4.1/bin/mris_preproc --s sub-100000 --target fsaverage --hemi lh --meas thickness --fwhm 5 --out sub-100000-lh-thickness.sm5.mgh --SUBJECTS_DIR .
Linux zislrds0024.zi.local 5.15.0-116-generic #126-Ubuntu SMP Mon Jul 1 10:14:24 UTC 2024 x86_64 x86_64 x86_64 GNU/Linux
mris_preproc 7.4.1
freesurfer-linux-ubuntu22_x86_64-7.4.1-20230614-7eb8460
tmpdir is ./tmp.mris_preproc.3680819
Src lh sphere.reg
Trg lh sphere.reg
---------------------------------------------------
#@# 1/1 sub-100000 Wed Oct 9 17:04:50 CEST 2024 --------------
-rwxr-xr-x 1 johannes.wiesner domänen-benutzer 537347 Oct 7 13:34 ./sub-100000/surf/lh.thickness
-----------------------
mri_surf2surf --srcsubject sub-100000 --srchemi lh --srcsurfreg sphere.reg --trgsubject fsaverage --trghemi lh --trgsurfreg sphere.reg --tval ./tmp.mris_preproc.3680819/sub-100000.1.mgh --sval ./sub-100000/surf/lh.thickness --sfmt curv --noreshape --cortex
Source registration surface changed to sphere.reg
Target registration surface changed to sphere.reg
7.4.1
setenv SUBJECTS_DIR .
cd /zi/home/johannes.wiesner/work/testing/data_for_johannes
mri_surf2surf --srcsubject sub-100000 --srchemi lh --srcsurfreg sphere.reg --trgsubject fsaverage --trghemi lh --trgsurfreg sphere.reg --tval ./tmp.mris_preproc.3680819/sub-100000.1.mgh --sval ./sub-100000/surf/lh.thickness --sfmt curv --noreshape --cortex
sysname Linux
hostname zislrds0024.zi.local
machine x86_64
user johannes.wiesner
srcsubject = sub-100000
srcval = ./sub-100000/surf/lh.thickness
srctype = curv
trgsubject = fsaverage
trgval = ./tmp.mris_preproc.3680819/sub-100000.1.mgh
trgtype =
srcsurfreg = sphere.reg
trgsurfreg = sphere.reg
srchemi = lh
trghemi = lh
frame = 0
fwhm-in = 0
fwhm-out = 0
label-src = lh.cortex.label
label-trg = lh.cortex.label
OKToRevFaceOrder = 1
UseDualHemi = 0
Reading source surface reg ./sub-100000/surf/lh.sphere.reg
Loading source data
Reading curvature file ./sub-100000/surf/lh.thickness
Reading source subject label mask lh.cortex.label
masking the input
Reading target surface reg ./fsaverage/surf/lh.sphere.reg
Done
Using surf2surf_nnfr()
Mapping Source Volume onto Source Subject Surface
surf2surf_nnfr: building source hash (res=16).
Surf2Surf: Forward Loop (163842)
surf2surf_nnfr: building target hash (res=16).
Surf2Surf: Reverse Loop (134333)
Reverse Loop had 13574 hits
Surf2Surf: Dividing by number of hits (163842)
INFO: nSrcLost = 0
nTrg121 = 151627, nTrgMulti = 12215, MnTrgMultiHits = 2.11126
nSrc121 = 97519, nSrcLost = 0, nSrcMulti = 36814, MnSrcMultiHits = 2.17029
Reading target space mask label lh.cortex.label
masking the input
Saving target data
NOT Copying ctab from source
Saving to ./tmp.mris_preproc.3680819/sub-100000.1.mgh
-----------------------
mri_concat ./tmp.mris_preproc.3680819/sub-100000.1.mgh --o sub-100000-lh-thickness.sm5.mgh
ninputs = 1
Checking inputs
nframestot = 1
Allocing output
Done allocing
nframes = 1
Writing to sub-100000-lh-thickness.sm5.mgh
-----------------------
mri_surf2surf --hemi lh --s fsaverage --sval sub-100000-lh-thickness.sm5.mgh --tval sub-100000-lh-thickness.sm5.mgh --fwhm-trg 5 --noreshape --cortex
7.4.1
setenv SUBJECTS_DIR .
cd /zi/home/johannes.wiesner/work/testing/data_for_johannes
mri_surf2surf --hemi lh --s fsaverage --sval sub-100000-lh-thickness.sm5.mgh --tval sub-100000-lh-thickness.sm5.mgh --fwhm-trg 5 --noreshape --cortex
sysname Linux
hostname zislrds0024.zi.local
machine x86_64
user johannes.wiesner
srcsubject = fsaverage
srcval = sub-100000-lh-thickness.sm5.mgh
srctype =
trgsubject = fsaverage
trgval = sub-100000-lh-thickness.sm5.mgh
trgtype =
srcsurfreg = sphere.reg
trgsurfreg = sphere.reg
srchemi = lh
trghemi = lh
frame = 0
fwhm-in = 0
fwhm-out = 5
label-src = lh.cortex.label
label-trg = lh.cortex.label
OKToRevFaceOrder = 1
UseDualHemi = 0
Reading source surface reg ./fsaverage/surf/lh.sphere.reg
Loading source data
Reading source subject label mask lh.cortex.label
masking the input
INFO: trgsubject = srcsubject
Approximating gaussian smoothing of target with fwhm = 5.000000,
std = 2.123305, with 18 iterations of nearest-neighbor smoothing
Reading target space mask label lh.cortex.label
NN smoothing output with n = 18
Saving target data
NOT Copying ctab from source
Saving to sub-100000-lh-thickness.sm5.mgh
Cleaning up
-----------------------
rm -r ./tmp.mris_preproc.3680819
Wed Oct 9 17:05:01 CEST 2024
mris_preproc done
Could I use mris_smooth and I apply it to each of the ?h.thickness.fwhm0.fsaverage.mgh files?
I tried to use mris_smooth on a normalized file but got the following error:
mris_smooth sub-100000-lh-thickness.sm0.mgh output.mgh
error: No such file or directory
error: ERROR: MRISread: file 'sub-100000-lh-thickness.sm0.mgh' has 0 vertices!
Probably trying to use a scalar data file as a surface!
You have a few options. You can use the post-recon-all script, which runs qcache given a recon-all output, or try something like the following (from ChatGPT):
Resample the thickness and surface area data to fsaverage space:
Use the FreeSurfer command mri_surf2surf to resample cortical thickness and surface area maps from the subject’s native space to fsaverage space. You would typically have these maps in surf/lh.thickness, surf/rh.thickness, surf/lh.area, and surf/rh.area. Example for thickness (left hemisphere):
Repeat this for the right hemisphere and for surface area using the lh.area and rh.area files.
Smooth the resampled maps:
You can smooth the data using the mri_surf2surf command with the --fwhm option to apply a smoothing kernel.Example for thickness smoothing (left hemisphere):
This applies a 10 mm FWHM smoothing to the cortical thickness map. You can adjust the FWHM value as needed. Repeat this for the right hemisphere and for surface area maps.
TLDR; In order to re-create what the -qcache option is doing, you can run mris_preproc. For that you need each subjects ?h.thickness, ?h.sphere.reg and ?h.cortex.label files. You can also use mri_surf2surf which only needs the first two files but the results will be slightly different.
Long answer: Thanks @Steven! So I ran both mris_preproc and mri_surf2surf in order to normalize and smooth one left-hemispheric thickness file of the same subject. I can already say that the output is very similar. Most of the 163.842 vertex values are identical, the only difference is that mris_preproc will in rare cases assign a 0 to a vertex (would be curious to know why). Since mris_preproc is the function that the -qcache option uses I will probably stick to that. In order to run this function these files are needed:
?h.thickness
?h.sphere.reg
?h.cortex.label
In other words: This is the minimum set of files that is needed if you don’t have access to a full recon-all output.
This was the code I used for mris_preproc:
#!/bin/bash
set -e
# set path to data dir and change to it
data_dir=/zi/home/johannes.wiesner/work/testing/test_creating_thickness_maps/data/
cd $data_dir
# set subjects dir to data dir
export SUBJECTS_DIR=$data_dir
mris_preproc \
--s sub-100000 \
--target fsaverage \
--hemi lh \
--meas thickness \
--fwhm 10 \
--out mris_preproc.mgh
and for mri_surf2surf:
#!/bin/bash
set -e
# set path to data dir and change to it
data_dir=/zi/home/johannes.wiesner/work/testing/test_creating_thickness_maps/data/
cd $data_dir
# set subjects dir to data dir
export SUBJECTS_DIR=$data_dir
# NOTE: The fsaverage folder is treated like a subject, so it has to be part of SUBJECTS_DIR!
# You have to copy this folder into the dataset folder before running function
mri_surf2surf \
--srcsubject sub-100000 \
--trgsubject fsaverage \
--hemi lh \
--sval sub-100000/surf/lh.thickness \
--tval sub-100000/surf/mri_surf2surf.mgh \
--fwhm 10
I then loaded the output files from both scripts using nilearn.load_img and checked if the data arrays are equal.