ieegio: IO toolbox for Intracranial Electroencephalography

The goal of ieegio is to provide integrated toolbox for common file formats used in intracranial Electroencephalography (iEEG) and deep-brain stimulation (DBS). Besides reading and writing files, ieegio also maps data between volumes and surfaces, chains coordinate transforms, and compares regions of interest. Currently the package is under active development.

Installation

Install the released version from CRAN:

install.packages("ieegio")

Or the development version from r-universe:

install.packages("ieegio", repos = "https://dipterix.r-universe.dev")

Or from GitHub:

# install.packages("pak")
pak::pak("dipterix/ieegio")

ieegio provides additional sample data. To download them:

# list all the sample data
sample_names <- ieegio::ieegio_sample_data()

lapply(sample_names, ieegio::ieegio_sample_data)

All examples below use these sample files.

Electrophysiology

Read a recording, then pull a channel out by number or by label:

library(ieegio)

edf <- read_edf(ieegio_sample_data("edfPlusD.edf"))

# channel metadata as a data frame
head(edf$get_channel_table())

# by number, or by label
channel <- edf$get_channel(1)
channel <- edf$get_channel("squarewave")

To write a recording, build the channels with as_edf_channel. A channel is either a signal or an annotation table:

signal <- sin(seq(0, 10, by = 0.01))

channels <- list(
  as_edf_channel(signal, channel_num = 1, sample_rate = 200, label = "sine"),
  as_edf_channel(
    data.frame(timestamp = c(0, 5), comments = c("start", "end")),
    channel_num = 2
  )
)

path <- tempfile(fileext = ".edf")
write_edf(channels = channels, con = path)

See the ephys article for more.

Volumes and surfaces

read_volume and read_surface are the universal entry points; the format is inferred from the file name. Volumes carry a header, the data, and a transforms list mapping voxel index to other coordinate systems.

volume <- read_volume(ieegio_sample_data("brain.demosubject.nii.gz"))

# index with `[`, plot anatomical slices
volume[128, , ]
plot(volume, position = c(-50, -10, 15), which = "coronal")

Surface files each hold one kind of data, and merge combines them into a single object:

geometry    <- read_surface(ieegio_sample_data("gifti/icosahedron3d/geometry.gii"))
measurement <- read_surface(ieegio_sample_data("gifti/icosahedron3d/rand.gii"))

merged <- merge(geometry, measurement)
plot(merged, name = list("measurements", 1))

See the imaging article for more.

Surface annotations from AFNI/SUMA

NIML datasets (*.niml.dset) are read by read_surface, which decides between labels and continuous values from the dataset itself. Sample data includes matching std.141 geometry and annotations:

geometry <- read_surface(
  ieegio_sample_data("gifti/std.141.lh.inf_200.gii"))

annotation <- read_surface(
  ieegio_sample_data("niml/std.141.lh.aparc.a2009s.annot.niml.dset"))

merge(geometry, annotation)
#> <ieegio Surface>
#>   Header class: basic_geometry
#>   Geometry :
#>     # of Vertex     : 198812
#>     # of Face index : 397620
#>   Annotations: `node label`
#>     # of labels: 76
#>
#> Contains: `geometry`, `annotations`

io_read_niml is the low-level reader if you need the raw element tree. All NIML storage forms are supported (plain text, binary and base64 in either byte order, optionally gzip compressed).

Color maps

read_colormap handles FreeSurfer lookup tables, .annot color tables, AFNI label tables, FSL color maps, and threeBrain/RAVE JSON:

cmap <- read_colormap(
  ieegio_sample_data("niml/rh.std.141.Glasser_HCP.lbl.niml.dset"))
cmap
#> <ieegio Colormap [discrete]>
#>   Colorspace   : RGB
#>   Color stops  : 181  (keys 0 to 180)
#>   Lookup       : 181 labels

# write as a FreeSurfer lookup table
write_colormap(cmap, tempfile(fileext = ".txt"), format = "fs_lut")

Regions of interest

A region of interest is described first and computed later. as_ieegio_roi records the criteria, resolve_roi_as carries them out and hands back geometry in world (RAS) coordinates, and detect_roi_overlap compares two regions.

atlas <- read_volume(ieegio_sample_data("atlases/YBA/YBA690.nii.gz"))

# describe: which voxels count as the region
roi <- as_ieegio_roi(atlas, threshold_lb = 1, threshold_ub = 5)

# compute: turn the description into geometry
resolve_roi_as(roi, "pointcloud")
#> <ieegio ROI [pointcloud]>
#> ... # of Vertex : 3495

Because both sides are resolved to world coordinates first, regions of different kinds can be compared directly. Here the whole atlas is tested against a facial-nerve tract bundle:

tracts <- read_streamlines(
  ieegio_sample_data("streamlines/CNVII_R.trk"), half_voxel_offset = TRUE)

detect_roi_overlap(as_ieegio_roi(atlas, threshold_lb = 1), tracts, radius = 2)
#> <ieegio ROI overlap: `x` and `y` overlap>
#>   x: volume
#>   y: streamlines (annotated)
#>   Units of `y` overlapping `x`: 53.7%
#>   Early stop: FALSE

The result also carries the annotated tracts back, so each streamline knows whether it hit the region.

Coordinate spaces and transforms

Name a space with new_space, then move a surface into it. Transforms can be read from ANTs or FSL FLIRT files, or chained with as_ieegio_transform.

mni <- new_space("MNI152", orientation = "RAS")

surface <- read_surface(ieegio_sample_data("gifti/icosahedron3d/geometry.gii"))
surface_to_surface(surface, space_from = "scanner", space_to = mni,
                   transform = diag(1, 4))
#> ...
#>     # of transforms : 1
#>       Transform Targets : MNI152

Mapping and editing volumes

# volume mask or atlas label to a smoothed mesh
volume_to_surface(atlas, threshold_lb = 1, threshold_ub = 5)

# down-sample or super-sample
resample_volume(volume, new_dim = c(128, 128, 128))

# burn electrode contacts at RAS positions
burn_volume(volume, ras_position = rbind(c(-50, -10, 15), c(-40, -10, 15)),
            col = c("red", "green"), radius = 2)

burn_curve does the same for a trajectory, such as a depth (sEEG) lead or a DBS shaft, given start and end positions.

Supported formats

The following formats have been (or will be) supported. Please feel free to make a wish-list by posting an issue in this repository.

Electrophysiology (examples)

* - NWB format only has low-level support due to its format nature. Please see examples with help(read_nwb)

Imaging (examples)

** - VTK surfaces can be read from all five extensions, but only written as legacy *.vtk, XML *.vtp, or HDF5 *.vtpb.

*** - read_streamlines infers *.vtk and *.vtp from the file name; call io_read_vtk_streamlines directly for *.pvtp and *.vtpb. write_streamlines also writes *.vtpb and *.h5.

The VTK readers and writers use the Python vtk package.

Color maps

Meta data

The HDF5 back-end is resolved at run time, choosing among readNSx, h5lite, hdf5r, and h5py, so no single one of them is required. LazyH5 and LazyFST read HDF5 and FST files lazily, and convert_fst_to_csv and convert_fst_to_hdf5 convert between the formats.