diff --git a/docs/workflows.rst b/docs/workflows.rst index b48f2476c..0b1b296fd 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -942,29 +942,27 @@ Amplitude of low-frequency fluctuation (ALFF) is a measure that ostensibly local spontaneous neural activity in resting-state BOLD data. It is calculated by the following: -1. The ``denoised, interpolated BOLD`` is passed along to the ALFF workflow. -2. If censoring+interpolation was performed, then the interpolated time series is censored at this - point. -3. Voxel-wise BOLD time series are normalized (mean-centered and scaled to unit standard deviation) +1. The ``denoised BOLD`` is passed along to the ALFF workflow. +2. Voxel-wise BOLD time series are normalized (mean-centered and scaled to unit standard deviation) over time. This will ensure that the power spectrum from ``periodogram`` and ``lombscargle`` are roughly equivalent. -4. The power spectrum and associated frequencies are estimated from the BOLD data. +3. The power spectrum and associated frequencies are estimated from the BOLD data. - If censoring+interpolation was not performed, then this uses :func:`scipy.signal.periodogram`. - If censoring+interpolation was performed, then this uses :func:`scipy.signal.lombscargle`. -5. The square root of the power spectrum is calculated. -6. The power spectrum values corresponding to the frequency range retained by the +4. The square root of the power spectrum is calculated. +5. The power spectrum values corresponding to the frequency range retained by the temporal filtering step are extracted from the full power spectrum. -7. The mean of the within-band power spectrum is calculated and multiplied by 2. -8. The ALFF value is multiplied by the standard deviation of the voxel-wise +6. The mean of the within-band power spectrum is calculated and multiplied by 2. +7. The ALFF value is multiplied by the standard deviation of the voxel-wise ``denoised, interpolated BOLD`` time series. This brings ALFF back to its original scale, as if the time series was not normalized. ALFF will only be calculated if the bandpass filter is enabled (i.e., if the ``--disable-bandpass-filter`` flag is not used). -Smoothed ALFF derivatives will also be generated if the ``--smoothing`` flag is used. +ALFF will also be calculated from the smoothed, denoised BOLD if the ``--smoothing`` flag is used. ReHo @@ -976,8 +974,12 @@ Regional Homogeneity (ReHo) is a measure of local temporal uniformity in the BOL Greater ReHo values correspond to greater synchrony among BOLD activity patterns measured in a local neighborhood of voxels, with neighborhood size determined by a user-specified radius of voxels. ReHo is calculated as the coefficient of concordance among all voxels in a sphere centered on the target voxel. -For NIfTIs, ReHo is always calculated via AFNI’s 3dReho with 27 voxels in each neighborhood, using Kendall's coefficient of concordance (KCC). -For CIFTIs, the left and right hemisphere are extracted into GIFTI format via Connectome Workbench’s CIFTISeparateMetric. Next, the mesh adjacency matrix is obtained,and Kendall's coefficient of concordance (KCC) is calculated, with each vertex having four neighbors. +For NIfTIs, ReHo is always calculated via AFNI's 3dReho with 27 voxels in each neighborhood, +using Kendall's coefficient of concordance (KCC). +For CIFTIs, the left and right hemisphere are extracted into GIFTI format via +Connectome Workbench's CIFTISeparateMetric. +Next, the mesh adjacency matrix is obtained, and Kendall's coefficient of concordance (KCC) is calculated, +with each vertex having 5-6 neighbors. For subcortical voxels in the CIFTIs, 3dReho is used with the same parameters that are used for NIfTIs. diff --git a/xcp_d/tests/test_utils_restingstate.py b/xcp_d/tests/test_utils_restingstate.py index 11efc2d14..6d67fa141 100644 --- a/xcp_d/tests/test_utils_restingstate.py +++ b/xcp_d/tests/test_utils_restingstate.py @@ -50,6 +50,7 @@ def test_compute_alff(ds001419_data): # Now try with a sample mask sample_mask = np.ones(bold_data.shape[1], dtype=bool) sample_mask[20:30] = False + bold_data = bold_data[:, sample_mask] alff2 = restingstate.compute_alff_chunk((bold_data, 0.1, 0.01, TR, sample_mask)) diff --git a/xcp_d/tests/test_workflows_metrics.py b/xcp_d/tests/test_workflows_metrics.py index 1783ad8d5..87ab598d7 100644 --- a/xcp_d/tests/test_workflows_metrics.py +++ b/xcp_d/tests/test_workflows_metrics.py @@ -52,6 +52,7 @@ def test_nifti_alff(ds001419_data, tmp_path_factory): alff_wf.base_dir = tempdir alff_wf.inputs.inputnode.bold_mask = bold_mask alff_wf.inputs.inputnode.denoised_bold = bold_file + alff_wf.inputs.inputnode.smoothed_bold = bold_file alff_wf = clean_datasinks(alff_wf) compute_alff_res = alff_wf.run() @@ -94,6 +95,7 @@ def test_nifti_alff(ds001419_data, tmp_path_factory): alff_wf.base_dir = tempdir alff_wf.inputs.inputnode.bold_mask = bold_mask alff_wf.inputs.inputnode.denoised_bold = filename + alff_wf.inputs.inputnode.smoothed_bold = filename alff_wf = clean_datasinks(alff_wf) compute_alff_res = alff_wf.run() nodes = get_nodes(compute_alff_res) @@ -140,6 +142,7 @@ def test_cifti_alff(ds001419_data, tmp_path_factory): alff_wf.base_dir = tempdir alff_wf.inputs.inputnode.bold_mask = bold_mask alff_wf.inputs.inputnode.denoised_bold = bold_file + alff_wf.inputs.inputnode.smoothed_bold = bold_file alff_wf = clean_datasinks(alff_wf) compute_alff_res = alff_wf.run() @@ -169,9 +172,16 @@ def test_cifti_alff(ds001419_data, tmp_path_factory): # Now let's compute ALFF for the new file and see how it compares tempdir = tmp_path_factory.mktemp('test_cifti_alff_02') + config.workflow.smoothing = 0 + alff_wf = metrics.init_alff_wf( + name_source=bold_file, + TR=TR, + mem_gb=mem_gbx, + ) alff_wf.base_dir = tempdir alff_wf.inputs.inputnode.bold_mask = bold_mask alff_wf.inputs.inputnode.denoised_bold = filename + alff_wf = clean_datasinks(alff_wf) compute_alff_res = alff_wf.run() nodes = get_nodes(compute_alff_res) diff --git a/xcp_d/utils/doc.py b/xcp_d/utils/doc.py index 779e4a06d..b16a6221e 100644 --- a/xcp_d/utils/doc.py +++ b/xcp_d/utils/doc.py @@ -113,7 +113,6 @@ smoothing : :obj:`float` The full width at half maximum (FWHM), in millimeters, of the Gaussian smoothing kernel that will be applied to the post-processed and denoised data. - ALFF and ReHo outputs will also be smoothing with this kernel. """ docdict['custom_confounds_folder'] = """ @@ -381,8 +380,8 @@ smoothed_denoised_bold : :obj:`str` Path to the censored, denoised, interpolated, filtered, re-censored, and smoothed BOLD file. This file is the result of denoising the censored preprocessed BOLD data, - followed by cubic spline interpolation, band-pass filtering, re-censoring, and spatial - smoothing. + followed by cubic spline interpolation, band-pass filtering, + re-censoring (if output_interpolated is False), and spatial smoothing. """ docdict['atlases'] = """ diff --git a/xcp_d/utils/restingstate.py b/xcp_d/utils/restingstate.py index 8fc9f570a..f1b6bc433 100644 --- a/xcp_d/utils/restingstate.py +++ b/xcp_d/utils/restingstate.py @@ -130,7 +130,7 @@ def compute_alff(*, data_matrix, low_pass, high_pass, TR, sample_mask): Parameters ---------- data_matrix : numpy.ndarray - data matrix points by timepoints + data matrix points by timepoints, after optional censoring low_pass : float low pass frequency in Hz high_pass : float @@ -161,6 +161,15 @@ def compute_alff(*, data_matrix, low_pass, high_pass, TR, sample_mask): fs = 1 / TR # sampling frequency n_voxels, n_volumes = data_matrix.shape + if sample_mask is not None: + if sample_mask.sum() != n_volumes: + raise ValueError(f'{sample_mask.sum()} != {n_volumes}') + + time_arr = np.arange(0, sample_mask.size * TR, TR) + time_arr = time_arr[sample_mask] + if time_arr.size != n_volumes: + raise ValueError(f'{time_arr.size} != {n_volumes}') + alff = np.zeros(n_voxels) for i_voxel in range(n_voxels): voxel_data = data_matrix[i_voxel, :] @@ -189,7 +198,7 @@ def compute_alff(*, data_matrix, low_pass, high_pass, TR, sample_mask): angular_frequencies = 2 * np.pi * frequencies_hz power_spectrum = signal.lombscargle( time_arr, - voxel_data_censored, + voxel_data, angular_frequencies, normalize=True, ) diff --git a/xcp_d/workflows/bold/cifti.py b/xcp_d/workflows/bold/cifti.py index 298e3641f..e2685de74 100644 --- a/xcp_d/workflows/bold/cifti.py +++ b/xcp_d/workflows/bold/cifti.py @@ -292,7 +292,8 @@ def init_postprocess_cifti_wf( ('outputnode.temporal_mask', 'inputnode.temporal_mask'), ]), (denoise_bold_wf, alff_wf, [ - ('outputnode.denoised_interpolated_bold', 'inputnode.denoised_bold'), + ('outputnode.censored_denoised_bold', 'inputnode.denoised_bold'), + ('outputnode.smoothed_denoised_bold', 'inputnode.smoothed_bold'), ]), ]) # fmt:skip diff --git a/xcp_d/workflows/bold/metrics.py b/xcp_d/workflows/bold/metrics.py index c74359e47..b4383b797 100644 --- a/xcp_d/workflows/bold/metrics.py +++ b/xcp_d/workflows/bold/metrics.py @@ -5,23 +5,18 @@ from nipype.interfaces import utility as niu from nipype.pipeline import engine as pe from niworkflows.engine.workflows import LiterateWorkflow as Workflow -from templateflow.api import get as get_template from xcp_d import config from xcp_d.config import dismiss_hash from xcp_d.interfaces.bids import DerivativesDataSink -from xcp_d.interfaces.nilearn import Smooth from xcp_d.interfaces.plotting import PlotDenseCifti, PlotNifti from xcp_d.interfaces.restingstate import ComputeALFF, ReHoNamePatch, SurfaceReHo from xcp_d.interfaces.workbench import ( CiftiCreateDenseFromTemplate, CiftiSeparateMetric, CiftiSeparateVolumeAll, - CiftiSmooth, - FixCiftiIntent, ) from xcp_d.utils.doc import fill_doc -from xcp_d.utils.utils import fwhm2sigma @fill_doc @@ -62,8 +57,9 @@ def init_alff_wf( Inputs ------ denoised_bold - This is the ``filtered, interpolated, denoised BOLD``, - although interpolation is not necessary if the data were not originally censored. + This is the denoised BOLD after optional re-censoring. + smoothed_bold + Smoothed BOLD after optional re-censoring. bold_mask bold mask if bold is nifti temporal_mask @@ -125,6 +121,7 @@ def init_alff_wf( niu.IdentityInterface( fields=[ 'denoised_bold', + 'smoothed_bold', 'bold_mask', 'temporal_mask', # only used for CIFTI data if the anatomical workflow is enabled @@ -195,61 +192,62 @@ def init_alff_wf( ]) # fmt:skip if smoothing: # If we want to smooth - if file_format == 'nifti': - workflow.__desc__ = workflow.__desc__ + ( - ' The ALFF maps were smoothed with Nilearn using a Gaussian kernel ' - f'(FWHM={str(smoothing)} mm).' - ) - # Smooth via Nilearn - smooth_data = pe.Node( - Smooth(fwhm=smoothing), - name='niftismoothing', - ) - workflow.connect([ - (alff_compt, smooth_data, [('alff', 'in_file')]), - (smooth_data, outputnode, [('out_file', 'smoothed_alff')]) - ]) # fmt:skip - - else: # If cifti - workflow.__desc__ = workflow.__desc__ + ( - ' The ALFF maps were smoothed with the Connectome Workbench using a Gaussian ' - f'kernel (FWHM={str(smoothing)} mm).' - ) + workflow.__desc__ += ' ALFF was also calculated from the smoothed, denoised BOLD data.' + + # compute alff + compute_smoothed_alff = pe.Node( + ComputeALFF( + TR=TR, + low_pass=low_pass, + high_pass=high_pass, + n_threads=config.nipype.omp_nthreads, + ), + mem_gb=mem_gb['bold'], + name='compute_smoothed_alff', + n_procs=config.nipype.omp_nthreads, + ) + workflow.connect([ + (inputnode, compute_smoothed_alff, [ + ('smoothed_bold', 'in_file'), + ('bold_mask', 'mask'), + ('temporal_mask', 'temporal_mask'), + ]), + (compute_smoothed_alff, outputnode, [('alff', 'smoothed_alff')]) + ]) # fmt:skip - # Smooth via Connectome Workbench - sigma_lx = fwhm2sigma(smoothing) # Convert fwhm to standard deviation - # Get templates for each hemisphere - lh_midthickness = str( - get_template('fsLR', hemi='L', suffix='sphere', density='32k', raise_empty=True)[0] - ) - rh_midthickness = str( - get_template('fsLR', hemi='R', suffix='sphere', density='32k', raise_empty=True)[0] - ) - smooth_data = pe.Node( - CiftiSmooth( - sigma_surf=sigma_lx, - sigma_vol=sigma_lx, - direction='COLUMN', - right_surf=rh_midthickness, - left_surf=lh_midthickness, - num_threads=config.nipype.omp_nthreads, - ), - name='ciftismoothing', - mem_gb=mem_gb['bold'], - n_procs=config.nipype.omp_nthreads, - ) + # Plot the ALFF map + ds_report_smoothed_alff = pe.Node( + DerivativesDataSink( + source_file=name_source, + ), + name='ds_report_smoothed_alff', + run_without_submitting=False, + ) - # Always check the intent code in CiftiSmooth's output file - fix_cifti_intent = pe.Node( - FixCiftiIntent(), - name='fix_cifti_intent', - mem_gb=mem_gb['bold'], + if file_format == 'cifti': + plot_smoothed_alff = pe.Node( + PlotDenseCifti(base_desc='alff'), + name='plot_smoothed_alff', ) workflow.connect([ - (alff_compt, smooth_data, [('alff', 'in_file')]), - (smooth_data, fix_cifti_intent, [('out_file', 'in_file')]), - (fix_cifti_intent, outputnode, [('out_file', 'smoothed_alff')]), + (inputnode, plot_smoothed_alff, [ + ('lh_midthickness', 'lh_underlay'), + ('rh_midthickness', 'rh_underlay'), + ]), + (plot_smoothed_alff, ds_report_smoothed_alff, [('desc', 'desc')]), ]) # fmt:skip + ds_report_smoothed_alff.inputs.desc = 'alffSmoothedSurfacePlot' + else: + plot_smoothed_alff = pe.Node( + PlotNifti(name_source=name_source), + name='plot_smoothed_alff', + ) + ds_report_smoothed_alff.inputs.desc = 'alffSmoothedVolumetricPlot' + + workflow.connect([ + (compute_smoothed_alff, plot_smoothed_alff, [('alff', 'in_file')]), + (plot_smoothed_alff, ds_report_smoothed_alff, [('out_file', 'in_file')]), + ]) # fmt:skip return workflow diff --git a/xcp_d/workflows/bold/nifti.py b/xcp_d/workflows/bold/nifti.py index 23531509a..f526516fd 100644 --- a/xcp_d/workflows/bold/nifti.py +++ b/xcp_d/workflows/bold/nifti.py @@ -303,7 +303,8 @@ def init_postprocess_nifti_wf( ('outputnode.temporal_mask', 'inputnode.temporal_mask'), ]), (denoise_bold_wf, alff_wf, [ - ('outputnode.denoised_interpolated_bold', 'inputnode.denoised_bold'), + ('outputnode.censored_denoised_bold', 'inputnode.denoised_bold'), + ('outputnode.smoothed_denoised_bold', 'inputnode.smoothed_bold'), ]), ]) # fmt:skip