diff --git a/aroma/aroma.py b/aroma/aroma.py index baf1cae..6c3e5e5 100644 --- a/aroma/aroma.py +++ b/aroma/aroma.py @@ -5,8 +5,8 @@ import os.path as op import shutil -import nibabel as nib import pandas as pd +from nilearn._utils import load_niimg from aroma import utils, features, _version @@ -29,6 +29,8 @@ def aroma_workflow( generate_plots=True, debug=False, quiet=False, + csf=None, + brain=None, mc_source="auto", ): """Run the AROMA workflow. @@ -169,7 +171,7 @@ def aroma_workflow( fsl_dir = op.join(os.environ["FSLDIR"], "bin", "") # Get TR of the fMRI data, if not specified if not TR: - in_img = nib.load(in_file) + in_img = load_niimg(in_file) TR = in_img.header.get_zooms()[3] # Check TR @@ -189,57 +191,37 @@ 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: - LGR.warning( - " - 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) + brain_img, csf_img, out_img, edge_img = utils.load_masks( + in_file, csf=csf, brain=brain + ) # Run ICA-AROMA LGR.info("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, brain_img, dim, TR + ) LGR.info("Step 2) Automatic classification of the components") LGR.info(" - 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) LGR.info(" - extracting the CSF & Edge fraction features") features_df = pd.DataFrame() features_df["edge_fract"], features_df["csf_fract"] = features.feature_spatial( - mel_IC_MNI + z_maps=mel_IC_MNI, + csf_mask=csf_img, + brain_mask=brain_img, + edge_mask=edge_img, + out_mask=out_img, ) LGR.info(" - extracting the Maximum RP correlation feature") - mel_mix = op.join(out_dir, "melodic.ica", "melodic_mix") mc = utils.load_motpars(mc, source=mc_source) - features_df["max_RP_corr"] = features.feature_time_series(mel_mix, mc) + features_df["max_RP_corr"] = features.feature_time_series(mixing, mc) LGR.info(" - extracting the High-frequency content feature") - mel_FT_mix = op.join(out_dir, "melodic.ica", "melodic_FTmix") - features_df["HFC"] = features.feature_frequency(mel_FT_mix, TR) + features_df["HFC"] = features.feature_frequency(mixing_FT, TR) LGR.info(" - classification") motion_ICs = utils.classification(features_df, out_dir) @@ -254,6 +236,6 @@ def aroma_workflow( if den_type != "no": LGR.info("Step 3) Data denoising") utils.denoising(fsl_dir, in_file, out_dir, - mel_mix, den_type, motion_ICs) + mixing, den_type, motion_ICs) LGR.info("Finished") diff --git a/aroma/cli/aroma.py b/aroma/cli/aroma.py index c8adf86..502e9a4 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", @@ -48,8 +52,9 @@ def _get_parser(): default=None, 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." ), ) @@ -130,6 +135,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", diff --git a/aroma/features.py b/aroma/features.py index 08f0764..fb1fd98 100644 --- a/aroma/features.py +++ b/aroma/features.py @@ -1,10 +1,9 @@ """Functions to calculate ICA-AROMA features for component classification.""" import logging -import os -import nibabel as nib import numpy as np from nilearn import image, masking +from nilearn._utils import load_niimg from . import utils @@ -157,7 +156,7 @@ def feature_frequency(mel_FT_mix, TR): return HFC -def feature_spatial(mel_IC): +def feature_spatial(*, z_maps, csf_mask, edge_mask, brain_mask, out_mask): """Extract the spatial feature scores. For each IC it determines the fraction of the mixture modeled thresholded @@ -169,6 +168,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 ------- @@ -180,20 +182,15 @@ def feature_spatial(mel_IC): mel_IC file """ # Get the number of ICs - mel_IC_img = nib.load(mel_IC) + mel_IC_img = load_niimg(z_maps) num_ICs = mel_IC_img.shape[3] - masks_dir = utils.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) for i in range(num_ICs): # Extract IC from the merged melodic_IC_thr2MNI2mm file - temp_IC = image.index_img(mel_IC, i) + temp_IC = image.index_img(mel_IC_img, i) # Change to absolute Z-values temp_IC = image.math_img("np.abs(img)", img=temp_IC) diff --git a/aroma/tests/test_features.py b/aroma/tests/test_features.py index efd9c63..5414146 100644 --- a/aroma/tests/test_features.py +++ b/aroma/tests/test_features.py @@ -1,8 +1,7 @@ +"""Tests for the aroma.features module.""" import numpy as np -import os.path as op -import pandas as pd -from aroma import features +from aroma import features, utils def test_feature_time_series(mel_mix, mc, max_correls): @@ -35,8 +34,19 @@ def test_feature_spatial(mel_IC, edgeFract, csfFract): np.random.seed(1) + # Get masks + brain_img, csf_img, out_img, edge_img = utils.load_masks( + mel_IC, csf=None, brain=None + ) + # Run feature_spatial - edge_fract, csf_fract = features.feature_spatial(mel_IC) + edge_fract, csf_fract = features.feature_spatial( + z_maps=mel_IC, + csf_mask=csf_img, + brain_mask=brain_img, + edge_mask=edge_img, + out_mask=out_img, + ) # Read features csv edgeFract = np.load(edgeFract) diff --git a/aroma/utils.py b/aroma/utils.py index fd5c0a4..2779b1e 100644 --- a/aroma/utils.py +++ b/aroma/utils.py @@ -1,4 +1,5 @@ """Utility functions for ICA-AROMA.""" +import datetime import logging import os import os.path as op @@ -8,10 +9,81 @@ import numpy as np import pandas as pd from nilearn import image, masking +from nilearn._utils import load_niimg +from scipy import ndimage LGR = logging.getLogger(__name__) +def load_masks(in_file, csf=None, brain=None): + """Load masks based on available inputs. + + Parameters + ---------- + in_file + csf + brain + + Returns + ------- + brain_img + csf_img + edge_img + out_img + """ + if csf is None: + LGR.info( + "No CSF TPM/mask provided. Assuming standard space data and 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: + if brain is None: + LGR.info("No brain mask provided. Using compute_epi_mask.") + brain_img = masking.compute_epi_mask(in_file) + + csf_img, edge_img, out_img = derive_native_masks(in_file, csf, brain_img) + + return brain_img, csf_img, edge_img, out_img + + +def derive_native_masks(in_file, csf, brain): + """Estimate or load necessary masks based on inputs. + + Parameters + ---------- + in_file : str + 4D EPI file. + csf : str or img_like + CSF tissue probability map or binary mask. + brain : str or img_like + Brain mask. + + Returns + ------- + masks : dict + Dictionary with the different masks as img_like objects. + """ + out_img = image.math_img("1 - mask", mask=brain) + csf_img = load_niimg(csf) + csf_data = csf_img.get_fdata() + if len(np.unique(csf_data)) == 2: + LGR.info("CSF mask provided. Inferring other masks.") + else: + LGR.info("CSF TPM provided. Inferring CSF and other masks.") + csf_img = image.math_img("csf >= 0.95", csf=csf_img) + gmwm_img = image.math_img("(brain - csf) > 0", brain=brain, 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) + + return csf_img, edge_img, out_img + + def runICA(fsl_dir, in_file, out_dir, mel_dir_in, mask, dim, TR): """Run MELODIC and merge the thresholded ICs into a single 4D nifti file. @@ -34,6 +106,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 @@ -41,10 +122,17 @@ 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 = op.join( + out_dir, + "temp_{}.nii.gz" + ).format(datetime.datetime.now().strftime('%Y-%m-%dT%H%M%S')) + 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") 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, @@ -102,11 +190,11 @@ def runICA(fsl_dir, in_file, out_dir, mel_dir_in, mask, dim, TR): "{0} --in={1} --outdir={2} --mask={3} --dim={4} " "--Ostats --nobet --mmthresh=0.5 --report " "--tr={5}" - ).format(op.join(fsl_dir, "melodic"), in_file, mel_dir, mask, dim, TR) + ).format(op.join(fsl_dir, "melodic"), in_file, mel_dir, temp_file, dim, TR) os.system(melodic_command) # Get number of components - mel_IC_img = nib.load(mel_IC) + mel_IC_img = load_niimg(mel_IC) nr_ICs = mel_IC_img.shape[3] # Merge mixture modeled thresholded spatial maps. Note! In case that @@ -119,7 +207,7 @@ def runICA(fsl_dir, in_file, out_dir, mel_dir_in, mask, dim, TR): z_temp = op.join(mel_dir, "stats", "thresh_zstat{0}.nii.gz".format(i)) # Get number of volumes in component's thresholded image - z_temp_img = nib.load(z_temp) + z_temp_img = load_niimg(z_temp) if z_temp_img.ndim == 4: len_IC = z_temp_img.shape[3] # Extract last spatial map within the thresh_zstat file @@ -138,120 +226,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) - - -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", - ] - ) - ) + os.remove(temp_file) + return mel_IC_thr, mel_IC_mix, mel_FT_mix def cross_correlation(a, b): @@ -377,7 +353,7 @@ def denoising(fsl_dir, in_file, out_dir, mixing, den_type, den_idx): motion_components = mixing[:, den_idx] # Create a fake mask to make it easier to reshape the full data to 2D - img = nib.load(in_file) + img = load_niimg(in_file) full_mask = nib.Nifti1Image(np.ones(img.shape[:3], int), img.affine) data = masking.apply_mask(img, full_mask) # T x S