local/fwe_tests.R

In UNCDEPENdLab/fmri.pipeline: Analysis of fMRI Data on HPC clusters

library(fmri.pipeline)

## ----- 3DFWHMX

# TESTS

# path <- "/proj/mnhallqlab/studies/MMClock/MR_Proc/10637_20140304/mni_5mm_aroma/sceptic_vchosen_ventropy_dauc_pemax_vtime_preconvolve/FEAT_LVL1_run1.feat"
# res4d_file <- file.path(path, "stats", "res4d.nii.gz")
# fwhmx_mask_file <- file.path(path, "mask.nii.gz")

# test <- afni_3dfwhmx$new(input_file = res4d_file, mask_file = fwhmx_mask_file, average = "geometric", ncpus = 2)
# test$run()
# test$is_complete()
# test$get_acf_params() # average ACF params
# test$get_fwhm_by_volume() # FWHM for each sub-brik (volume)
# test$get_acf_by_radius() # FWHM for each sub-brik (volume)
# test$delete_outputs(prompt = FALSE) # delete files generated by this input

## 3dFWHMx list ======

# # get ACF parameters for 8 first-level residual files for testing
path <- "/proj/mnhallqlab/studies/MMClock/MR_Proc/10637_20140304/mni_5mm_aroma/sceptic_vchosen_ventropy_dauc_pemax_vtime_preconvolve"
res4d_files <- list.files(pattern = "res4d.nii.gz", path = path, full.names = T, recursive = T)
fwhmx_mask_files <- list.files(pattern = "mask.nii.gz", path = path, full.names = T, recursive = T)
test <- afni_3dfwhmx_list$new(input_files = res4d_files, mask_files = fwhmx_mask_files, scheduler = "slurm")

# # submit 3dFWHMx jobs to cluster -- the object will divide the list into smaller jobs if needed
# # if 3dFWHMx has completed for all inputs already, the $submit method will just return immediately and note
# # that there is nothing to submit. If you want to force re-estimation, use $submit(force = TRUE).
# test$submit()

# # has 3dFWHMx completed for all datasets provided to the list?
# test$is_complete()

# # average a, b, c parameters across datasets -- good for input to 3dClustSim -acf
# test$get_acf_average()

# # average effective FWHM estimate using -ACF
# test$get_effective_fwhm()

# # per-dataset ACF summary if you want to see what got averaged
# test$get_acf_df()


## 3DCLUSTERIZE
# straight z-stat clusterizing at arbitrary threshold
cobj <- afni_3dclusterize$new(
  threshold_file = "/proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-pe/L2m-l2_l2c-emotion.happy/L3m-age_sex/FEAT_l1c-EV_pe.gfeat/cope1.feat/stats/zstat6.nii.gz",
  lower_thresh = -3, upper_thresh = 3, bisided = TRUE, NN = 1, clust_nvox = 35, pref_map = "zstat6_clusterized.nii.gz"
)

cobj$generate_subclusters(break_nvox = 250, min_subclust_nvox = 25)
cobj$get_clust_df()

#cobj$run() # execute 3dClusterize (use force = TRUE to force re-estimation)
cobj$generate_subclusters(break_nvox = 200, min_subclust_nvox=40) #, max_subclust_nvox=160)



cobj$get_call() # shows command line for this 3dClusterize operation

cobj$is_complete() # TRUE if 3dClusterize is already done
cobj$has_clusters()
cobj$get_outputs() # vector of output filenames from 3dClusterize
clust_nifti <- cobj$get_cluster_map_nifti() # nifti object with clusters.


vv <- cobj$get_clust_df() # return the clusters as a data.frame (see attributes )
cobj$add_whereami() # add a whereami objet to this clusterize and run whereami with the calculated clusters
ww <- cobj$whereami()$get_whereami_df() # data.frame with anatomical labels, roi_num matches vv
cobj$get_clust_df(include_whereami = FALSE)


#ptfce example
cobj <- afni_3dclusterize$new(
  threshold_file = "/proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-abspe/L2m-l2_l2c-overall/L3m-age_sex/FEAT_l1c-EV_abspe.gfeat/cope1.feat/stats/zstat5_ptfce.nii.gz",
  lower_thresh = -5.1246, upper_thresh = 5.1246, bisided = TRUE, NN = 1, clust_nvox = 1, pref_map = "cluster_mask.nii.gz"
)
cobj$run(force=T)
cobj$get_clust_df()


## WHEREAMI


# where_test <- afni_whereami$new(
#   coord_file = "/proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-pe/L2m-l2_l2c-emotion.happy/L3m-age_sex/FEAT_l1c-EV_pe.gfeat/cope1.feat/stats/clusters.1D",
#   coord_file_columns = 1:3,
#   coord_space = "MNI",
#   coord_orientation = "LPI",
#   atlases = c("CA_ML_18_MNI", "Brainnetome_1.0", "CA_N27_LR"),
#   omask = "/proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-pe/L2m-l2_l2c-emotion.happy/L3m-age_sex/FEAT_l1c-EV_pe.gfeat/cope1.feat/stats/zstat_clusterize.nii.gz"
# )

# where_test$get_call()
# where_test$run()
# where_test$get_whereami_df()



# test from 3dClusterize input objet
# where_test_fromobj <- afni_whereami$new(
#   afni_3dclusterize_obj = x,
#   atlases = c("CA_ML_18_MNI", "Brainnetome_1.0", "CA_N27_LR")
# )

