From d2d1d70f05157f9fcbaca4b194d5a4a45625c433 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 7 Nov 2020 11:53:30 -0500 Subject: [PATCH 1/9] Add function for deriving masks. --- aroma/utils.py | 124 ++++++++++++++----------------------------------- 1 file changed, 34 insertions(+), 90 deletions(-) diff --git a/aroma/utils.py b/aroma/utils.py index bc4feb1..86189e9 100644 --- a/aroma/utils.py +++ b/aroma/utils.py @@ -6,6 +6,40 @@ import nibabel as nib import numpy as np from nilearn import image, masking +from scipy import ndimage + + +def derive_masks(in_file, csf=None): + if csf is None: + print("No CSF TPM/mask provided. Using packaged masks.") + mask_dir = get_resource_path() + csf_img = nib.load(op.join(mask_dir, 'mask_csf.nii.gz')) + out_img = nib.load(op.join(mask_dir, 'mask_out.nii.gz')) + brain_img = image.math_img('1 - mask', mask=out_img) + edge_img = nib.load(op.join(mask_dir, "mask_edge.nii.gz")) + else: + brain_img = masking.compute_epi_mask(in_file) + out_img = image.math_img('1 - mask', mask=brain_img) + csf_img = nib.load(csf) + csf_data = csf_img.get_fdata() + if len(np.unique(csf_data)) == 2: + print("CSF mask provided. Inferring other masks.") + else: + print("CSF TPM provided. Inferring CSF and other masks.") + csf_img = image.math_img('csf >= 0.3', csf=csf_img) + gmwm_img = image.math_img('(brain - csf) > 0', brain=brain_img, csf=csf_img) + gmwm_data = gmwm_img.get_fdata() + eroded_data = ndimage.binary_erosion(gmwm_data, iterations=4) + edge_data = gmwm_data - eroded_data + edge_img = nib.Nifti1Image(edge_data, gmwm_img.affine, header=gmwm_img.header) + + masks = { + "brain": brain_img, + "csf": csf_img, + "edge": edge_img, + "out": out_img, + } + return masks def runICA(fsl_dir, in_file, out_dir, mel_dir_in, mask, dim, TR): @@ -139,96 +173,6 @@ def runICA(fsl_dir, in_file, out_dir, mel_dir_in, mask, dim, TR): zstat_4d_img.to_filename(mel_IC_thr) -def register2MNI(fsl_dir, in_file, out_file, affmat, warp): - """Register an image (or time-series of images) to MNI152 T1 2mm. - - If no affmat is defined, it only warps (i.e. it assumes that the data has - been registered to the structural scan associated with the warp-file - already). If no warp is defined either, it only resamples the data to 2mm - isotropic if needed (i.e. it assumes that the data has been registered to - a MNI152 template). In case only an affmat file is defined, it assumes that - the data has to be linearly registered to MNI152 (i.e. the user has a - reason not to use non-linear registration on the data). - - Parameters - ---------- - fsl_dir : str - Full path of the bin-directory of FSL - in_file : str - Full path to the data file (nii.gz) which has to be registerd to - MNI152 T1 2mm - out_file : str - Full path of the output file - affmat : str - Full path of the mat file describing the linear registration (if data - is still in native space) - warp : str - Full path of the warp file describing the non-linear registration (if - data has not been registered to MNI152 space yet) - - Output - ------ - melodic_IC_mm_MNI2mm.nii.gz : merged file containing the mixture modeling - thresholded Z-statistical maps registered to - MNI152 2mm - """ - # Define the MNI152 T1 2mm template - fslnobin = fsl_dir.rsplit('/', 2)[0] - ref = op.join(fslnobin, 'data', 'standard', 'MNI152_T1_2mm_brain.nii.gz') - - # If the no affmat- or warp-file has been specified, assume that the data - # is already in MNI152 space. In that case only check if resampling to - # 2mm is needed - if not affmat and not warp: - in_img = nib.load(in_file) - # Get 3D voxel size - pixdim1, pixdim2, pixdim3 = in_img.header.get_zooms()[:3] - - # If voxel size is not 2mm isotropic, resample the data, otherwise - # copy the file - if (pixdim1 != 2) or (pixdim2 != 2) or (pixdim3 != 2): - os.system(' '.join([op.join(fsl_dir, 'flirt'), - ' -ref ' + ref, - ' -in ' + in_file, - ' -out ' + out_file, - ' -applyisoxfm 2 -interp trilinear'])) - else: - os.copyfile(in_file, out_file) - - # If only a warp-file has been specified, assume that the data has already - # been registered to the structural scan. In that case apply the warping - # without a affmat - elif not affmat and warp: - # Apply warp - os.system(' '.join([op.join(fsl_dir, 'applywarp'), - '--ref=' + ref, - '--in=' + in_file, - '--out=' + out_file, - '--warp=' + warp, - '--interp=trilinear'])) - - # If only a affmat-file has been specified perform affine registration to - # MNI - elif affmat and not warp: - os.system(' '.join([op.join(fsl_dir, 'flirt'), - '-ref ' + ref, - '-in ' + in_file, - '-out ' + out_file, - '-applyxfm -init ' + affmat, - '-interp trilinear'])) - - # If both a affmat- and warp-file have been defined, apply the warping - # accordingly - else: - os.system(' '.join([op.join(fsl_dir, 'applywarp'), - '--ref=' + ref, - '--in=' + in_file, - '--out=' + out_file, - '--warp=' + warp, - '--premat=' + affmat, - '--interp=trilinear'])) - - def cross_correlation(a, b): """Perform cross-correlations between columns of two matrices. From d6ac4fe43c1621586c3800ead489548818c86d4b Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 7 Nov 2020 11:56:59 -0500 Subject: [PATCH 2/9] Run black. --- aroma/utils.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/aroma/utils.py b/aroma/utils.py index 86189e9..ddf5a19 100644 --- a/aroma/utils.py +++ b/aroma/utils.py @@ -13,21 +13,21 @@ def derive_masks(in_file, csf=None): if csf is None: print("No CSF TPM/mask provided. Using packaged masks.") mask_dir = get_resource_path() - csf_img = nib.load(op.join(mask_dir, 'mask_csf.nii.gz')) - out_img = nib.load(op.join(mask_dir, 'mask_out.nii.gz')) - brain_img = image.math_img('1 - mask', mask=out_img) + csf_img = nib.load(op.join(mask_dir, "mask_csf.nii.gz")) + out_img = nib.load(op.join(mask_dir, "mask_out.nii.gz")) + brain_img = image.math_img("1 - mask", mask=out_img) edge_img = nib.load(op.join(mask_dir, "mask_edge.nii.gz")) else: brain_img = masking.compute_epi_mask(in_file) - out_img = image.math_img('1 - mask', mask=brain_img) + out_img = image.math_img("1 - mask", mask=brain_img) csf_img = nib.load(csf) csf_data = csf_img.get_fdata() if len(np.unique(csf_data)) == 2: print("CSF mask provided. Inferring other masks.") else: print("CSF TPM provided. Inferring CSF and other masks.") - csf_img = image.math_img('csf >= 0.3', csf=csf_img) - gmwm_img = image.math_img('(brain - csf) > 0', brain=brain_img, csf=csf_img) + csf_img = image.math_img("csf >= 0.3", csf=csf_img) + gmwm_img = image.math_img("(brain - csf) > 0", brain=brain_img, csf=csf_img) gmwm_data = gmwm_img.get_fdata() eroded_data = ndimage.binary_erosion(gmwm_data, iterations=4) edge_data = gmwm_data - eroded_data From e382be78a1bd731232d549cc24f3d8e0886ec1bd Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 7 Nov 2020 12:01:39 -0500 Subject: [PATCH 3/9] Use 95% probability threshold for CSF mask. --- aroma/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/aroma/utils.py b/aroma/utils.py index ddf5a19..b3b3396 100644 --- a/aroma/utils.py +++ b/aroma/utils.py @@ -26,7 +26,7 @@ def derive_masks(in_file, csf=None): print("CSF mask provided. Inferring other masks.") else: print("CSF TPM provided. Inferring CSF and other masks.") - csf_img = image.math_img("csf >= 0.3", csf=csf_img) + csf_img = image.math_img("csf >= 0.95", csf=csf_img) gmwm_img = image.math_img("(brain - csf) > 0", brain=brain_img, csf=csf_img) gmwm_data = gmwm_img.get_fdata() eroded_data = ndimage.binary_erosion(gmwm_data, iterations=4) From 31c4c741627387bde1e02abcd01c4820a338e670 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 7 Nov 2020 12:05:40 -0500 Subject: [PATCH 4/9] Add docstring. --- aroma/utils.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/aroma/utils.py b/aroma/utils.py index b3b3396..079a22a 100644 --- a/aroma/utils.py +++ b/aroma/utils.py @@ -10,6 +10,21 @@ def derive_masks(in_file, csf=None): + """Estimate or load necessary masks based on inputs. + + Parameters + ---------- + in_file : str + 4D EPI file. + csf : None or str, optional + CSF tissue probability map or binary mask. + If None, use precomputed, standard space masks packaged with AROMA. + + Returns + ------- + masks : dict + Dictionary with the different masks as img_like objects. + """ if csf is None: print("No CSF TPM/mask provided. Using packaged masks.") mask_dir = get_resource_path() From 63abc6e0283415443104ed028afb455e2a5b33c4 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 7 Nov 2020 12:34:57 -0500 Subject: [PATCH 5/9] Use masks. --- aroma/aroma.py | 42 ++++++++---------------------------------- aroma/features.py | 16 +++++++--------- 2 files changed, 15 insertions(+), 43 deletions(-) diff --git a/aroma/aroma.py b/aroma/aroma.py index 32f9444..518f01b 100644 --- a/aroma/aroma.py +++ b/aroma/aroma.py @@ -151,53 +151,27 @@ def aroma_workflow( # Define/create mask. Either by making a copy of the specified mask, or by # creating a new one. - new_mask = op.join(out_dir, "mask.nii.gz") - if mask: - shutil.copyfile(mask, new_mask) - elif in_feat and op.isfile(op.join(in_feat, "example_func.nii.gz")): - # If a Feat directory is specified, and an example_func is present use - # example_func to create a mask - bet_command = "{0} {1} {2} -f 0.3 -n -m -R".format( - op.join(fsl_dir, "bet"), - op.join(in_feat, "example_func.nii.gz"), - op.join(out_dir, "bet"), - ) - os.system(bet_command) - os.rename(op.join(out_dir, "bet_mask.nii.gz"), new_mask) - if op.isfile(op.join(out_dir, "bet.nii.gz")): - os.remove(op.join(out_dir, "bet.nii.gz")) - else: - if in_feat: - print( - " - No example_func was found in the Feat directory. " - "A mask will be created including all voxels with varying " - "intensity over time in the fMRI data. Please check!\n" - ) - math_command = "{0} {1} -Tstd -bin {2}".format( - op.join(fsl_dir, "fslmaths"), in_file, new_mask - ) - os.system(math_command) + masks = utils.derive_masks(in_file, csf=None) # Run ICA-AROMA print("Step 1) MELODIC") - utils.runICA(fsl_dir, in_file, out_dir, mel_dir, new_mask, dim, TR) + component_maps, mixing, mixing_FT = utils.runICA( + fsl_dir, in_file, out_dir, mel_dir, masks["brain"], dim, TR + ) print("Step 2) Automatic classification of the components") print(" - registering the spatial maps to MNI") - mel_IC = op.join(out_dir, "melodic_IC_thr.nii.gz") mel_IC_MNI = op.join(out_dir, "melodic_IC_thr_MNI2mm.nii.gz") - utils.register2MNI(fsl_dir, mel_IC, mel_IC_MNI, affmat, warp) + utils.register2MNI(fsl_dir, component_maps, mel_IC_MNI, affmat, warp) print(" - extracting the CSF & Edge fraction features") - edge_fract, csf_fract = features.feature_spatial(mel_IC_MNI) + edge_fract, csf_fract = features.feature_spatial(mel_IC_MNI, masks) print(" - extracting the Maximum RP correlation feature") - mel_mix = op.join(out_dir, "melodic.ica", "melodic_mix") - max_RP_corr = features.feature_time_series(mel_mix, mc) + max_RP_corr = features.feature_time_series(mixing, mc) print(" - extracting the High-frequency content feature") - mel_FT_mix = op.join(out_dir, "melodic.ica", "melodic_FTmix") - HFC = features.feature_frequency(mel_FT_mix, TR) + HFC = features.feature_frequency(mixing_FT, TR) print(" - classification") motion_ICs = utils.classification(out_dir, max_RP_corr, edge_fract, HFC, csf_fract) diff --git a/aroma/features.py b/aroma/features.py index 9257635..95c99e0 100644 --- a/aroma/features.py +++ b/aroma/features.py @@ -151,7 +151,7 @@ def feature_frequency(mel_FT_mix, TR): return HFC -def feature_spatial(mel_IC): +def feature_spatial(mel_IC, masks): """Extract the spatial feature scores. For each IC it determines the fraction of the mixture modeled thresholded @@ -163,6 +163,9 @@ def feature_spatial(mel_IC): mel_IC : str Full path of the nii.gz file containing mixture-modeled thresholded (p<0.5) Z-maps, registered to the MNI152 2mm template + masks : dict + Dictionary of masks. Keys are mask names and values are img_like + objects. Returns ------- @@ -177,11 +180,6 @@ def feature_spatial(mel_IC): mel_IC_img = nib.load(mel_IC) num_ICs = mel_IC_img.shape[3] - masks_dir = get_resource_path() - csf_mask = os.path.join(masks_dir, "mask_csf.nii.gz") - edge_mask = os.path.join(masks_dir, "mask_edge.nii.gz") - out_mask = os.path.join(masks_dir, "mask_out.nii.gz") - # Loop over ICs edge_fract = np.zeros(num_ICs) csf_fract = np.zeros(num_ICs) @@ -203,17 +201,17 @@ def feature_spatial(mel_IC): # Get sum of Z-values of the voxels located within the CSF # (calculate via the mean and number of non-zero voxels) - csf_data = masking.apply_mask(temp_IC, csf_mask) + csf_data = masking.apply_mask(temp_IC, masks["csf"]) csf_sum = np.sum(csf_data) # Get sum of Z-values of the voxels located within the Edge # (calculate via the mean and number of non-zero voxels) - edge_data = masking.apply_mask(temp_IC, edge_mask) + edge_data = masking.apply_mask(temp_IC, masks["edge"]) edge_sum = np.sum(edge_data) # Get sum of Z-values of the voxels located outside the brain # (calculate via the mean and number of non-zero voxels) - out_data = masking.apply_mask(temp_IC, out_mask) + out_data = masking.apply_mask(temp_IC, masks["out"]) out_sum = np.sum(out_data) # Determine edge and CSF fraction From 6d45c321350b7825be337b9155c2127dcfdbd447 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 7 Nov 2020 12:35:13 -0500 Subject: [PATCH 6/9] Return ICA outputs. --- aroma/utils.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/aroma/utils.py b/aroma/utils.py index 079a22a..48226af 100644 --- a/aroma/utils.py +++ b/aroma/utils.py @@ -79,6 +79,15 @@ def runICA(fsl_dir, in_file, out_dir, mel_dir_in, mask, dim, TR): TR : float TR (in seconds) of the fMRI data + Returns + ------- + mel_IC_thr : str + Path to 4D nifti file containing thresholded component maps. + mel_IC_mix : str + Path to mixing matrix. + mel_FT_mix : str + Path to component power spectrum file. + Output ------ melodic.ica/: MELODIC directory @@ -90,6 +99,7 @@ def runICA(fsl_dir, in_file, out_dir, mel_dir_in, mask, dim, TR): mel_dir = op.join(out_dir, 'melodic.ica') mel_IC = op.join(mel_dir, 'melodic_IC.nii.gz') mel_IC_mix = op.join(mel_dir, 'melodic_mix') + mel_FT_mix = op.join(mel_dir, "melodic_FTmix") mel_IC_thr = op.join(out_dir, 'melodic_IC_thr.nii.gz') # When a MELODIC directory is specified, @@ -186,6 +196,7 @@ def runICA(fsl_dir, in_file, out_dir, mel_dir_in, mask, dim, TR): "stat * mask[:, :, :, None]", stat=zstat_4d_img, mask=mask ) zstat_4d_img.to_filename(mel_IC_thr) + return mel_IC_thr, mel_IC_mix, mel_FT_mix def cross_correlation(a, b): From 14761a197854e7fa1bebc3020a47456f95ab6f49 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 7 Nov 2020 12:45:00 -0500 Subject: [PATCH 7/9] Fix bugs. --- aroma/aroma.py | 2 +- aroma/features.py | 4 +--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/aroma/aroma.py b/aroma/aroma.py index 518f01b..6106af3 100644 --- a/aroma/aroma.py +++ b/aroma/aroma.py @@ -185,6 +185,6 @@ def aroma_workflow( if den_type != "no": print("Step 3) Data denoising") - utils.denoising(fsl_dir, in_file, out_dir, mel_mix, den_type, motion_ICs) + utils.denoising(fsl_dir, in_file, out_dir, mixing, den_type, motion_ICs) print("Finished") diff --git a/aroma/features.py b/aroma/features.py index 95c99e0..400acb5 100644 --- a/aroma/features.py +++ b/aroma/features.py @@ -1,11 +1,9 @@ """Functions to calculate ICA-AROMA features for component classification.""" -import os - import nibabel as nib import numpy as np from nilearn import image, masking -from .utils import cross_correlation, get_resource_path +from .utils import cross_correlation def feature_time_series(mel_mix, mc): From 6a468369bfa91bf1bc2166d1307665ea0736a273 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sun, 8 Nov 2020 09:34:43 -0500 Subject: [PATCH 8/9] Write mask to temporary file for MELODIC. --- aroma/utils.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/aroma/utils.py b/aroma/utils.py index 48226af..5501dd0 100644 --- a/aroma/utils.py +++ b/aroma/utils.py @@ -2,6 +2,7 @@ import os import os.path as op import shutil +from tempfile import mkstemp import nibabel as nib import numpy as np @@ -95,6 +96,9 @@ def runICA(fsl_dir, in_file, out_dir, mel_dir_in, mask, dim, TR): thresholded Z-statistical maps located in melodic.ica/stats/ """ + temp_file = mkstemp(suffix=".nii.gz") + mask.to_filename(temp_file) + # Define the 'new' MELODIC directory and predefine some associated files mel_dir = op.join(out_dir, 'melodic.ica') mel_IC = op.join(mel_dir, 'melodic_IC.nii.gz') @@ -157,7 +161,7 @@ def runICA(fsl_dir, in_file, out_dir, mel_dir_in, mask, dim, TR): op.join(fsl_dir, 'melodic'), in_file, mel_dir, - mask, + temp_file, dim, TR ) @@ -196,6 +200,8 @@ def runICA(fsl_dir, in_file, out_dir, mel_dir_in, mask, dim, TR): "stat * mask[:, :, :, None]", stat=zstat_4d_img, mask=mask ) zstat_4d_img.to_filename(mel_IC_thr) + + os.remove(temp_file) return mel_IC_thr, mel_IC_mix, mel_FT_mix From ab8c6041b4d2197ef16e6f76330601743d209527 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Mon, 9 Nov 2020 09:11:25 -0500 Subject: [PATCH 9/9] Add CSF argument. --- aroma/aroma.py | 1 + aroma/cli/aroma.py | 22 +++++++++++++++++++--- 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/aroma/aroma.py b/aroma/aroma.py index 6106af3..2c2c456 100644 --- a/aroma/aroma.py +++ b/aroma/aroma.py @@ -22,6 +22,7 @@ def aroma_workflow( TR=None, overwrite=False, generate_plots=True, + csf=None, ): """Run the AROMA workflow. diff --git a/aroma/cli/aroma.py b/aroma/cli/aroma.py index 0cd6bc7..9472719 100644 --- a/aroma/cli/aroma.py +++ b/aroma/cli/aroma.py @@ -38,7 +38,11 @@ def _get_parser(): dest="in_file", type=lambda x: is_valid_file(parser, x), required=False, - help="Input file name of fMRI data (.nii.gz)", + help=( + "Input file name of fMRI data (.nii.gz). " + "This file may be in standard (MNI152) or native space, with some restrictions. " + "The data should be smoothed prior to running AROMA." + ), ) inputs.add_argument( "-f", @@ -47,8 +51,9 @@ def _get_parser(): required=False, type=lambda x: is_valid_path(parser, x), help=( - "Feat directory name (Feat should have been run without temporal " - "filtering and including registration to MNI152)" + "Path to FSL FEAT directory. " + "FEAT should have been run without temporal filtering, " + "but with registration to MNI152 space." ), ) @@ -113,6 +118,17 @@ def _get_parser(): # Optional options optoptions = parser.add_argument_group("Optional arguments") optoptions.add_argument("-tr", dest="TR", help="TR in seconds", type=float) + optoptions.add_argument( + "--csf", + dest="csf", + type=lambda x: is_valid_file(parser, x), + default=None, + help=( + "Path to a cerebrospinal fluid (CSF) mask or tissue probability map. " + "If this file is not provided, then data are assumed to be in standard space, " + "and prepackaged masks will be used instead." + ), + ) optoptions.add_argument( "-den", dest="den_type",