Working with ROIs (Freesurfer)

The .label files created during autorecon3, or by saving an overlay are plaintext lists of vertices. It should be straightforward to find the intersection or union of the vertices in two or more .label files using common Unix shell commands.

.label Syntax

Below is the first few lines of the lh.banksts.label file for a participant. This file was produced during the Lausanne parcellation procedure.

#!ascii label  , from subject FS_0043 vox2ras=TkReg
1457
32992  -47.867  -62.891  33.007 0.0000000000
32993  -48.881  -62.907  32.963 0.0000000000
34014  -47.298  -61.696  33.538 0.0000000000
34015  -47.972  -61.259  33.648 0.0000000000
34024  -47.634  -62.264  33.251 0.0000000000
34025  -48.286  -61.775  33.303 0.0000000000
...etc

The first line is a header line. The second line indicates the number of vertices in the label file. The remaining lines indicate the vertex number, x,y,z coordinates of the vertex, and I have no idea what that last column indicates but it's not terribly important for our purposes.

Working From Individual-Level Results

Running selxavg3-sess will produce some number of statistical contrast maps that can be loaded as an overlay in tksurfer. Significant clusters can be turned into ROIs for that particular individual in a series of steps. This has the advantage of being applicable to a single dataset -- there is no need to wait until sufficient fMRI data has been collected to produce significant group-level contrasts. The downside is that this is process has a good number of steps.

Load Overlay

The first step is to load up the participant's surface data and the desired overlay. Care must be taken in selecting a contrast overlay that is appropriate for the type of analysis that you eventually plan to do. As a rule of thumb, the overlay contrast should be orthogonal (i.e., completely independent) of the contrast that you are ultimately interested in. For example, if your goal is to compare mean signal strength for low- vs. high-familiar trials, it would be inappropriate to use the contrast of high-familiar.vs.low-familiar to generate ROIs from which to calculate the mean signal strength. The reason for this is that this contrast will be selecting voxels that have already been found to demonstrate a significant difference. Naturally, any analysis of familiarity effects on these voxels will be biased! A more appropriate contrast to use would be something like task.vs.rest, which isn't biased towards either high or low-familiarity items, but is instead just showing the voxels that were more highly active during the task trials (which presumably contain equal numbers of high and low familiarity items).

After the overlay is loaded, you may want to play around with the overlay configuration parameters: for example you might choose to truncate your results to only show positive t-values, or modify your significance threshold.

Save Clusters as Labels

Now comes the somewhat tedious part of individually converting each cluster into a .label file. This simple but repetitive process is described here and you can get pretty quick at it:

1. If you wish to subdivide a large cluster, first create one or more paths that completely cross (e.g., bisect) the cluster.
2. Use the cursor to select a point anywhere in the desired cluster
3. Click the Custom Fill button. A new window will popup:
   • Choose Up to functional values below threshold
   • Optionally choose Up to and including paths if had subdivided the cluster. You will have to repeat these steps on the remainder of the cluster
   • The region will be colored in white before reverting to a yellow outlined region
4. On the Tools window, go to 'File > Label > Save Selected Label
   • A box will open with a path to save the ROI in the subject's label directory by default
   • Give it a meaningful name that includes the hemisphere and uses the .label filename extension (e.g. lh.TASK_001.label)
   • Click OK

Merge .label Files

If you wish to work with the resulting .label files on an individual level, your work is done. However, if you wish to merge multiple clusters into a single .annot file, you have to first create a color table which is called by mris_label2annot.

Create a CTAB File

Existing instructions for ctab files are a little hard to follow, so I'll do my best. The CTAB files are plaintext files with 6 columns and 1 row for each label. The structure of the columns is as follows:

Index

Label

R

G

B

A

Left and right hemispheres should have their own CTAB files, and the indices for each hemisphere should be numbered independently. For example, if there are 11 labels in the lh and 8 in the rh, then one ctab file should have indices 1 to 11 and the other should have indices 1 to 8.

The label names should correspond to the ?h.*.label files you created. For example, suppose you created the following 5 label files in the previous step:

lh.FAM_001.label
lh.FAM_002.label
rh.FAM_003.label
rh.FAM_004.label
rh.FAM_005.label

The lh ctab file should contain the following labels:

0    unknown    128    128    128    0
0    FAM_001    110    112    100    0
0    FAM_002    128    140    130    0

and the rh ctab file should contain the following labels: The lh ctab file should contain the following labels:

0    unknown    128    128    128    0
0    FAM_003    110    112    100    0
0    FAM_004    190    140    105    0
0    FAM_005    200    150    100    0

Sets of RGB values can be generated according to your whims, or you can use an external website to produce a set of n distinguishable colors. For example, I went to http://tools.medialab.sciences-po.fr/iwanthue/ to generate 34 distinct colors. Note that the final column, A, refers to the alpha transparency value, and is always set to 0 as far as I can tell.

