2.0 Unable to deal with onsets like older versions

[smeisler]

Hello,

I am using the most recent version of BayesfMRI 2.0 (commit f0939d1f5b6f47af0c72a7e93713db5ae575d4ee).

I am working with HCP data and defined the onsets the same way as in the HCP vignette. I get the following error:

SETTING UP DATA 
 .. reading in data for session 1
 .. reading in data for session 2
 MAKING DESIGN MATRICES 
Error in rep(0, nsec * downsample): invalid 'times' argument
Traceback:

1. BayesGLM_cifti(cifti_fname = c(fname1_ts, fname2_ts), surfL_fname = fname_gifti_left, 
 .     surfR_fname = fname_gifti_right, onsets = list(onsets1, onsets2), 
 .     TR = TR, nuisance = list(motion1, motion2), dHRF = 1, DCT = 4, 
 .     session_names = fmri_acqs, resamp_res = NULL, num.threads = parallel::detectCores() - 
 .         2, verbose = TRUE, outfile = "/nese/mit/group/gablab/data/HCP/HCP_style/code/bayesfmri/testFullRes", 
 .     avg_sessions = TRUE)
2. make_HRFs(onsets[[ss]], TR = TR, duration = ntime, deriv = dHRF)

Additionally, when I do not explicitly specify dHRF = 1,

