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)
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].
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
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
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)
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)
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.
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:
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)
Patient 26 seizure 1 (p 1481 Extended Fig.4)