Annotation Coordinates

When visualizing spatial relationships among network nodes representing regions in surface space, it is necessary to assign each node a position in 3D space. The XYZ coordinates of each vertex associated with a label can be found in the corresponding .label files in the label/regenerated_?h_* directories that were generated when carrying out the Lausanne parcellation. For example, to calculate representative spatial locations for the labels in the 60 label parcellation, the associated label files are in label/regenerated_?h_60.

The .label files are plaintext, and are formatted as follows:

#!ascii label  , from subject FS_T1_501 vox2ras=TkReg
1858
50139  -32.785  -24.134  19.744 0.0000000000
50148  -32.758  -24.176  19.440 0.0000000000
50149  -33.035  -24.591  18.914 0.0000000000
50157  -33.293  -24.616  18.023 0.0000000000
51159  -32.615  -23.587  20.632 0.0000000000
... etc.

After the header row and the number of vertices in the file, the left-most column is the vertex number, and the three following columns are the x,y,z coordinates for that vertex. A reasonable value for the representative spatial location of that node can be obtained by calculating the mean x,y,z coordinate.

But my coordinates are embedded in .annot files!

This will require a preparatory step: We must first decompile the .annot files into a series of .label files. Suppose we have lh and rh .annot files for some functionally-defined labels (e.g., from a group-wise GLM contrast), and the contrast results are saved in cache.th13.neg.sig.ocn.annot for each hemisphere (the exact names of the contrast annotation files will depend on the thresholding and tails applied to the group-level analysis).

The first thing you will need to do is copy them to the same temporary working directory, renaming them in the process (otherwise they will overwrite each other because the files have the same names!). In this example, I am going to rename my .annot files to functional_mask

mkdir ${SUBJECTS_DIR}/blobs
cp ${SUBJECTS_DIR}/RFX/my_analysis.lh/my_contrast/glm.wls/osgm/cache.th13.neg.sig.ocn.annot ${SUBJECTS_DIR}/blobs/lh.functional_mask.annot
cp ${SUBJECTS_DIR}/RFX/my_analysis.rh/my_contrast/glm.wls/osgm/cache.th13.neg.sig.ocn.annot ${SUBJECTS_DIR}/blobs/rh.functional_mask.annot

Now, we go into our temporary working directory to convert the annotation files into a series of .label files:

cd ${SUBJECTS_DIR}/blobs
ANNOTATION=functional_mask
for HEMI in ( lh rh ); do
  mri_annotation2label --subject fsaverage --hemi ${HEMI} --outdir ${HEMI} --annotation ${ANNOTATION} 
done

This will create lh and rh output directories in your current working directory (we use separate directories to make bookkeeping and the next step a little easier). In each of the output directories will be a series of .label files for the next step.

getSegCoords.m

A MATLAB script has been written to facilitate the calculation of mean segment coordinates. This script, getSegCoords.m can be found in the ubfs Scripts/Matlab folder. It takes a number of optional parameters, but if called without providing any parameters, you will be prompted to select specific files and folders that the script needs to complete its task.

[xyz, hemi, segnames]=getSegCoords

Called this way, without any parameters, the script will assume that the paired lh and rh files and directories have names that differ only with respect to the '?h' characters (by default, the automatch option is set to true). You will be prompted to first identify the label directory (i.e., the directory that contains the ?h.segment_id.label files) for the left hemisphere. The right hemisphere directory will be figured out automatically. The function then calls parseFSSegments.m for each hemisphere, which requires you to identify the '*.sum.txt' files that contain the individual segment names.

Source code for this MATLAB function can be found here.

vertex_key.m

Another MATLAB function that you might find useful for working with vertices and labels from your decompiled .annot files, vertex_key.m, will parse all the .label files in a given directory, and concatenate all the vertex information into a single .csv file. This file may be useful for finding the overlap between two sets of labels, so long as both sets of labels represent the same surface space (e.g., fsaverage). There may be some other useful applications, but i intend to use this to determine which Lausanne labels are overlapping with significant clusters from a GLM analysis. The MATLAB file should be in the ubfs/Scripts/Matlab folder, but here's the code as yet another programming example:

 function vertex_key(directory, outfilename)
