Surface parcellation in SUMA


I am trying to convert a publicly-available surface parcellation to standard mesh in order to extract mean timeseries for connectivity analysis. My current method does not work as expected (explained below), but I’m also wondering if there is a more straight forward implementation.

I am using 3dROIstats to extract the timeseries, but the number of columns in the resulting files (~30) does not equal the number of regions that should be there (~200). Also, the indices of the columns look weird (eg Mean_-31930). Like, how is there a negative index?

So I am confused and wondering how 3dROIstats labels the regions from a dataset derived using FSread_annot. It seems like the .cmap may give me some pointers, but the label numbers are huge and make no sense to me. In the .cmap generated from FSread_annot, the comment line has (surf annotation), but none of those numbers are negative.

Sorry for the wall of text below, but I wanted to also see if someone can see a simpler solution or catch a mistake I made.



The parcellation can be found here.
Specifically, I’m interested in using /Parcellations/Freesurfer5.3/fsaverage6//label/?h.Schaefer2018_400Parcels_17Networks_order.annot.

As input into this workflow, I have:

  1. resting state data preprocessed in the volume
  2. freesurfer data on each subject
  3. a freesurfer distribution
  4. the .annot files from the above repo

I specify the freesurfer distribution as different from the .annot files because @SUMA_Make_Spec_FS will not work solely on the fsaverage6 directory in the parcellation distribution (because only the surf and label directories are supplied). The parcellation was derived on fsaversage6, so I’m starting there.

Here is my workflow:

  1. Generate SUMA-readable mesh from the fsaverage6 directory supplied in the freesurfer distribution, specifying only -ld 60 because that’s all I need.

cd ${dir}/fsaverage6
@SUMA_Make_Spec_FS -sid fsaverage6 -ld 60 -NIFTI

  1. Convert the .annot files in the parcellation to SUMA-readable formats, then resample to std.60 mesh.

for hh in "lh" "rh"
	# convert the .annot to SUMA
	FSread_annot \
		-input $dir/Parcellations/FreeSurfer5.3/fsaverage6/label/${hh}.${dset}.annot \
		-hemi ${hh} \
		-FScmap $dir/Parcellations/roject_to_individual/${dset}_LUT.txt \
		-col_1D ./${dset}.${hh}.col \
		-roi_1D ./${dset}.${hh}.roi \
		-cmap_1D ./${dset}.${hh}.cmap \
		-dset ./${dset}.${hh}.niml.dset

	# resample the SUMAarized .annot to std.60
	SurfToSurf \
		-i_gii std.60.${hh}.smoothwm.gii \
		-i_gii ${hh}.smoothwm.gii \
		-prefix std.60.${hh} \
		-mapfile std.60.fsaverage6_${hh}.niml.M2M \
		-dset ${dset}.${hh}.niml.dset \
		-output_params NearestNode

Among the intermediary files, these steps result in std.60.?h.Schaefer2018_400Parcels_17Networks_order.?h.niml.dset, which I can open in SUMA and they look reasonable.
Side note: I couldn’t get SurfToSurf to work without having to have filenames that specify the hemisphere twice. One I add manually using -prefix, and another the program adds. The SurfToSurf help file says that only .1D outputs are available, but this method gives me .niml.dsets as well.

  1. Generate SUMA-readable mesh for each subject

cd $dir/sub-01
@SUMA_Make_Spec_FS -sid sub-01 -ld 60 -NIFTI

  1. Project the functional data to the surface

for hh in "lh" "rh"
	3dvol2surf \
		-spec ./SUMA/std.60.sub-01_${hh}.spec \
		-surf_A smoothwm \
		-surf_B pial \
		-sv ./SUMA/sub-01_SurfVol.nii \
		-grid_parent ./func/filtered_func_data.nii.gz \
		-map_func ave \
		-f_steps 15 \
		-f_index nodes \
		-outcols_NSD_format \
		-out_niml ./func/vol2surf_func.${hh}.niml.dset

  1. Extract the mean timeseries from the parcellation

for hh in "lh" "rh"
	3dROIstats \
		-nobriklab \
		-1DRformat \
		-mask ./SUMA/std.60.${hh}.Schaefer2018_400Parcels_17Networks_order.${hh}.niml.dset \
		./func/vol2surf_func.${hh}.niml.dset > ./func/Schaefer2018_400Parcels.ts.${hh}.1D

As an addendum:

I should be getting 200 regions per hemisphere, but instead 3dROIstats is giving me 39.

I have attached an image of the ROI file. It looks like 3dROIstats might be averaging over the colored regions, but, under the hood each of those regions are further subdivded.

