Reading FreeSurfer stats files
In ggseg: Plotting Tool for Brain Atlases

#| label: setup
#| include: false
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = FALSE
)

FreeSurfer's recon-all produces stats files with region-level measures -- cortical thickness, surface area, volume. ggseg includes readers that pull these files directly into data frames shaped for plotting.

All three functions below come from the ggseg.formats package, which ggseg re-exports.

#| label: load-packages
library(ggseg)
library(ggplot2)
library(dplyr)

Reading a single stats file

After recon-all finishes, each subject has a stats/ folder with parcellation data. read_freesurfer_stats() reads one file at a time:

#| label: read-single-stats
subjects_dir <- Sys.getenv("SUBJECTS_DIR")
stats_file <- file.path(subjects_dir, "bert/stats/lh.aparc.stats")

data <- read_freesurfer_stats(stats_file)
data

The hemisphere is encoded in the filename (lh. or rh.). To match the atlas, prepend it to the label column:

#| label: fig-single-stats
#| fig-cap: "Brain plot of average cortical thickness from a single stats file."
ggplot(data |> mutate(label = paste0("lh_", label))) +
  geom_brain(atlas = dk(), mapping = aes(fill = ThickAvg))

Reading stats files across subjects

read_atlas_files() reads all matching stats files from a subjects directory. Pass a regular expression to select the right files:

#| label: read-atlas-files
dat <- read_atlas_files(subjects_dir, "aparc.stats$")
dat

The trailing $ matters -- it ensures you get aparc.stats but not aparc.a2009s.stats.

The function adds hemisphere prefixes automatically, so the data is ready for plotting:

#| label: fig-atlas-files
#| fig-cap: "Brain plot of cortical thickness variability across subjects."
ggplot(dat) +
  geom_brain(atlas = dk(), mapping = aes(fill = ThickStd))

To show all metrics at once, pivot and facet:

#| label: fig-all-metrics
#| fig-cap: "All FreeSurfer metrics shown in faceted brain plots."
#| fig-width: 10
library(tidyr)

dat |>
  pivot_longer(-c(subject, label), names_to = "stat", values_to = "val") |>
  ggplot() +
  geom_brain(atlas = dk(), mapping = aes(fill = val)) +
  facet_wrap(~stat)

Reading FreeSurfer stats tables

FreeSurfer's aparcstats2table and asegstats2table commands create summary tables across subjects. Read these with read_freesurfer_table():

#| label: read-fs-table
table_path <- "path/to/aparc.volume.table"

read_freesurfer_table(table_path)

The output has three columns: subject, label, and value.

Cleaning label suffixes

Stats tables append the measure name to each label (e.g. lh_bankssts_volume). The measure argument strips that suffix and renames the value column:

#| label: read-fs-table-measure
dat <- read_freesurfer_table(table_path, measure = "volume")
dat

Filtering non-region rows

FreeSurfer tables include summary measures (total volume, ICV) that don't correspond to atlas regions. Filter them before plotting:

#| label: fig-filtered-table
#| fig-cap: "Brain plot of volume data after filtering non-region rows."
ggplot(dat |> filter(grepl("lh|rh", label))) +
  geom_brain(atlas = dk(), mapping = aes(fill = volume))