% function vertex_key(DIRECTORY, OUTFILENAME)
% Optional parameters: 
%   DIRECTORY (string): 
%        The source directory to parse (default: pwd())
%   OUTFILENAME (string): the name of the output text file produced
%       (default: 'vertices.csv')
%
% No return values
%
% Reads through all .label files in current directory and gathers vertices
% and coordinates into a single csv text file with format:
% Vertex, X, Y, Z, hemi, labelname
%
% Useful for finding overlap or differences between sets of labels
%
% Sample usage:
% vertex_key('lh', 'lh_vertices.csv');
%

%Enforce default outfile name if none provided
if isempty(outfilename)
   outfilename='vertices.csv';
end

%Enforce default directory (PWD) if none provided
if isempty(directory)
   directory=pwd();
end

outfile=fopen(outfilename, 'w');
d=dir([directory filesep '*.label']);
flist={d.name};
for f=1:length(flist)
   %open the region label file
   fname=flist{f};
   fp=fopen(fullfile(directory,fname));
   %first 2 lines are junk
   tline=fgetl(fp);
   tline=fgetl(fp); %except this second line tells us how many vertices there are
   nvertices=str2double(tline);
   vdata=nan(nvertices,4);
   region=strsplit(fname,'.');
   hemi=region{1};
   region=region{2};
   ctr=0;
   while 1
       %read through each line of the .label file
       tline=fgetl(fp);
       %terminate loop if we hit the end
       if ~ischar(tline), break, end
       %otherwise, parse the line: vertex, x, y, z, nonsense we ignore
       ctr=ctr+1;
       lastline=tline;
       C=strsplit(lastline);
       for c=1:4
           %store the stuff we care about
           vdata(ctr,c)=str2double(C{c});
       end
   end
   fclose(fp);
   
   %write the vertex info for this region
   for n=1:nvertices
       fprintf(outfile, '%d,%.3f,%.3f,%.3f,%s,%s\n', ...
           vdata(n,1), vdata(n,2), vdata(n,3), vdata(n,4), ...
           hemi, region);
   end 
end
fclose(outfile);
end

Finding Label Overlap

Problem: I had a bunch of labels from the Lausanne parcellation, and I had a significance-thresholded volume from a group-level analysis. I wished to see which of the labels overlapped with the significant clusters. I decided to use a combination of MATLAB, FreeSurfer tools and a MySQL database:

1. Extract the vertices using vertex_key.m from the ?h.lausanne.annot files
2. The mri_glmfitsim output produces .annot files that cover the regions in the thresholded volumes. Extract the vertices from those files in the same way.
3. Import the .csv vertex information into two tables in MySQL (I'm obviously glossing over the steps of downloading and installing MySQL and MySQL Workbench on your computer, and having a basic understanding of SQL queries)
4. Use an INNER JOIN query to find all vertices in both tables, and I recommend using COUNT and GROUP BY to reduce the number of rows returned and get more useful results
5. Optionally determine how many vertices are present in the annotation of interest (in this example, how many vertices are present in each of the lausanne labels)
   • You can use these values to determine the proportion of overlap. In my experiment, I chose to keep only regions where 0.75 or more of the total number of vertices for a label were shared with the significant cluster labels.

Cool. Now what? Well, suppose you've identified a list of labels that you want to keep because they overlap sufficiently with your significant regions, and you've got them in a MATLAB cell array named keylabels.

keylabels{1}
ans=
      'lh.superiorfrontal_27'

Now you might want to use these labels to determine the indices (e.g., rows, columns) that you want to keep for some sort of analysis (e.g., reduce an adjacency matrix to include only those regions). The MATLAB ismember function will help here:

nodes_to_keep=find(ismember(lausanne_labels, keylabels));

This will return a vector of the node indices where the lausanne_label name matches the keylabel name. These are the node indices that you will want to retain:

ReducedAdjacency=ADJ(nodes_to_keep, nodes_to_keep); %The 1000 x 1000 ADJ matrix is now a 255 x 255 matrix!