Error in BayesGLM_cifti(cifti_fname = c(fname1_ts, fname2_ts), surfL_fname = fname_gifti_left, : length(dHRF) == 1 & dHRF %in% c(0, 1, 2) are not all TRUE
Traceback:

1. BayesGLM_cifti(cifti_fname = c(fname1_ts, fname2_ts), surfL_fname = fname_gifti_left, 
 .     surfR_fname = fname_gifti_right, onsets = list(onsets1, onsets2), 
 .     TR = TR, nuisance = list(motion1, motion2), DCT = 4, session_names = fmri_acqs, 
 .     resamp_res = NULL, num.threads = parallel::detectCores() - 
 .         2, verbose = TRUE, outfile = "/nese/mit/group/gablab/data/HCP/HCP_style/code/bayesfmri/testFullRes", 
 .     avg_sessions = TRUE)
2. stopifnot(length(dHRF) == 1 & dHRF %in% c(0, 1, 2))

Is there a new format for importing onsets?

Thanks, Steven

[smeisler]

Strangely enough, if I run make_HRFs(onsets2, TR=TR, 405) (where 405 is the number of volumes in that session), I am able to get the HRFs just fine, it is just not running in the main compute script.

[smeisler]

I get a different error when trying just a single session:

SETTING UP DATA 
 MAKING DESIGN MATRICES 
 RUNNING MODELS 

 .. LEFT CORTEX ANALYSIS 
Error in BayesGLM(data = session_data, beta_names = beta_names, mesh = NULL, : I detect 16 task based on the design matrix, but the length of beta_names is 8.  Please fix beta_names.
Traceback:

1. BayesGLM_cifti(cifti_fname = fname1_ts, surfL_fname = fname_gifti_left, 
 .     surfR_fname = fname_gifti_right, onsets = onsets1, TR = TR, 
 .     nuisance = motion1, dHRF = 1, resamp_res = NULL, num.threads = parallel::detectCores() - 
 .         2, verbose = TRUE, outfile = "/nese/mit/group/gablab/data/HCP/HCP_style/code/bayesfmri/testFullRes", 
 .     avg_sessions = TRUE)
2. BayesGLM(data = session_data, beta_names = beta_names, mesh = NULL, 
 .     vertices = surf_list[[each_hem]]$vertices, faces = surf_list[[each_hem]]$faces, 
 .     scale_BOLD = scale_BOLD, scale_design = FALSE, Bayes = do_Bayesian, 
 .     EM = do_EM, ar_order = ar_order, ar_smooth = ar_smooth, aic = aic, 
 .     num.threads = num.threads, return_INLA_result = return_INLA_result, 
 .     outfile = outfile_name, verbose = verbose, avg_sessions = avg_sessions, 
 .     meanTol = meanTol, varTol = varTol, emTol = emTol, trim_INLA = trim_INLA)
3. stop(paste0("I detect ", K, " task based on the design matrix, but the length of beta_names is ", 
 .     length(beta_names), ".  Please fix beta_names."))

For reference, here is what onsets1 looks like:

[damondpham]

Hi @smeisler , the onset file input requirements should not have changed. Could you send over your data for me to take a look?

Also, are you certain that you are using the latest 2.0 branch? I'm confused that you encountered this problem:

Additionally, when I do not explicitly specify dHRF = 1,

Error in BayesGLM_cifti(cifti_fname = c(fname1_ts, fname2_ts), surfL_fname = fname_gifti_left, : length(dHRF) == 1 & dHRF %in% c(0, 1, 2) are not all TRUE
Traceback:

1. BayesGLM_cifti(cifti_fname = c(fname1_ts, fname2_ts), surfL_fname = fname_gifti_left, 
 .     surfR_fname = fname_gifti_right, onsets = list(onsets1, onsets2), 
 .     TR = TR, nuisance = list(motion1, motion2), DCT = 4, session_names = fmri_acqs, 
 .     resamp_res = NULL, num.threads = parallel::detectCores() - 
 .         2, verbose = TRUE, outfile = "/nese/mit/group/gablab/data/HCP/HCP_style/code/bayesfmri/testFullRes", 
 .     avg_sessions = TRUE)
2. stopifnot(length(dHRF) == 1 & dHRF %in% c(0, 1, 2))

It's familiar to me but in an older version--I thought I since fixed it on the latest branch.

[smeisler]

I was jumping around versions yesterday; just reupgraded to most recent and am no longer getting the dHRF error, but still am getting the onset error.

Data are in the following link. The code I use to set up my files is below:

data_dir <- '/PATH/TO/HCP/DIR'
subject <- 100307
task='WM'
TR=0.72
tasks <- c('0bk_body','0bk_tools','0bk_places','0bk_faces',
           '2bk_body','2bk_tools','2bk_places','2bk_faces') # Task event filenames
names_tasks <- c('0-Back Body','0-Back Tools','0-Back Places','0-Back Faces',
                 '2-Back Body','2-Back Tools','2-Back Places','2-Back Faces')
K <- length(tasks)
fmri_acqs <- c('LR','RL')

# Load surfaces
dir_s <- file.path(data_dir, subject, 'T1w', 'fsaverage_LR32k') # Note this is in sub space, not MNI, but I get the error in both spaces
fname_gifti_left <- file.path(dir_s, paste0(subject,'.L.midthickness_MSMAll.32k_fs_LR.surf.gii'))
fname_gifti_right <- file.path(dir_s, paste0(subject,'.R.midthickness_MSMAll.32k_fs_LR.surf.gii'))

# Get fMRI images and event info
for(i in 1:length(fmri_acqs)){
    acq = fmri_acqs[i]
    func_dir <- file.path(data_dir,subject, 'MNINonLinear', 'Results', paste('tfMRI',task,acq,sep="_"))
    assign(paste0('fname',i,'_ts'), file.path(func_dir,paste('tfMRI',task,acq,'Atlas_MSMAll.dtseries.nii',sep="_")))
    assign(paste0('motion',i),as.matrix(read.table(file.path(func_dir,'Movement_Regressors.txt'), header=FALSE)))
    onsets <- vector('list', length=K)
    names(onsets) <- names_tasks
    for(t in tasks){
        ind_t <- which(tasks==t)
        fname_t <- file.path(func_dir,'EVs',paste0(t,'.txt'))
        onsets[[ind_t]] <- read.table(fname_t, header=FALSE)
    }
    assign(paste0('onsets',i), onsets)
}

# Run the GLM
thetas <- NULL # No starting values for precision parameters
results <- BayesGLM_cifti(
    cifti_fname = c(fname1_ts, fname2_ts),
    surfL_fname = fname_gifti_left,
    surfR_fname = fname_gifti_right,
    onsets = list(onsets1, onsets2),
    TR = TR,
    nuisance = list(motion1, motion2),
    session_names = fmri_acqs,
    resamp_res = NULL,
    num.threads =  parallel::detectCores() - 2,
    verbose = TRUE,
    outfile = '/PATH/TO/OUTFILE',
    avg_sessions = TRUE)

[damondpham]

Thanks @smeisler ! One issue I think I see is that for multisession input, onsets should be a session-length list of task-length lists, something like this: onsets <- list(sess1=list(ta1, tb1, tc1), sess2=list(ta2, tb2, tc2). In the code above, it looks like onsets ends up being just list(ta2, tb2, tc2).

I found an issue in the code for multisession input, and I just pushed a fix on 2.0.

Unfortunately, there's still the CIFTI mask-mismatch error that you saw. I will be looking into that error later today.

[smeisler]

Thanks for your quick investigation and fixes, really appreciate it!

[damondpham]

Actually, I see there is onsets which looks like list(ta1, tb1, tc1), but then in the call to BayesfMRI there is onsets = list(onsets1, onsets2) which looks potentially correct, though I can't say for sure since I don't know what onsets1 and onsets2 are.

[smeisler]

Yes, sorry for the confusion, I wanted to put everything in a loop that loops over the "RL" and "LR" acquisitions.

Onsets are loaded the same way as in the HCP vignette, but are assigned a variable name within the loop with the line:

assign(paste0('onsets',i), onsets), so it is just naming onset1, onset2, onsetX dynamically within the loop. onset1 and onset2 should match the same format as the vignette, as noted in one of my comments above.

These onsets did work for a previous version of BayesfMRI.

[damondpham]

Ohh right I forgot about assign here! My bad.

I think I pushed a fix to ciftiTools 11.0 which you can get from GitHub. Could you try again with the patches to ciftiTools@11.0 and BayesfMRI@2.0?

[smeisler]

When I do not do any resampling, I still get the mask mismatch error, but when resampling to 10000 points it seems to proceed (at least no immediate error yet). No onset/design matrix errors yet, which is encouraging, especially since I am using the multisession command.

[damondpham]

Okay cool! I thought the recent fix to ciftiTools should fix the mask mismatch error. Are you sure you are using the recently updated 2.0? If you do ciftiTools:::make_xifti to see the source code and look for where HCP_32k_auto_mwall appears, does it say something like

    if (HCP_32k_auto_mwall & nrow(xifti$data$cortex_left) == 29696) {

or

    if ((HCP_32k_auto_mwall & is.null(cortexL_mwall)) & nrow(xifti$data$cortex_left) == 29696) {

? It should look like the latter, not the former.

[smeisler]

The second: if ((HCP_32k_auto_mwall & is.null(cortexR_mwall)) & nrow(xifti$data$cortex_right) == 29716) and

if ((HCP_32k_auto_mwall & is.null(cortexL_mwall)) & nrow(xifti$data$cortex_left) == 
            29696)

[damondpham]

oh okay, I forgot I needed to add a patch to BayesfMRI 2.0 as well, to go along with the new change! I just added it. could you try again after re-installing BayesfMRI 2.0 ?

[smeisler]

So far no errors on the NULL resampling! Will update later if it finishes or not (might be a while)

[damondpham]

Okay great, let me know! Thanks as usual for your patience helping me sort out these problems!

[smeisler]

Looks like it works (job ended up crashing later but it was due to memory restraints), so I'll close my two issues now!