Update / consolidation on this issue:

I believe the issue is with how 3dROIstats is handling surface annotation datasets (or this surface annotation in particular).

I have checked the dataset created from FSread_annot, and everything looks to be in order. The regions are colorized as they should, and each node has a label that corresponds to freesurfer’s “annotation value”: (B * 256^2) + (G * 256) + (R) (which is why the values are so big).

While the image looks like a small number of regions relative to 200 per hemi, the LUT varies the colors but 1 at a time, thus giving each of the 200 regions a unique annotation value.

However, the output from 3dROIstats is still confusing. I am only getting 39 columns in the left hemisphere, whereas I should be getting 201. If 3dROIstats was averaging over the regions that are similarly colored in the image in my previous post then I should be getting 57 columns of timeseries.

The expected behavior of 3dROIstats, based on my previous experience, is that if you have a region with a label (eg 5193422 - a surface annotation value corresponding to r: 206 g: 62 b: 79), then 3dROIstats would have a corresponding Mean_5193422 column.

Also, the output labels of 3dROIstats are still confusing (eg Mean_XXXX). Here are the column names I’m getting:

Mean_-32442  	Mean_-32186  	Mean_-31930  	Mean_-30854  	Mean_-29577  	Mean_-27418  	Mean_-27417  	Mean_-27162  	Mean_-26625  	Mean_-26369  	Mean_-26113  	Mean_-25782  	Mean_-25526  	Mean_-13269  	Mean_-13013  	Mean_-1828  	Mean_-257  	Mean_-2  	Mean_-1  	Mean_256  	Mean_257  	Mean_511  	Mean_17Networks_LH_VisPeri_ExStrSup_1  	Mean_4728  	Mean_4984  	Mean_10505  	Mean_12300  	Mean_12935  	Mean_12936  	Mean_13191  	Mean_13505  	Mean_15044  	Mean_15045  	Mean_15300  	Mean_16077  	Mean_16078  	Mean_16333  	Mean_30208  	Mean_30464

Of note:

  • none of those numbers are annotation or node label values.
  • Many of the numbers are negative
  • There is a weird random string in there (Mean_17Networks_LH_VisPeri_ExStrSup_1)
  • I have no idea how 39 regions corresponds to the annotation file

This might have something to do with line 587 in 3dROIstats

AFNI_get_dset_val_label(mask_dset, (double)(i - 32768), sklab);

I’ll keep looking, but if anyone has thoughts let me know.

Further update / solution, in case there are others trying to convert a Freesurfer annotation to SUMA surfaces.

The issue arose because Freesurfer’s annotation values are so massive (almost all are greater than 32767). I used 3dRank to lower the annotation values to something more reasonable, which gave me the expected 3dROIstats output.

Thanks to Pete Molfese for the suggestion!


for hh in "lh" "rh"
	3dRank \
		-prefix std.60.Schaefer2018_400Parcels_17Networks_ranked.${hh}.niml.dset \
		-input std.60.${hh}.Schaefer2018_400Parcels_17Networks_order.${hh}.niml.dset

I am reviving this thread because 3dROIstats is behaving weirdly using my surface parcellation. I think I’ve found a solution, but also maybe a bug?

(In narrator voice) Previously on this thread:

I converted a Freesurfer annotation to .niml.dset using FS_read_annot, then I resample it to -ld 60 and ranked the node values so that 3dROIstats could handle the parcel labels.

I am now extracting mean time series from each parcel using time series projected to the surface and the results from 3dRank as the mask.

I have examined the .rankmap.1D file and the mapping from rank to annotation values seems to be in order and I have the correct number of labels. However, when extract the time series using the ranked parcellation as a mask, I’m getting weird results.

What I expect: I have 201 labels, 200 parcels of interest + the medial wall for each hemisphere. I was expecting 201 columns of timeseries from 3dROIstats with labels of Mean_0 to Mean_200, which should correspond to the rank values in .rankmap.1D

What I’m getting: 200 columns with the following values:

