manuscript_reproducible

Anne-Cecile Lesage, PhD, Oliver Zhou, BA, Jiefei Wang, PhD

2025-09-07

This vignette shows how to reproduce results from the companion software paper of the EZFragility CRAN R package for the first seizure of patients 01 and 26. This allows to check how well the EZFragility package reproduces results from the original neural fragility paper (see Figure 4 and extended seizure 4)

Load the package and data

library(EZFragility)
library(Epoch)
library(ggplot2)
library(ggtext)
library(gsignal)
#> 
#> Attaching package: 'gsignal'
#> The following object is masked from 'package:Epoch':
#> 
#>     resample
#> The following objects are masked from 'package:stats':
#> 
#>     filter, gaussian, poly

Data preprocess (Not included in the package)

To allow fast and easy exploration of fragility results, we preprocessed the multipatient data from the OpenNeuro Fragility Data Set with RAVE 2.0 . Main steps included Notch filtering to remove power line interference, re-referencing, and epoching based on seizure onset time to isolate relevant segments of [-30:30s].

Download Data

To access the preprocessed data, you can use the EpochDownloader class from the Epoch package. The following code downloads the data and lists the available datasets.

dl<-EpochDownloader(progress = FALSE)
names(dl)
#>  [1] "FragilityData_subpt2_3"   "FragilityData_subpt13_2" 
#>  [3] "FragilityData_subpt15_1"  "FragilityData_subjh105_2"
#>  [5] "FragilityData_subpt12_1"  "FragilityData_subjh103_1"
#>  [7] "FragilityData_subpt7_2"   "FragilityData_subpt11_2" 
#>  [9] "FragilityData_subpt13_1"  "FragilityData_subpt8_2"  
#> [11] "FragilityData_subpt2_2"   "FragilityData_subjh103_2"
#> [13] "FragilityData_subjh105_1" "FragilityData_subpt6_3"  
#> [15] "FragilityData_subpt6_2"   "FragilityData_subpt8_1"  
#> [17] "FragilityData_subjh101_2" "FragilityData_subjh103_3"
#> [19] "FragilityData_subpt16_1"  "FragilityData_subpt01_1" 
#> [21] "FragilityData_subpt01_4"  "FragilityData_subpt12_2" 
#> [23] "FragilityData_subpt3_1"   "FragilityData_subpt10_1" 
#> [25] "FragilityData_subpt11_4"  "FragilityData_subpt2_1"  
#> [27] "FragilityData_subpt11_3"  "FragilityData_subpt11_1" 
#> [29] "FragilityData_subpt14_2"  "FragilityData_subpt15_2" 
#> [31] "FragilityData_subpt10_3"  "FragilityData_subpt3_2"  
#> [33] "FragilityData_subpt01_3"  "FragilityData_subjh101_1"
#> [35] "FragilityData_subjh105_3" "FragilityData_subpt10_2" 
#> [37] "FragilityData_subpt15_4"  "FragilityData_subpt8_3"  
#> [39] "FragilityData_subpt01_2"  "FragilityData_subpt6_1"  
#> [41] "FragilityData_subpt15_3"  "FragilityData_subpt16_3" 
#> [43] "FragilityData_subjh105_4" "FragilityData_subpt14_1" 
#> [45] "FragilityData_subpt16_2"  "FragilityData_subpt7_3"  
#> [47] "FragilityData_subpt17_2"

The preprocessed voltage data from patient pt01 seizure 1 and pt26 seizure 1 can be loaded by