Working From Group-Level Results

Group-level contrast maps, generated by Mri_glmfit mri_glmfit are well-suited for generating a set of ROIs that are appropriate for use across the set of participants entered into the group-level analysis. They have the useful property of being spatially consistent across all participants. There are a couple things to note, however. The first is that there must first be a sufficient number of participants to produce statistically significant contrast clusters. The second is that the group-level analysis must be done in fsaverage space, however you may need to map the ROIs back to a subject's native surface space.

Preparation of .label Files

The group analysis used the fsaverage surface, and since mris_divide_parcellation (which you may also be using, or use later) assumes we're subdividing a particular subject's annotation, it makes sense to copy the source .annot files into the fsaverage label/ folder (giving the appropriate ?h prefix)

ANNOT_FILE=cache.th30.abs.sign.ocn.annot
HEMI=rh 
cp ${HEMI}.${ANNOT_FILE}  $SUBJECTS_DIR\fsaverage\label\${HEMI}.${ANNOT_FILE}

The steps described below do not work on .annot files because they are not plaintext and easily manipulated. However, mri_annotation2label will convert a .annot file to a series of .label files.

e.g. mri_annotation2label --annotation ${ANNOT_FILE} --hemi ${HEMI} --subject fsaverage --surf pial --labelbase ${HEMI}.T1

Set Operations on .label Files

Intersecting .label Files (Intersection)

Given what we know about the contents of a label file, the unix comm command should suffice to find the common entries between two label files, which is the set intersection of vertices appearing in each.

https://mail.nmr.mgh.harvard.edu/pipermail//freesurfer/2009-August/011743.html

Merging .label Files (Union)

