CCN Wiki - User contributions [en]

Detrending FreeSurfer Data

2018-04-12T15:44:16Z

Connor:

Over the course of a run, there can be a linear drift in the signal in different regions of the brain. There are many possible causes for this that have nothing to do with any interesting aspect of your data -- in other words, this linear drift is a nuisance artifact. The second step is to remove this signal drift from the data because it can introduce spurious correlations between two unrelated time series. You can see this for yourself in a quick experiment you could whip up in Excel: take two vectors of 100 randomly generated numbers (e.g., randbetween(1,99)). They should be uncorrelated. Now add 1, 2, 3, ... , 99, 100 to the values in each vector. This simulates a linear trend in the data. You shouldn't be surprised to find that the two vectors are now highly and positively correlated!

A script has been written called detrend.sh that removes the linear trend in your BOLD data:
==== detrend.sh ====
#!/bin/bash
USAGE="Usage: detrend.sh filepattern surface sub1 ... subN"
if [ "$#" == "0" ]; then
echo "$USAGE"
exit 1
fi

#first parameter is the filepattern for the .nii.gz time series to be detrended, up to the surface indicator
#e.g., fmcpr.sm6.$surf.?h.nii.gz would use fmcpr.sm6.[self|fsaverage] as the filepattern
filepat="$1"
shift

#second parameter should be specified either as '''self''' or '''fsaverage'''
surf="$1"
shift

#after the shift command, all the arguments are shifted down one place and the first two arguments
#(the filepattern and surface)
#fall off the list. The remaining arguments should be subject_ids
subs=( "$@" );
hemis=( "lh" "rh" );

for sub in "${subs[@]}"; do
source_dir=${SUBJECTS_DIR}/${sub}/bold
echo ${source_dir}
if [ ! -d ${source_dir} ]; then
#The subject_id does not exist
echo "${source_dir} does not exist!"
else
cd ${source_dir}
readarray -t runs < runs
for r in "${runs[@]}"; do
if [ -n "${r}" ]; then
#the -n test makes sure that the run number is not an empty string
#caused by a trailing newline in the runs file
for hemi in "${hemis[@]}"; do
cd ${source_dir}/${r}
pwd
#subject_id does exist. Detrend
if [ "${surf}" = "self" ]; then

SURFTOUSE=${sub}
else
SURFTOUSE=fsaverage
fi
mri_glmfit --y ${source_dir}/${r}/${filepat}.${surf}.${hemi}.nii.gz \
--glmdir ${source_dir}/${r}/${hemi}.detrend \
--qa --save-yhat --eres-save \
--surf ${SURFTOUSE} ${hemi}
mv ${source_dir}/${r}/${hemi}.detrend/eres.mgh ${source_dir}/${r}/${filepat}.${hemi}.mgh
done
fi
done
fi
done

== Running the Script ==
Before running this script, you will need to create a text file called 'runs' in the bold/ directory for each subject's dataset, e.g.,
*FS_T1_501/
**bold/
***runs
***005/
***006/

The <code>runs</code> file simply lists each run folder on its own line:
005
006
The detrend.sh script uses this file to determine the folders containing the data to be detrended. If this file doesn't already exist, you can manually generate it in any text editor (e.g., <code>nano runs</code> or <code>gedit runs</code>), but the quickest method takes advantage of the fact that the run folders all start with 0, and uses common command-line utilities:

SUBJECT=FS_501
cd ${SUBJECTS_DIR}/${SUBJECT}/bold
#if there are fewer than 10 runs of BOLD data, then the run directories will probably have 2 leading zeros
ls -1 | grep "^00*" > runs

Assuming all your subject folders have the same run folders to detrend, you would detrend multiple subjects using detrend.sh, specifying a file pattern for the source data (i.e., the name of the preprocessed files generated by FS-FAST, omitting anything after the '?h' hemisphere identifier), followed by a list of subject IDs:
#A SPECIFIC EXAMPLE: (note these parameters may differ '''substantially''' from what you would be typing in)
detrend.sh fmcpr.sm6 self FS_T1_501 FS_T2_501 FS_T1_505 FS_T2_505
#For a more generalizable example of how you should call this function, see the section below using variables
The gist is that it calls the mri_glmfit function and saves the residuals after the linear trend has been removed from the data. Multiple files are generated in ?h.detrend/ directories in each run directory. The detrended data is subsequently copied back to the run directory as a new file called ${filepat}.?h.mgh, where ${filepat} is whatever file pattern you provided to the script (note that the source data are .nii.gz files, whereas the detrended data are .mgh files).