pt01sz1 <- dl$FragilityData_subpt01_1
pt01sz1
#> Epoch Object @ 1000 Hz: 
#>  Time -30       -29.999   -29.998   -29.997   -29.996   -29.995   -29.994   
#>  G1   -340001.1 -337468.1 -360006   -371602.5 -381722.9 -388522.7 -409074   ... 
#>  G2   -258822.7 -260679   -262803.7 -258883.7 -254481.7 -252272.8 -258038.6 ... 
#>  G3   -239279.3 -233117.6 -223820.4 -206524.9 -194307.3 -175271.1 -163560.5 ... 
#>  G4   123333.6  116130.2  121380.4  133184.7  156112.5  191592.9  235748.3  ... 
#>  G7   63804.24  59270.13  51457.29  36245.07  17488.84  4907.711  -14826.07 ... 
#>  G8   145233.7  132424.9  126525    115485.8  112664.6  106142.3  95139.39  ... 
#>  G9   -245320.8 -238745.9 -241513.8 -253517.7 -264246.5 -271382   -270364.3 ... 
#>  G10  -495154.8 -480099.8 -473566.2 -468136   -469964.3 -477490.8 -485846.7 ... 
#>  G13  209178.2  211959.4  214446.7  230023.3  251683.6  261378.3  274732.9  ... 
#>  G14  170761.3  165173.8  156977.7  156472.7  160976.6  162013.1  158694.1  ... 
#> ...
#> [82 rows x 60001 cols]
#> rowData: [3 vars] name, soz, resected 
#> metaData: [14 vars] patient, type, timeWindow, seizureTime, totalRuns, run, ... 
#> Use tblData, rowData, colData, metaData to get the data

pt26sz1 <- dl$FragilityData_subjh103_1
pt26sz1
#> Epoch Object @ 1000 Hz: 
#>  Time -30       -29.999  -29.998   -29.997   -29.996   -29.995   -29.994   
#>  ABT1 1201953   1185980  1173963   1165929   1161175   1154228   1145490   ... 
#>  ABT2 420150    415039.8 407836.8  399440    392928.7  384582.9  375171.1  ... 
#>  ABT3 362938.1  361439.3 358149.9  351630.6  346282.1  339042.7  331386.2  ... 
#>  ABT4 -98432.64 -96252.9 -96281.68 -98979.26 -98652.86 -104023.3 -114738.3 ... 
#>  MBT1 378468.9  377128.7 375202.9  374722.7  376949.6  376394.8  377590.7  ... 
#>  MBT2 517111.9  508906.5 499356.3  491326    484962.9  477177.7  468958.4  ... 
#>  MBT3 500315.6  496029.1 489656.4  481150    474300.6  463975.6  450947.4  ... 
#>  MBT4 79912.36  75784.27 69818.25  65409     65214.41  67384.24  67453.96  ... 
#>  PBT3 612871.8  605579.9 596700.3  587817.1  579319.7  568749.6  557960.1  ... 
#>  PBT4 267135.9  261987.7 258530.2  255491.5  256614.3  256436.3  256095.7  ... 
#> ...
#> [87 rows x 60001 cols]
#> rowData: [3 vars] name, soz, resected 
#> metaData: [14 vars] patient, type, timeWindow, seizureTime, totalRuns, run, ... 
#> Use tblData, rowData, colData, metaData to get the data

Remove Artifacts

The following function applies a band-pass filter voltage between 0.5 and Nyquist frequency with fourth-order Butterworth filter to the Epoch objects to remove artifacts.

butterworthFilter <- function(epoch) {
    order <- 4

    sampling_freq <- metaData(epoch)$samplingRate
    nyquist_freq <- sampling_freq / 2

    lowpass <- 0.5
    highpass <- nyquist_freq * 0.99
    normalized_freqs <- c(lowpass, highpass) / nyquist_freq
    filter_type <- "pass"
    butter_filter <- gsignal::butter(
        n = order, 
        w = normalized_freqs, 
        type = filter_type)

    # Apply filter to epoch data
    mat <- tblData(epoch)
    
    # Apply zero-phase filtering (filtfilt) to each row
    filtered_data <- gsignal::filtfilt(
        filt = butter_filter,
        x = t(mat))

    filtered_data <- t(filtered_data)
    tblData(epoch) <- filtered_data
    
    epoch
}

We apply the filter to the epochs and clip the data to the relevant time window of [-10:10s] around seizure onset:

pt01sz1Clipped <- pt01sz1 |>
    crop(start=-10, end=10) |>
    butterworthFilter()
