-
Notifications
You must be signed in to change notification settings - Fork 5
resources afni
The sub-package afni contains a number of modules:
- copy
- deconvolve
- group
- masks
- motion
- process
- submit
These modules are organized, largely, according to the type of work performed. That is, any functions that deal with mask generation are located in the module masks. Together, these modules are capable of conducting resting state and task-based EPI analyses.
copy contains a single function, copy_data. It is the role of this module to (a) copy structural, functional, and motion files needed by AFNI, and (b) to start and populate the afni_data dictionary that will be used to pass files to various modules.
The deconvolve module currently contains five functions used for preparing and conducting deconvolutions of the EPI time series:
- write_decon
- write_new_decon
- run_reml
- timing_files
- regress_resting
write_decon constructs a 3dDeconvolve command that incorporates censor files, motion (mean and derivative) baseline regressors, and task timing files for the various pre-processed (scaled) EPI runs.
Currently, only a two gamma basis function is supported, AFNI is allowed to determine the number of polynomial regressors (using local times), and output names follow traditional AFNI formatting. The AFNI naming convention is due to the fact that 3dDeconvolve is employed to generate (a) the deconvolution matrices and (b) the 3dREMLfit command.
The generated 3dDeconvolve command is saved to a script in the subject's func directory, for review, and then the command is executed in order to generate the files listed above.
The write_new_decon function is similar to the write_decon above, save for how motion censoring is handled. Rather than removing volumes from the analysis, the behavioral timing files are converted into binary vectors (1 = behavior, 0 = not) of length equal to the number of volumes. The binary motion censor vector (1 = good data, 0 = motion) is then used to censor any behavior that co-occurs with a motion event via matrix multiplication. The resulting censored behavior vector is then convolved with a two-gamma basis function and incorporated in the deconvolution via the -stim_file option. An inverted binary motion censor file (0 = good data, 1 = motion) is also included as a baseline regressor, prior to the behavior regressors.
run_reml has two parts. First, a white matter time series is constructed from an eroded white matter mask. Then the foo.REML_cmd generated by write_decon is executed with the addition of the white matter time series as a nuissance regressor.
The timing_files function is specific to the EMU project, but could easily be updated for other projects. When user-specified timing files/deconvolution are supplied, the timing_files module is not used. timing_files will orient to the foo_events.tsv BIDS files that accompany each EPI run in func. From these events files all unique behaviors will be identified (ref. the "trial_type" column). These unique behaviors will then be used to construct AFNI-styled timing files using onset times only and the asterisk "*" holding for runs in which a certain behavior does not occur. This could be updated to construct onset:duration, but is not supported currently.
The decon_plan nested dictionary is constructed by the timing_files module, where the parent key is the deconvolution name, child key is the behavior (<8 characters!), and value is the path to the timing file. The behavior name is embedded in the file name as a description. For example:
{
"UniqueBehs": {
"behA": "/path/to/foo_desc-behA_events.1D",
"behB": "/path/to/foo_desc-behB_events.1D",
"behC": "/path/to/foo_desc-behC_events.1D",
}
}
With the decon_plan above, an output file "decon_task-foo_UniqueBehs_stats_REML+tlrc" will be generated by modeling three behaviors, behA-C, and "behA" etc will be used to label the coefficient and statistics sub-bricks of the deconvolved file.
The regress_resting function is used to project a regression matrix for resting state analyses. It is currently written for only a single pre-processed (scaled) EPI file, but this could be updated. There are three main steps.
First, a principal components analysis is conducted on an eroded CSF mask in order to identify a nuissance time series. Second, a deconvolution via 3dDeconvolve is conducted. Here, only the censor file and nuissance regressors (CSF timeseries, mean motion, derivative motion) are included in the model, and so 3dDeconvolve is actually being used to "clean" the data, producing a set of nuissance models on the main output file (foo_stats+tlrc) and the remaining "clean" data as a residual (foo_errts+tlrc). As in write_decon, the deconvolution command is written as a shell script to the participant's func directory for review.
Third, the "no censor" matrix produced by 3dDeconvolve is used to project a correlation matrix using the pre-processed (scaled) EPI file as input. By default, the function will use the "anaticor" method, where the time series of an eroded white matter mask will be included in the projection as a nuissance regressor.
The group module currently contains three functions used for preparing and conducting group-level analyses:
- int_mask
- resting_etac
- task_etac
The int_mask function constructs a group-level gray matter intersection mask. That is, each participants' EPI-anat intersection mask is combined (union/frac=1) to construct a group intersection mask. This mask indicates which voxels have meaningful signal across all participants. The group intersection mask is then multiplied by a template gray matter mask. This helps focus analyses to cortical and sub-cortical regions (excluding white matter, CSF) and also serves to reduce the multiple comparison correction as fewer voxels are being included in models of spurious results.
The resting_etac function constructs an "A vs not-A" comparison, or a test against zero. Z-transformed correlation matrices are used as input, and analyses are constrained to gray matter voxels via the group-level gray matter intersection mask (above). With ETAC (equitable thresholding and clustering), a variety of p-value threshold are used in order to detect both large clusters of a "smaller" statistic and smaller clusters with a "larger" statistic. Here, default options are NN=2, two-sided, p-values = .01, .005, .002, .0001.
Additionally, a number of blurs/smooths are supported by ETAC. It is recommended the blurring be conducted by ETAC if it is planned to blur your data. Here, two blur sizes are hardcoded, 4 and 8, and these are selected used if data was not blurred as part of pre-processing. A script is written to the analyses directory for review.
task_etac first determines which participants have both of the supplied behaviors, and determines which sub-brick of the deconvolved file was assigned the behavior coefficient. Next, and "A vs B" paired 3dttest++ command is constructed, using the group gray matter intersection mask, noise simulations, and the same ETAC options as resting_etac. The command is written as a shell script to the output directory for researcher review, and finally the command is executed is the output is not detected.
The masks module currently contains three functions used for construct anatomic and functional masks:
- make_intersection_mask
- make_tissue_masks
- make_minimum_masks
make_intersection_mask constructs a binary mask of voxels where (a) sufficient fMRI signal is detected and (b) where an anatomic mask occurs. First, 3dAutomask is used to detect regions of >minimum EPI signal for each run, and default options are supplied for future utility. Next, these run EPI masks are used to construct a union mask, or mask of voxels which contained sufficient EPI signal across all runs (this could be made more liberal via the -frac option, but fMRI does not need to be more liberal). Finally this unionized EPI mask is multiplied by the brain mask to construct an EPI-anatomy intersection mask.
Construct ATROPOS-class tissue masks, particularly for white matter and CSF via make_tissue_masks. First, start by thresholding the probabilistic fMRIprep tissue amsks (label-foo_probseg). Next erode the tissue masks via the -dilate 1 option of 3dmask_tool. It would be useful to check that the eroded CSF mask still contains voxels.
The function make_minimum_masks constructs a simple binary mask for voxels which contain >0 signal across all runs. This mask is used to avoid biasing the scaling step by a series of zeros.
The motion module contains the single function mot_files used to make the various motion files required by 3dDeconvolve. Using the fMRIprep confounds tsv file, the mean and derivative motion columns for 6dof are extracted as pandas dataframes. Mean and derivative motion columns are merely concatenated across runs and then written to 1D files for use with the -ortvec option. For censoring, fMRIprep records each motion event as a 1 in a unique binary column. These columns are extracted as their own dataframe and then summed in order to get a single column of volumes that require censoring. The values are inverted (0 = censor, 1 = good data), and volumes that precede a motion event are also excluded. The column is then written as the motion censor binary vector required by 3dDeconvolve. Finally, the binary values are re-inverted (verted?) to produce a binary vector (0 = good data, 1 = censor) for potential use as a baseline regressor. An info_<task>_censored_volumes.json file is written to the subject's directory which records the number of volumes excluded due to motion.
The process module contains functions that work with manipulating/processing MRI data. The functions are:
- blur_epi
- scale_epi
- reface
- resting_seed
- resting_metrics
The blur_epi function is an optional function used to blur or smooth EPI data. It is recommended for task-based analyses where ETAC is not employed as the group-level test (let ETAC control blurring). This function will calculated 1.5 x voxel size (dimension K), rounded up to the nearest integer for use as a FWHM Gaussian kernel. Each pre-processed EPI run data will be blurred by said kernel.
scale_epi will scale the timeseries of each voxel such that it is centered at the "MRI unit" of 100. This is done individually for each run.
The reface function is not incorporated in the data processing stream, but is instead used to anonymize T1-weighted files that will be hosted. This function supports defacing, refacing, and "refacing_plus" methods of dealing with face information via AFNI's drive @afni_refacer_run. Refacing is the default usage. Instead of writing output to project_dir/derivatives/afni, these output will be written to project_dir/derivatives/<reface_method>.
The resting_seed function creates a seed of size -srad 2 from supplied coordinates in MNI space. It then extracts the averaged timeseries for the seed. Next, a correlation matrix of resting state data with the averaged seed timeseries is produced via 3dTcorr1D, and this correlation matrix is converted for group-level testing via Fisher's Z-transform.
A series of metrics for resting state data are calculated via the resting_metrics function, because I felt nice. These include total SNR, global correlation, and Monte Carlo simulations of spurious clusters (ACF method).
The submit module contains two functions used for submitting sub-processes, because I cannot avoid writing bash, particularly with AFNI. It is mainly used by modules in the afni package, but also used by resources.fmriprep.preprocess and resources.ashs.hipseg. By default, both functions are hard-coded to load the modules afni-20.2.06 and c3d-1.0.0-gcc-8.2.0, but others can be added (e.g. resources.fmriprep.preprocess.run_freesurfer).
- submit_hpc_subprocess
- submit_hpc_sbatch
The function submit_hpc_subprocess submits bash syntax as a sub-process on the same scheduled resources as the parent job. subprocess.wait() is employed so the function returns stdout, stderr only once the process has finished.
The function submit_hpc_sbatch schedules a child process via sbatch, and the option --wait is employed so stdout/err are only returned after the job has finished. Child job names will be titled "1234foo", with "foo" being a short description of the step, so they are easily identifiable from the parent job (p1234). Stdout/err are named after the child job name, and saved to out_dir, usually the participants' working directory in /scratch or an sbatch_out directory within the subjects' session directory. Resource allotment is controlled via options.