A union operation will find all vertices appearing in either .label file. You will want to make sure that only unique values are retained (i.e., you don't want to list the same vertex multiple times in the output .label file).

Or check out:

mri_mergelabels -i label1 -i label2 ... -o outputlabel 

This is a simple script that will merge two or more label files. It does this by catting the label files together (after removing the first two lines). It inserts a new header (the first two lines). The number of entries in the new file (ie, the number on the second line), is computed by summing those from the input files.

ROI 'VennDiagram' Matlab Function

ROIVennDiagram.m uses intersect and setdiff functions to extract intersecting and unique vertices. Variable names suggest that it's used for T1 and T2 data, but you can use it for any paired data. The label files generated above (using mri_annotation2label) need to be copied into a lh and rh directory within a parent dir. A simple modification would allow these files to remain in your fsaverage or other subject's dir, if you wanted to clutter those up.

function ROIvenndiagram = ROIvenndiagram(hemi)
%% function ROIvenndiagram = ROIvenndiagram('hemi');
%Pull out intersecting and different vertices for FS label files,
%and store in one of 3 text files (T1, T1T2, T2)
%
%Currently takes one string input indicating lh or rh.
%
%You'll need to adjust the file names on lines 42 and 43 to match your
%file names.
%
%Note - The headers in the output text files have quotes around them, as
%the headers begin with matlab unfriendly characters (#). So, you'll need to
%go into those 3 output text files and delete the quotes.
%
%Intended for use in parent dir containing lh and rh dirs, which
%contain label files.
% _____________________________________
% G.J. Smith
% gjsmith4@buffalo.edu
% _____________________________________
parentdir=pwd;
hemidir = sprintf('%s/%s', pwd, hemi);
cd(hemidir);

%% Pre-initialize array for cating labels
T1Agg = cell(1,5);
T2Agg = cell(1,5);

% Begin cycle through label files
for labeln = 1:3000 %arbirary number higher than any label number
    %% Apply appropriate number of 0's before number - eg. 001, 010, 100
    numcorrection=num2str(labeln);
    if labeln < 10
        labelnum=strcat('00',numcorrection);
    elseif labeln < 100
        labelnum=strcat('0',numcorrection);
    else
        labelnum=strcat(numcorrection);
    end
    
    %% T1 and T2 label file names -- Edit these to match your file naming convention
    T1_filename = sprintf('%s.T1-%s.label', hemi, labelnum);
    T2_filename = sprintf('%s.T2-%s.label', hemi, labelnum);
    
    %% Build one large aggregate file with all vertices from all label files
    header = 2; % lines 1 and 2 of label file stored seperately
    delimiter = ' ';
    if exist([pwd filesep T1_filename], 'file')
        %load data
        T1file = importdata(T1_filename,delimiter,header);
        %Store in one array
        T1Agg = vertcat(T1Agg, num2cell(T1file.data));
    end
    if exist([pwd filesep T2_filename], 'file')
        %load data
        T2file = importdata(T2_filename,delimiter,header);
        %Store in one array
        T2Agg = vertcat(T2Agg, num2cell(T2file.data));
    end
end

%% Check that label files were found
exist T1file;
if ans == 1
    display('Labels found, pulling out intersects and setdiffs');
else
    display('No label files found, check your dir or file names');
    return;
end

%% Find intersection and differences
T1Agg = cell2mat(T1Agg);
T2Agg = cell2mat(T2Agg);
%compare
%union = union(T1file.data(:,1), T2file.data(:,1));
intersection = intersect(T1Agg, T2Agg, 'rows'); %'rows' works with same # of columns
T1_diffs = setdiff(T1Agg, T2Agg, 'rows');
T2_diffs = setdiff(T2Agg, T1Agg, 'rows');

%% Build new label files
%T1
T1_differences = cell(50000, 5); %hard coded arbitrarily high number of rows for pre-init
T1_differences(1,1) = T1file.textdata(1,1); % header
T1_differences{2,1} = length(T1_diffs); % header line 2, new number of rows
T1_differences(3:length(T1_diffs) + 2,:) = num2cell(T1_diffs); % data

%T2
T2_differences = cell(50000, 5); %hard coded arbitrarily high number of rows for pre-init
T2_differences(1,1) = T2file.textdata(1,1); % header
T2_differences{2,1} = length(T2_diffs); % header line 2, new number of rows
T2_differences(3:length(T2_diffs) + 2,:) = num2cell(T2_diffs); % data

%Match
Match_labelfile = cell(50000, 5); %hard coded arbitrarily high number of rows for pre-init
Match_labelfile(1,1) = T1file.textdata(1,1); % header T1 and T2 for same subject is the same
Match_labelfile{2,1} = length(intersection); % header line 2, new number of rows
Match_labelfile(3:length(intersection) + 2, :) = num2cell(intersection); % data

%% Remove blank rows from oversized pre-alo arrays
T1_differences(all(cellfun(@isempty,T1_differences),2),:) = [];
T2_differences(all(cellfun(@isempty,T2_differences),2),:) = [];
Match_labelfile(all(cellfun(@isempty,Match_labelfile),2),:) = [];

%% Write files
cd(parentdir);
T1_fname = sprintf('%s.T1.label.txt', hemi);
writetable(cell2table(T1_differences), T1_fname,'Delimiter',' ','WriteVariableNames',false);

%T2
T2_fname = sprintf('%s.T2.label.txt', hemi);
writetable(cell2table(T2_differences), T2_fname,'Delimiter',' ','WriteVariableNames',false);

%Match
Match_fname = sprintf('%s.T1T2.label.txt', hemi);
writetable(cell2table(Match_labelfile), Match_fname,'Delimiter',' ','WriteVariableNames',false);

end

Visualize and Cluster (Intersected Label File)

• in progress*

Open tksurfer.

tksurfer fsaverage ?h inflated

Import ?h.aparc.annot files for a reference

Load the intersected label file.

Use cut line to make cluster boundaries (If necessary, otherwise see below).

Erase aparc.annot labels that overlap with a desired label.

Select area you want to cluster and use custom fill to create new label. Up to and including path, up to other labels, and up to unlabeled, filling from last clicked vertex. If you left adjacent aparc.annot labels these can be used as boundaries.

Save each label file using save selected label. -- must premake text files as cant specify file names. -- must be a better way of doing this, or getting export annotation to work.

Pro tip: Colour your label file red, when creating new labels, they will be blue. (see below)

CTAB and Annot

Now you'll need to create a ctab file for the labels to make a .annot file.

If you don't have a ton of labels, this can easily be done by hand by copying and editing the 'FreeSurferColorLUT.txt' file found in /usr/local/freesurfer, and simply renaming the regions(label names) to match your labelling convention. Then save it as a .annot.ctab file, e.g. 'RH.annot.ctab'.

Then you'll want to use the mris_label2annot command found here. Example:

mris_label2annot --l rh.label-001.txt --l rh.label-002.txt --l rh.label-003.txt --l rh.label-004.txt --l rh.label-005.txt --l rh.label-006.txt --l rh.label-007.txt --l rh.label-008.txt --l rh.label-009.txt --l rh.label-010.txt --l   rh.label-011.txt --l rh.label-012.txt --l rh.label-013.txt --l rh.label-014.txt --l rh.label-015.txt --l rh.label-016.txt --l rh.label-017.txt --l rh.label-018.txt --l rh.label-019.txt --l rh.label-020.txt --l rh.label-021.txt --l  rh.label-022.txt --l rh.label-023.txt --l rh.label-024.txt --l rh.label-025.txt --l rh.label-026.txt --ctab RH.annot.ctab --s fsaverage --h rh --annot T1T2intersected

Now that you have .annot files for your intersecting vertices, you can subparcellate and/or continue with network analyses.