# where_test_fromobj$get_call()
# where_test_fromobj$run()
# where_test_fromobj$get_whereami_df()

# where_test_fromobj$get_omask_call()
# where_test_fromobj$get_outputs()


#  -coord_file /proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-pe/L2m-l2_l2c-emotion.happy/L3m-age_sex/FEAT_l1c-EV_pe.gfeat/cope1.feat/stats/clusters.1D'[1,2,3]'

#  -lpi -space MNI
#  -atlas CA_ML_18_MNI -atlas Brainnetome_1.0 -atlas CA_N27_LR
#  -omask /proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-pe/L2m-l2_l2c-emotion.happy/L3m-age_sex/FEAT_l1c-EV_pe.gfeat/cope1.feat/stats/zstat_clusterize.nii.gz -tab


## 3DCLUSTSIM

# for testing
x <- simulate_null_3dttest$new(
  residuals_file = "/proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-abspe/L2m-l2_l2c-overall/L3m-int_only/FEAT_l1c-EV_abspe.gfeat/cope1.feat/stats/res4d.nii.gz",
  mask_file = "/proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-abspe/L2m-l2_l2c-overall/L3m-int_only/FEAT_l1c-EV_abspe.gfeat/cope1.feat/mask.nii.gz",
  n_permutations = 100000,
  njobs = 32
)
x$get_3dttest_calls(include_complete = TRUE)
x$submit()
x$get_permutation_files()


mytest <- afni_3dclustsim$new(
  insdat_file = x$get_permutation_files()["permutation_file"], insdat_mask_file = x$get_permutation_files()["mask_file"],
  scheduler = "slurm", prefix = "test_sdat", out_dir = "/proj/mnhallqlab/users/michael/fmri.pipeline/local",
  clustsim_mask = "/proj/mnhallqlab/lab_resources/standard/mni_icbm152_nlin_asym_09c/mni_icbm152_t1_tal_nlin_asym_09c_mask_2.3mm.nii", ncpus = 8
)
mytest$submit()


mytest <- afni_3dclustsim$new(
  residuals_file = "/proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-pe/L2m-l2_l2c-emotion.happy/L3m-age_sex/FEAT_l1c-EV_pe.gfeat/cope1.feat/stats/res4d.nii.gz",
  residuals_mask_file = "/proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-pe/L2m-l2_l2c-emotion.happy/L3m-age_sex/FEAT_l1c-EV_pe.gfeat/cope1.feat/mask.nii.gz",
  residuals_njobs = 32,  scheduler = "slurm", prefix = "res4d",
  clustsim_mask = "/proj/mnhallqlab/lab_resources/standard/mni_icbm152_nlin_asym_09c/mni_icbm152_t1_tal_nlin_asym_09c_mask_2.3mm.nii", ncpus = 8
)


cobj <- mytest$apply_clustsim(
  statistic_nifti = "/proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-abspe/L2m-l2_l2c-overall/L3m-int_only/FEAT_l1c-EV_abspe.gfeat/cope1.feat/stats/zstat1.nii.gz"
)

cobj$get_clust_df()

cobj2 <- mytest$apply_clustsim(athr = .05, pthr = .001, NN = 1, 
  statistic_nifti = "/proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-abspe/L2m-l2_l2c-overall/L3m-age_sex/FEAT_l1c-EV_abspe.gfeat/cope1.feat/stats/tstat3.nii.gz",
  output_cluster_mask = TRUE, output_thresholded_image = FALSE
)



cobj2$get_clust_df()

mytest$submit()


# --- PTFCE ---


# for testing
zf <- list.files(
  pattern = "zstat[0-9]+\\.nii\\.gz", path = "/proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-abspe/L2m-l2_l2c-overall/L3m-age_sex/FEAT_l1c-EV_abspe.gfeat/cope1.feat/stats",
  full.names = TRUE
)

x <- ptfce_spec$new(
  # z_files = "/proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-abspe/L2m-l2_l2c-overall/L3m-int_only/FEAT_l1c-EV_abspe.gfeat/cope1.feat/stats/zstat1.nii.gz",
  z_files = zf,
  mask_files = "/proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-abspe/L2m-l2_l2c-overall/L3m-int_only/FEAT_l1c-EV_abspe.gfeat/cope1.feat/mask.nii.gz",
  fsl_smoothest_file = "/proj/mnhallqlab/users/michael/mmclock_pe/mmclock_nov2021/feat_l3/L1m-abspe/L2m-l2_l2c-overall/L3m-int_only/FEAT_l1c-EV_abspe.gfeat/cope1.feat/stats/smoothness",
  fwe_p = c(.05, .01)
)

# y <- ptfce_spec$new(
#   gfeat_dir = c(
#     "/proj/mnhallqlab/studies/MMClock/group_analyses/MMClock_aroma_preconvolve_fse_groupfixed/sceptic-clock-feedback-v_entropy-preconvolve_fse_groupfixed/v_entropy/v_entropy-Intercept-Age.gfeat"
#     "/proj/mnhallqlab/studies/MMClock/group_analyses/MMClock_aroma_preconvolve_fse_groupfixed/sceptic-clock-feedback-v_entropy-preconvolve_fse_groupfixed/v_entropy/v_entropy-Intercept.gfeat"
#   )
#   ,
#   fwe_p = c(.05, .01)
# )
# y$submit()