| SOCR ≫ | TCIU Website ≫ | TCIU GitHub ≫ |
Functions & Workflow of TCIU Analytics
Validation using 4D Complex (mag, phase) and Real (solely-mag reconstruction) fMRI data
¶ SOCR Team ¶
09/14/2024
1 Background
This document is set up for package TCIU which is developed in the R environment.
In this vignette, we introduce the detailed usage of all the functions with examples in the TCIU R package. Besides, we focus in the use of fmri_stimulus_detect() function on our simulated brain fMRI data (generated by fmri_simulate_func()) to get the p-value to locate the stimulated parts in the brain as well as the post-hoc process and result visualization. User can also refer to the concise steps in our website TCIU Predictive Analytics.
require(TCIU)
require(DT)2 Function list
2.1 FMRI data simulation
- fmri_simulate_func: a real-valued fMRI data simulation function, used to simply generate a 3D fMRI data associated with brain area with activated parts inside.
2.2 Visualization of the fMRI data in time-series
fmri_time_series: create four interactive time series graphs for the real, imaginary, magnitude, and phase parts for the fMRI spacetime data.
fmri_kimesurface: transform the fMRI time-series data at a fixed voxel location into a kimesurface.
fmri_image: create images for front view, side view, and top view of the fMRI image.
fmri_ts_forecast: forecast the fMRI data based on the time series
2.3 Activated areas detection
fmri_stimulus_detect: take a real/complex valued fMRI data and detects locations where stimulus occurs
fmri_post_hoc: conduct the post-hoc process (i.e. FDR correction and spatial clustering) for a 3-dimensional p-value array.
2.4 Activated areas visualization
fmri_2dvisual: a visualization method, use ggplot2 to draw the brain from axial, sagittal and coronal view with activated area identified by p-values.
fmri_pval_comparison_2d: a plots arrangement method, which uses gridExtra to combine multiple 2d plots of the fMRI data together. This can bring convenience for users to compare the performance of different statistical tests based on the p-values they provide.
fmri_3dvisual: a visualization method, use plotly to draw the 3D plot of the brain with the activated areas determined by p-values.
fmri_pval_comparison_3d: a visualization method, use plotly to compare the activated parts inside in the brain using two sets of color palette. The activated parts are decided by p-values.
2.5 Activated areas detection by three-phase analysis
fmri_3dvisual_label: a visualization method to display 3d p-values region by region
fmri_ROI_phase1: calculate p-values on region of interest(ROI) of the brain
fmri_ROI_phase2: calculate tensor-on-tensor regression on region of interest(ROI) of the brain
3 Analysis of simulated fMRI data
Here we simulate our 4D fMRI data with its first three dimension as 64 times 64 times 40 and its fourth dimension as the 160-time length. We specify the dimension of the data, starting time points when the fMRI data receives the stimulus and the duration of each stimulated time periods to get it. As we want our data shaped like a brain, we provide a brain-shaped mask for it. Otherwise, it will generate a sphere mask automatically.
Details/Parameters
The function fmri_simulate_func() is used to simulate fMRI data with specified dimension and total time points. The fMRI data can be brain-shaped by using the mask data provided in our package, if the dimension fits the same as our data (c(64, 64. 40)). Otherwise, the function will generate a 3D sphere data with multiple activated part inside. The activated parts can be detected based on the p values.
Parameters in this function:
dim_data : a vector of length 3 to identify the dimension of fMRI data user wants to simulate
mask : a 3D array of 1’s and 0’s or NULL. To specify the area inside the brain shell. One may use the mask data provided by this package, or generate a 3D array of 1’s and 0’s of the same dimension with the fMRI data to be generated. If NULL, then the function would generate a 3D sphere mask.
ons : a vector of the start time points of the time period when the fMRI data receives stimulation
dur : a vector of the time period when the fMRI data receives stimulation.
fmri_generate = fmri_simulate_func(dim_data = c(64, 64, 40), mask = mask,
ons = c(1, 21, 41, 61, 81, 101, 121, 141),
dur = c(10, 10, 10, 10, 10, 10, 10, 10))
# the outputs include simulated fMRI data, its mask,
# the starting time points of the stimulated period and its duration
# as well as all the stimulated time points
dim(fmri_generate$fmri_data)## [1] 64 64 40 1604 Visualization of the fMRI data in time-series
4.1 Function “fmri_time_series”
Details/Parameters
We first provide an interactive time-series visualization for a fixed voxel. Since the input fMRI data needs to be complex-valued, we provide a sample complex voxel’s time-series vector and a sample reference plotly object.
The function fmri_time_series() is used to create four interactive time series graphs for the real, imaginary, magnitude, and phase parts for the fMRI spacetime data.
Parameters in this function:
fmridata : a 4D array which contains the spatial and temporal record of fMRI result or a single complex valued vector
voxel_location : a 3D array indicating the spatial location of the brain. If is.4d is false, set the voxel_location as NULL.
is.4d : The default is TRUE. If change to false, input a vector instead of a 4D array.
ref : The default is NULL. User can input an outside extra reference object to include in the final result.
Example
fmri_time_series(sample_save[[9]], voxel_location = NULL, is.4d = FALSE, ref = sample_save[[8]])4.2 Function “fmri_kimesurface”
Details/Parameters
We second provide a function to transform the fMRI time-series at a fixed voxel location into a kimesurface (kime-series).
The function fmri_kimesurface() is used to display in 3D the kime-series as 2D manifolds (kimesurface) over the Cartesian domain. User can choose to provide the 4D array of the fMRI spacetime image and the voxel_location or a single TS vector.
Parameters in this function:
fmridata : a 4D array which contains the spatial and temporal record of fMRI result or a single real valued vector.
voxel_location : a 3D array indicating the spatial location of the brain.
is.4d : The default is true. If change to false, need to input a vector instead of array.
Example
# a data-frame with 160 rows and 4 columns: time (1:10), phases (8), states (2), and fMRI data (Complex or Real intensity)
datatable(fmri_kimesurface(fmri_generate$fmri_data, c(44,30,33))[[1]])# ON Kime-Surface
fmri_kimesurface(fmri_generate$fmri_data, c(44,30,33))[[2]]# User can try themself to plot the on / off / on&off figure
# OFF Kime-Surface
# fmri_kimesurface(fmri_generate$fmri_data, c(44,30,33))[[3]]
# ON&OFF Kime-Surface
# fmri_kimesurface(fmri_generate$fmri_data, c(44,30,33))[[4]]4.3 Function “fmri_image”
Details/Parameters
The function fmri_image() is used to create images for front view, side view, and top view of the fMRI image. When providing the 4D array of the fMRI spacetime image and input the x,y,z position of the voxel, three views of the fMRI image and the time series image of the voxel will be shown.
Parameters in this function:
fmridata : a 4D array contains information for the fMRI spacetime image. The data should only contain the magnitude for the fMRI image.
option : The default is ‘manually’. If choose ‘auto’, then this function will lead you to key in the space (x,y,z) parameters and time (time) parameter for this function to generate graphs.
voxel_location : a 3D array indicating the spatial location of the brain. If option is auto, set the voxel_location as NULL.
time : time location for the voxel
Example
fmri_image(fmri_generate$fmri_data, option="manually", voxel_location = c(40,22,33), time=4)4.4 Function “fmri_ts_forecast”
Details/Parameters
The function fmri_ts_forecast() is used to forecast the fMRI data based on the time series.
Parameters in this function:
fmridata : a 4D array contains information for the fMRI spacetime image. The data should only contain the magnitude for the fMRI image.
voxel_location : a 3D array indicating the voxel location of the brain
cut : breaking point of the time-series data. The default is 10.
Example
smoothmod<-GaussSmoothArray(fmri_generate$fmri_data, sigma = diag(3,3))
fmri_ts_forecast(smoothmod, voxel_location=c(41,44,33))5 Activated areas detection
5.1 Load and display the data
Here we apply fmri_stimulus_detect() function on our simulated fMRI data to get the 3D array p value data. We use t-test here for simplicity. But we also have Wilcoxon test, HotellingT2 test, gLRT test for complex data, generalized likelihood test, generalized linear model and on-off intensity difference method to generate the p value to detect the stimulated parts in the fMRI data. The smaller the p value is, the higher level the stimulated part is.
p_simulate_t_test =
fmri_stimulus_detect(fmridata= fmri_generate$fmri_data,
mask = fmri_generate$mask,
stimulus_idx = fmri_generate$on_time,
method = "t-test" ,
ons = fmri_generate$ons,
dur = fmri_generate$dur)
dim(p_simulate_t_test)## [1] 64 64 40summary(p_simulate_t_test)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.9125 1.0000 0.8625 1.0000 1.0000From the dimension of the generated p_simulate_t_test, we know it is a 3D array p value with its dimension as 64 times 64 times 40.
5.2 Analysis of real fMRI data
As we save a 3D array p value data generated from the real fMRI brain data, here we can apply post-hoc test on it to improve the performance. We can apply FDR correction and spatial clustering to our 3D array p value data to prevent the high proportion of false positive voxels and keep only spatially connected significant voxels. Notice that when doing spatial clustering, the threshold(spatial_cluster.thr) and cluster size(spatial_cluster.size) are needed to filter contiguous clusters of locations in a 3D array that are below some threshold and with some minimum size.
# do the FDR correction
pval_fdr = fmri_post_hoc(phase2_pval , fdr_corr = "fdr",
spatial_cluster.thr = NULL,
spatial_cluster.size = NULL,
show_comparison = FALSE)
# do the spatial clustering
pval_posthoc = fmri_post_hoc(pval_fdr, fdr_corr = NULL,
spatial_cluster.thr = 0.05,
spatial_cluster.size = 5,
show_comparison = FALSE)5.3 Visualization and Comparison of fMRI data
5.3.1 2D and 3D visualization
We generate the 2D plot and 3D interactive of our simulated fMRI data based on the p value generated above from sagittal, coronal and axial view.
# the output figure is hidden
for(axis in c("x", "y", "z")){
axis_i = switch(axis,
"x" = {35},
"y" = {30},
"z" = {22})
print(fmri_2dvisual(p_simulate_t_test, list(axis, axis_i),
hemody_data=NULL, mask=fmri_generate$mask,
p_threshold = 0.05, legend_show = TRUE,
method = "scale_p",
color_pal = "YlOrRd", multi_pranges=TRUE))
}fmri_3dvisual(p_simulate_t_test, fmri_generate$mask,
p_threshold = 0.05, method="scale_p",
multi_pranges=TRUE)$plot5.3.2 Comparison of performance of different methods on the same fMRI data
We can use fmri_pval_comparison_3d() to visualize two p value data to see how they perform differently on detecting stimulated parts by 3D plot. Here we compare the difference of stimulated parts of two different fMRI data with the same mask, since our original fMRI is too big to use here for different statistical test. The output figure is similar to 3D visualization, so we hide them here.
# the two p value are the p value generated based on the simulated fMRI
# and the p value saved in the package and finished post hoc test
# the output figure is hidden
fmri_pval_comparison_3d(list(p_simulate_t_test, phase3_pval), mask,
list(0.05, 0.05), list("scale_p", "scale_p"),
multi_pranges=FALSE)We can also use fmri_pval_comparison_2d() to visualize whatever number of p value (generated by different statistical tests on the same fMRI data) to see their difference by 2D plot. For simplicity here we only compare two different 3D array p value data.
fmri_pval_comparison_2d(list(p_simulate_t_test, phase3_pval),
list('pval_simulated', 'pval_posthoc'),
list(list(35, 33, 22), list(40, 26, 33)),
hemody_data = NULL,
mask = mask, p_threshold = 0.05,
legend_show = FALSE, method = 'scale_p',
color_pal = "YlOrRd", multi_pranges=FALSE)
6 Activated areas detection by three-phase ROI Analysis
6.1 Phase1: Detect Potential Activated ROI
Details/Parameters
The function fmri_ROI_phase1 is used to calculate p-values of ROIs for a given real-valued fMRI data, which is usually the 1st phase of a ROI 3-phase analysis.
fmridata : a 4D array contains information for the fMRI spacetime image. The data should only contain the magnitude for the fMRI image.
label_mask : a 3D nifti or 3D array of data to indicates the corresponding indices of the ROIs
label_dict : a dataframe which contains the name of ROIs and their corresponding index
stimulus_idx : a vector that specifies when motion happens
rest_idx : a vector that specifies when study participant does not move
Example
ROI_phase1 = fmri_ROI_phase1(fmri_generate$fmri_data, mask_label, mask_dict, stimulus_idx = fmri_generate$on_time)6.2 Phase2: ROI-Based Tensor-on-Tensor Regression
Details/Parameters
The function fmri_ROI_tensor_regress is used to detect locations where stimulus occurs by calculating the p-values of the ROI-based tensor-on-tensor regression. Two methods, namely ‘t_test’ and ‘corrected_t_test’,can be chosen to calculate the p-values from the regression coefficients.
fmridata : a 4d array which contains the spatial and temporal record of fmri result.
label_mask : a 3d nifti or 3d array of data that shows the labeled brain atlas.
ROI_label_dict a dataframe or array or matrix to specify the indices and corresponding names of the ROI. The input of this parameter could take one of the list outputs of the fmri_ROI_detect function as a following step.
stimulus_idx : a vector of the start time points of the time period when the fMRI data receives stimulation.
stimulus_dur : a vector of the time period when the fMRI data receives stimulation.
fmri.design_order : a parameter to specify the order of the polynomial drift terms in the fmri.design function.
fmri.stimulus_TR : a parameter to specify the time between scans in seconds in the fmri.stimulus function.
rrr_rank : a parameter to specify the assumed rank of the coefficient array in the rrr function.
method : a string that represents method for calculating p-values from tensor-on-tensor regression coefficients. There are 2 options: ‘t_test’ and ‘corrected_t_test’. The default is ‘t_test’. ‘t_test’ is to calculate the test statistics ‘t-value’ across all voxels in the bounding box of ROI; ‘corrected_t_test’ is to calculate the test statistics ‘t-value’ by first across each voxel on a temporal basis, and then across all voxels in the bounding box of ROI.
parallel_computing : a logical parameter to determine whether to use parallel computing to speed up the function or not. The default is FALSE.
ncor : number of cores for parallel computing. The default is the number of cores of the computer minus 2.
Example
ROI_phase2 = fmri_ROI_phase2(fmridata = fmri_generate$fmridata, label_mask = mask_label,
ROI_label_dict = mask_dict, stimulus_idx = fmri_generate$on_time,
stimulus_dur = fmri_generate$dur, rrr_rank = 3,
fmri.design_order = 2, fmri.stimulus_TR = 3,
method = "t_test", parallel_computing = TRUE, max(detectCores()-2,1))6.3 Phase3: FDR Correction and Spatial Clustering
Example
# do the FDR correction
# do the spatial clustering
ROI_phase3 = fmri_post_hoc(ROI_phase2 , fdr_corr = "fdr",
spatial_cluster.thr = 0.05,
spatial_cluster.size = 5,
show_comparison = FALSE)6.4 3D visualization based on the activated areas by regions
Details/Parameters
The function fmri_3dvisual_region is used to visualize the 3D plot of the brain with activated parts region by region. The parameters are very similar to fmri_3dvisual with some additional parameters to assign the region information.
Example
# the output figure is hidden due to the size of vignettes
label_index = mask_dict$index
label_name = as.character(mask_dict$name)
label_mask = mask_label
fmri_3dvisual_region(phase1_pval, mask_label, label_index, label_name, title = "phase1 p-values")# the output figure is hidden due to the size of vignettes
fmri_3dvisual_region(list(phase2_pval,phase3_pval), mask_label,
label_index, label_name, title = "phase2&3 p-values")