pt01sz1Clipped
#> Epoch Object @ 1000 Hz: 
#>  Time -10       -9.999    -9.998    -9.997    -9.996    -9.995    -9.994    
#>  G1   -97434.26 -89359.77 -88724.66 -93025.75 -98619.75 -97673.33 -102924   ... 
#>  G2   -55989.79 -36619.11 -26922.03 -26327.91 -33203.6  -40211.48 -44676.64 ... 
#>  G3   -15242.09 -15774.45 -20890.98 -14905.51 -13759.81 -9549.735 -4954.656 ... 
#>  G4   -49189.35 -45764.11 -55338.7  -61854.82 -68740.64 -79799.13 -100319.3 ... 
#>  G7   -37633.84 -32203.94 -27773.13 -22065.52 -23003.35 -23840.33 -28567.06 ... 
#>  G8   52171.65  40409.96  26612.52  23772.93  28690.68  33429.71  44366.07  ... 
#>  G9   -18672.4  -32423.55 -48413.74 -65187    -73692.1  -75345.62 -69933.66 ... 
#>  G10  -51215.27 -46577.72 -50841.41 -58716.16 -72175.26 -79183.82 -90255.72 ... 
#>  G13  69088.62  66227.37  58157.27  49200.32  33599.13  22246.74  14703.11  ... 
#>  G14  146541.5  134900.6  125193.6  123374.1  124504.1  131296.2  133883    ... 
#> ...
#> [82 rows x 20001 cols]
#> rowData: [3 vars] name, soz, resected 
#> metaData: [14 vars] patient, type, timeWindow, seizureTime, totalRuns, run, ... 
#> Use tblData, rowData, colData, metaData to get the data

pt26sz1Clipped <- pt26sz1 |>
    crop(start=-10, end=10) |>
    butterworthFilter()
pt26sz1Clipped
#> Epoch Object @ 1000 Hz: 
#>  Time -10       -9.999    -9.998    -9.997    -9.996    -9.995    -9.994    
#>  ABT1 27109.7   24202.82  30718.95  39224.46  43013.2   41856.09  41087.12  ... 
#>  ABT2 -39879.77 -43843.39 -44896.01 -43413.08 -44747.49 -46216.34 -45935.69 ... 
#>  ABT3 -43985.83 -50632.48 -57497.43 -64032.84 -71726.81 -79050.59 -83367.55 ... 
#>  ABT4 -2436.194 -6294.849 -5752.691 -3796.922 -1623.323 -269.9269 1914.486  ... 
#>  MBT1 -19074.69 -24714.02 -28192.73 -28945.06 -29617.81 -31179.98 -34303.31 ... 
#>  MBT2 -12829.63 -15129.68 -15765.97 -16955.52 -17905.59 -17673.86 -18355.75 ... 
#>  MBT3 1276.242  -2763.672 -4415.659 -7185.243 -9707.351 -11835.44 -12031.04 ... 
#>  MBT4 -44319.07 -48397.84 -48219.48 -46570.68 -46580.42 -47798.36 -48300.71 ... 
#>  PBT3 -12234.64 -13255.32 -8687.257 -3274.97  -3107.564 -2905.417 -2566.788 ... 
#>  PBT4 -2841.745 -1082.851 3841.862  8728.639  11724.8   10871.25  7927.77   ... 
#> ...
#> [87 rows x 20001 cols]
#> rowData: [3 vars] name, soz, resected 
#> metaData: [14 vars] patient, type, timeWindow, seizureTime, totalRuns, run, ... 
#> Use tblData, rowData, colData, metaData to get the data

Visualize Voltage Plot

The following function reproduce the voltage plot for a display subset of electrodes in the original Figure 4. The SOZ electrode names are highlighted in red as in the original paper.

visualSOZ <- function(epoch, sozNames) {
    p <- plot(epoch)

    elecColor <- rep("black", nrow(epoch))
    elecColor[rownames(epoch) %in% sozNames] <- 'red'
    elecColor <- rev(elecColor) # match the electrode order in the plot

    p + 
    geom_vline(xintercept = 0, color = "black", linetype = "dashed", linewidth = 1)+ 
    theme(axis.text.y = element_markdown(colour = elecColor))
}

To visualize patient 01 seizure 1:

pt01sozName <- rownames(pt01sz1Clipped)[rowData(pt01sz1Clipped)$soz]
pt01Display <- c(pt01sozName, "MLT1", "MLT2", "MLT3", "MLT4")
pt01sz1Reordered <- pt01sz1Clipped[pt01Display, ]
visualSOZ(pt01sz1Reordered, pt01sozName)