Mean_1  	Mean_2  	Mean_3  	Mean_4  	Mean_5  	Mean_6  	Mean_7  	Mean_8  	Mean_9  	Mean_10  	Mean_11  	Mean_12  	Mean_13  	Mean_14  	Mean_15  	Mean_16  	Mean_17  	Mean_18  	Mean_19  	Mean_20  	Mean_21  	Mean_22  	Mean_23  	Mean_24  	Mean_25  	Mean_26  	Mean_27  	Mean_28  	Mean_29  	Mean_30  	Mean_31  	Mean_32  	Mean_33  	Mean_34  	Mean_35  	Mean_36  	Mean_37  	Mean_38  	Mean_39  	Mean_40  	Mean_41  	Mean_42  	Mean_43  	Mean_44  	Mean_45  	Mean_46  	Mean_47  	Mean_48  	Mean_49  	Mean_50  	Mean_51  	Mean_52  	Mean_53  	Mean_54  	Mean_55  	Mean_56  	Mean_57  	Mean_58  	Mean_59  	Mean_60  	Mean_61  	Mean_62  	Mean_63  	Mean_64  	Mean_65  	Mean_66  	Mean_67  	Mean_68  	Mean_69  	Mean_70  	Mean_71  	Mean_72  	Mean_73  	Mean_74  	Mean_75  	Mean_76  	Mean_77  	Mean_78  	Mean_79  	Mean_80  	Mean_81  	Mean_82  	Mean_83  	Mean_84  	Mean_85  	Mean_86  	Mean_87  	Mean_88  	Mean_89  	Mean_90  	Mean_91  	Mean_92  	Mean_93  	Mean_94  	Mean_95  	Mean_96  	Mean_97  	Mean_98  	Mean_99  	Mean_100  	Mean_101  	Mean_102  	Mean_103  	Mean_104  	Mean_105  	Mean_106  	Mean_107  	Mean_108  	Mean_109  	Mean_110  	Mean_111  	Mean_112  	Mean_113  	Mean_114  	Mean_115  	Mean_116  	Mean_117  	Mean_118  	Mean_119  	Mean_120  	Mean_121  	Mean_122  	Mean_123  	Mean_124  	Mean_125  	Mean_126  	Mean_127  	Mean_molecular_layer_HP  	Mean_129  	Mean_130  	Mean_131  	Mean_132  	Mean_133  	Mean_134  	Mean_135  	Mean_136  	Mean_137  	Mean_138  	Mean_139  	Mean_140  	Mean_141  	Mean_142  	Mean_143  	Mean_144  	Mean_145  	Mean_146  	Mean_147  	Mean_148  	Mean_149  	Mean_150  	Mean_151  	Mean_152  	Mean_153  	Mean_154  	Mean_155  	Mean_156  	Mean_157  	Mean_158  	Mean_159  	Mean_160  	Mean_161  	Mean_162  	Mean_163  	Mean_164  	Mean_165  	Mean_166  	Mean_167  	Mean_168  	Mean_169  	Mean_170  	Mean_171  	Mean_172  	Mean_173  	Mean_174  	Mean_175  	Mean_176  	Mean_177  	Mean_178  	Mean_179  	Mean_180  	Mean_181  	Mean_182  	Mean_183  	Mean_184  	Mean_185  	Mean_186  	Mean_187  	Mean_188  	Mean_189  	Mean_190  	Mean_191  	Mean_192  	Mean_193  	Mean_194  	Mean_195  	Mean_196  	Mean_197  	Mean_198  	Mean_199  	Mean_200

Note that there is no Mean_0, which I guess makes sense. Oddly though, look at what should be Mean_128. It says Mean_molecular_layer_HP. I have no idea where that came from because the node values are stored as float:

Number of values stored at each pixel = 1
  -- At sub-brick #0 'rank' datum type is float

As a solution, I am using 3dcalc to just add 1 to all nodes. This gives me the expected result of Mean_1 to Mean_201 and I can just subtract 1 to map it back to the annotation value. I just wanted to bring this up because I can’t imagine that this is the intended behavior of 3dROIstats. I understand that it skips over nodes with a 0 value, but why Mean_molecular_layer_HP?

Is there a reason you’re staying in the surface? Are you okay with using the parcellation in the volume space (which will allow you to look at subcortical structures)?

Sorry for the slow response. Your workaround looks reasonable. Because not all labels are found for all subjects, you might consider using a renumbering script, @SUMA_renumber_FS, instead of 3dRank. Usually that script is now called with @SUMA_Make_Spec_FS. It makes it possible to use the same index across all subjects and to use the various output segmentation masks.

That label for mean_molecular_layer_HP is from FreeSurfer It is listed here as FreeSurfer’s index 214. The reranking must be putting that exactly at an integer value, while the others are slightly off. If the ROI value exists, and there is a label for that index exists exactly, then all the labels would appear in the header line, not just that one region. For ROIs, you should use integer values, maybe with :

3dcalc -a mydset -expr a -datum short -nscale -prefix mydset_int

Also it’s better to display the datasets automatically with an integer color scale for ROIs

3drefit -cmap INT_CMAP mydset_int