== Using variables ==
As I mentioned in the comments in the sample code snippet above, what you type depends entirely on the filenames, which in turn depend entirely on how the data were preprocessed. You can use environment variables to help walk you through figuring out the correct file patterns. Also, a handy shortcut exists if you happen to have a <code>subjects</code> file in <code>$SUBJECTS_DIR</code>. Putting these two techniques together:
FILEPATTERN=fmcpr #the preprocessed files will almost always be called '''fmcpr'''
SMOOTHING=".sm4" #how much smoothing did you use when you ran preproc-sess?
SLICETIME=".down" #[".up" | ".down" | OMIT ]
SURFACE=self #[self | fsaverage]

detrend.sh ${FILEPATTERN}${SLICETIME}${SMOOTHING} $SURFACE `cat ${subjects}`
This will execute the detrend.sh script on all the subjects listed in the subjects text file, using the subject's self surface. If the data were preprocessed with the fsaverage surface, you would specify ''fsaverage'', in place of ''self'' I just figured this out today.

[[Category: Time Series]]

Time Series Data in Surface Space

2018-04-12T15:43:37Z

Connor:

We can calculate mean time course vectors calculated across all voxels within regions defined in a FreeSurfer annotation (.annot) file (we can also extract the time courses on a vertex-by-vertex basis -- [[#Voxelwise_Time_Series | see below]]). Though any annotation file can be used for this purpose, we have been working at the scale of the Lausanne parcellation. A script exists, <code>gettimecourses.sh</code>, that will handle this:
=== gettimecourses.sh ===
Update Jan 11, 2018: I added a feature that uses the filepattern string to determine whether to use the self or fsaverage annotation.

#!/bin/bash
USAGE="Usage: gettimecourses.sh annot filepattern sub1 ... subN"

if [ "$#" == "0" ]; then
echo "$USAGE"
exit 1
fi

#first two parameters are is the annot files and filepatterns for the \
#.nii.gz time series to be detrended, up to the hemisphere indicator
#e.g., detrend.fmcpr.sm6.self.?h.nii.gz would use detrend.fmcpr.sm6.self as the filepattern
annot="$1"
shift
filepat="$1"
shift

#subjects
subs=( "$@" );
#hemispheres
hemis=( lh rh )

for sub in "${subs[@]}"; do
source_dir=${SUBJECTS_DIR}/${sub}/bold

#Which annotation? Look for "fsaverage" or "self" in the filepat
if [[ "${filepat}" =~ "fsaverage" ]]
then
SURF="fsaverage"
elif [[ "${filepat}" =~ "self" ]]
then
SURF=${sub}
fi

if [ ! -d ${source_dir} ]; then
#The subject_id does not exist
echo "${source_dir} does not exist!"
else
cd ${source_dir}
readarray -t runs < runs
for hemi in "${hemis[@]}"; do
for r in "${runs[@]}"; do
mri_segstats \
--annot ${SURF} ${hemi} ${annot} \
--i ${source_dir}/${r}/${filepat}.${hemi}.mgh \
--sum ${sub}_${hemi}_${annot}_${r}.${filepat}.sum.txt \
--avgwf ${sub}_${hemi}_${annot}_${r}.${filepat}.wav.txt
done
done
fi
done

== Running the Script ==
This script has two preconditions (conditions that must be satisfied prior to running the script)
# The first is that the data have been detrended following the [[Detrending_FreeSurfer_Data | procedure for detrending FreeSurfer data]].
# The second precondition is that there should be a file called 'runs' in the bold/ directory.
If you have followed the instructions for detrending your data, both conditions should already be met.

The script is run in a terminal by specifying the name of an annot file that can be found in the $SUBJECT_ID/label/ directory, a file pattern, followed by a list of subject IDs:
gettimecourses.sh lausanne fmcpr.siemens.sm6.self FS_T1_501
The above command would look in the /label directory for subject ''FS_T1_501'' for ''lh.lausanne.annot'' and ''rh.lausanne.annot'', and extract the mean time series for each defined region within each .annot file for the bold data found in the file ''fmcpr.sm6.self.?h.mgz'' in the run folders indicated in the ''runs'' file.
A series of output files is produced in the subject's bold/ directory. The time series files are named ${sub}_?h_${annot}_${run}.${filepattern}.wav.txt

==== A Trick ====
Similar to the shortcut described for detrending the data, you can take advantage of a <code>subjects</code> file in your SUBJECTS_DIR:

FILEPATTERN=fmcpr.siemens.sm4.self
ANNOT=lausanne
gettimecourses.sh $ANNOT $FILEPATTERN `cat ${subjects}`

=== Decoding the Regions ===
The time series output files are plain text files with rows=time points and columns=regions. In the case of the Lausanne 2008 parcellation, there are approximately 500 columns per file. Interpreting these data will require some sort of reference for the identity of each column.

As you might expect, there is a relationship between the column order and the segmentation ID, as confirmed by this archived email exchange involving yours truly: [https://mail.nmr.mgh.harvard.edu/pipermail/freesurfer/2012-December/026853.html Extracting resting state time series from surface space using Lausanne 2008 parcellation]

When you ran the gettimecourses.sh script, two files were created for each run folder: one containing the time series data (*.wav.txt) and one with some summary data (*.sum.txt). You can use the information from the 5th column of the .sum files to determine the segment labels corresponding to the columns of your time series data. A MATLAB function, '''[[ParseFSSegments.m]]''', has been written (ubfs /Scripts/Matlab folder) that will extract the data from these summary files.
%if called without providing filename will launch a dialog box to select the .sum.txt to parse
lh_summary=parseFSSegments();
%If you already know the filename, you can provide it as a parameter
rh_summary=parseFSSegments('FS_T1_501_lh_myaparc_60_005.fmcpr.sm6.sum.txt');

lh_region_names=lh_summary.segname;
rh_region_names=rh_summary.segname;

Note that the segmentation names will be consistent across participants for a given segmentation scheme. In practical terms, this means you '''should''' only have to do this once for a given .annot file, though it's probably worth double checking a few .sum.txt files to make sure the regions are listed in the same order because Murphy's law dictates that sometimes things get screwed up.

== Other Information ==
The FreeSurfer command mris_anatomical_stats uses an .annot file to query properties of a surface. The output of this program can be used to determine anatomical characteristics of each region. These values might be useful for modeling purposes.
SUBJECT_ID=FS_T1_501
ANNOT=lausanne
HEMIS=( "lh" "rh" )
for hemi in "${HEMIS[@]}"; do
mris_anatomical_stats -a ${SUBJECT_ID}/label/${hemi}.${ANNOT}.annot -f ${SUBJECTS_DIR}/${SUBJECT_ID}/${hemi}.${ANNOT}.stats.txt ${SUBJECT_ID} ${hemi}
done
The above BASH code snippet will report characteristics for subject FS_T1_501, using the ?h.lausanne.annot files in his/her labels/ directory. The for-loop is a concise way of ensuring that the lh and rh hemispheres are measured, with the results being written to ''lh.lausanne.stats.txt'' and ''rh.lausanne.stats.txt'' in his/her subject directory (i.e., in ''$SUBJECTS_DIR/FS_T1_501/'')

== Voxelwise Time Series ==
The above methods extract the mean time series across all the voxels contained in the regions indexed by the vertices in a given region or set of regions. Some techniques, such as MVPA, however, work on the single-voxel level. I [http://www.mail-archive.com/freesurfer%40nmr.mgh.harvard.edu/msg54885.html asked the Freesurfer mailing list] about this and got the following suggestion for using MATLAB to obtain 1 time series per vertex:
waveforms = fast_vol2mat(MRIread('waveform.nii.gz'));

== Now What? ==
Your time series data can be used to carry out a functional connectivity analysis (the [[ Functional_Connectivity_(Cross-Correlation_Method) | correlational ]] or neural network approach) or machine learning classification.
[[Category: Time Series]]

Highlight Reel

2018-02-19T16:40:16Z

Connor:

== Rebecca ==
*Came up with an amazing nickname for Chris "Norgie" McNorgan.

==Elizabeth==
*Agreed to coordinate the sorting, tagging and filtering of 20K concept features. That's a really huge quantity of data.

== James ==
*Got stumped with Matlab. Went and got a haircut. Came back ready to get things done.
*Made an awesome shell file that contains a bunch of code needed for functional analysis (which would be a pain to type in by hand).

== Mandy ==
*Reminded me that I had a package to pick up from the main office. They were document covers, so that was really exciting.

== Erica ==
*Painted original artwork for the lab.
**Also makes really nice looking powerpoint slides that can explain waveform filtering
*Whizzed through a whole pile of statistical analyses in R like "no big". So second-nature to her that she doesn't even realize her skill level.
*Got a coffee maker
*Always catches the little details

== Greg ==
*Made a funny remark about how writing a guide about writing guides was "pretty meta"
*Realizes we should have just been issuing electronic Amazon gift cards all this time
*Enabled login-controlled edits to the wiki

== Connor ==
*Showed up