To visualize patient 26 seizure 1, we explicitly exclude the electrodes RTG29, RTG30, RTG31, and RTG32 from the SOZ electrodes, as they are not part of the SOZ in the original paper. The electrode names are reordered to match the original figure:

pt26sozName <- rownames(pt26sz1Clipped)[rowData(pt26sz1Clipped)$soz]
excludedElectrodes <- c("RTG29", "RTG30", "RTG31", "RTG32")
pt26sozName <- pt26sozName[!pt26sozName %in% excludedElectrodes]
pt26Display <- c(
    "ABT1", "ABT2",
    pt26sozName[1:16], 
    excludedElectrodes,
    pt26sozName[17:18])
pt26sz1Reordered <- pt26sz1Clipped[pt26Display, ]
visualSOZ(pt26sz1Reordered, pt26sozName)

Compute the Fragility Matrix

The following code computes the fragility matrix using all electrodes and store the results in the Fragility class object pt01sz1Frag

library(doSNOW)
# compute fragility

cl <- makeCluster(parallel::detectCores(), type = "SOCK")
registerDoSNOW(cl)

windowNum <- 250
step <- 125
pt01sz1Frag <- calcAdjFrag(epoch = pt01sz1Clipped, window = windowNum, step = step, parallel = TRUE, nSearch=100L, progress = FALSE)

pt26sz1Frag <- calcAdjFrag(epoch = pt26sz1Clipped, window = windowNum, step = step, parallel = TRUE, nSearch=100L, progress = FALSE)

# Stop the parallel backend
stopCluster(cl)

Fragility Heatmap

The following function plots the fragility heatmap with the same display options as the previous voltage plot. Looking at both plots allows to check correlation between soz patterns

fragHeatmap <- function(frag, sozNames, ranked=FALSE) {
    startTimes <- frag$startTimes

    indexsz <- which(abs(startTimes)<=0.01)
    elecColor <- rep("black", length(frag$electrodes))
    elecColor[frag$electrodes%in% sozNames] <- 'red'
    elecColor <- rev(elecColor)

    plotFragHeatmap(frag = frag, ranked=ranked) + 
    geom_vline(xintercept = indexsz, color = "black", linetype = "dashed", linewidth = 1) + 
    theme(
        axis.text.y = element_markdown(colour = elecColor)
    )
}

PT01 was a surgical success, with the high-fragility electrodes identified by the algorithm matching those identified by clinical epileptologists.

pt01sz1FragReordered <- pt01sz1Frag[pt01Display]
fragHeatmap(pt01sz1FragReordered, pt01sozName, ranked=FALSE)

Setting the “ranked” variable to TRUE in the plotFragHeatmap allows more contrast and render the heatmap closer to original figure 4.

fragHeatmap(pt01sz1FragReordered, pt01sozName, ranked=TRUE)

PT26 was a surgical failure, and the high-fragility electrodes were not consistent with clinical EEG interpretation

pt26sz1FragReordered <- pt26sz1Frag[pt26Display]
fragHeatmap(pt26sz1FragReordered, pt26sozName, ranked=FALSE)

Use the ranked option to get a more contrasted heatmap:

fragHeatmap(pt26sz1FragReordered, pt26sozName, ranked=TRUE)

Compute mean and standard deviation statistics for SOZ and non SOZ electrode group

The following function allows to reproduce the mean and standard visualization statistics of two electrode groups (soz labeled electrodes and non soz labeled electrodes) as in extended figure 4 from the original paper. This plot shows that the fragility biomarker statistics are respectively significantly/not significantly higher in the soz labeled group correlated with the ground truth success/failure outcome for patient 01 and patient 26.

fragDist <- function(frag, sozNames) {
    timeWindows <- frag$startTimes
    timeIdx <- which(timeWindows >= -5 & timeWindows <= 10)
    frag <- frag[, timeIdx]
    plotFragDistribution(frag = frag, groupIndex = sozNames, bandType="SEM", rollingWindow = 1) + 
    geom_vline(xintercept = 0, color = "black", linetype = "dashed", linewidth = 1)
}

Patient 01 seizure 1 (p 1481 Extended Fig.4)

fragDist(pt01sz1Frag[pt01Display], pt01sozName)

Patient 26 seizure 1 (p 1481 Extended Fig.4)

fragDist(pt26sz1Frag[pt26Display], pt26sozName)