OSCR

MEPrep: A robust pipeline for multi-echo fMRI denoising and preprocessing.

Code ↔ Paper

6 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 6 matches
  1. [1] § Methods › Pre-denoising multi-echo fMRI data with probabilistic ICA ↔ MEPrep_Source/MEPrep_ICA-AROMA/ICA_AROMA_functions_preICA.py, lines 217–287 · score 0.92 · shifted versions, realignment parameter, maximum absolute, maxRPcorr, maximum correlations, motion parameters
  2. [2] § Methods › Pre-denoising multi-echo fMRI data with probabilistic ICA ↔ MEPrep_Source/MEPrep_ICA-AROMA/ICA_AROMA_functions_pureMotion.py, lines 213–283 · score 0.92 · shifted versions, realignment parameter, maximum absolute, maxRPcorr, maximum correlations, motion parameters
  3. [3] § Methods › Pre-denoising multi-echo fMRI data with probabilistic ICA ↔ MEPrep_Source/MEPrep_ICA-AROMA/ICA_AROMA_functions_preICA.py, lines 461–535 · score 0.68 · maxRPcorr, CSF fractions, ICA AROMA, HFC, score, edge
  4. [4] § Methods › Pre-denoising multi-echo fMRI data with probabilistic ICA ↔ MEPrep_Source/MEPrep_ICA-AROMA/ICA_AROMA_functions_pureMotion.py, lines 457–533 · score 0.68 · maxRPcorr, CSF fractions, ICA AROMA, HFC, score, edge
  5. [5] § Methods › The framework and implementation of MEPrep ↔ MEPrep_Source/fmriprep/workflows/bold/confounds.py, lines 626–769 · score 0.51 · CIFTI formats, alignment, FD, confound, WM, NIfTI
  6. [6] § Methods › The framework and implementation of MEPrep ↔ MEPrep_Source/meprep/src/workflows/bold/confounds.py, lines 670–745 · score 0.50 · CIFTI formats, pre denoising, alignment, confound, WM, multiecho

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 591 lines · 27 KB · no license · 2 matches

  1. #!/usr/bin/env python
  3. # Functions for ICA-AROMA v0.3 beta
  5. from __future__ import division
  6. from __future__ import print_function
  7. from future import standard_library
  8. standard_library.install_aliases()
  9. from builtins import str
  10. from builtins import range
  11. from past.utils import old_div
  12. import numpy as np
  13. import shutil
  15. def runICA(fslDir, inFile, outDir, melDirIn, mask, dim, TR, preICA_scme=False):
  16. """ This function runs MELODIC and merges the mixture modeled thresholded ICs into a single 4D nifti file
  18. Parameters
  19. ---------------------------------------------------------------------------------
  20. fslDir: Full path of the bin-directory of FSL
  21. inFile: Full path to the fMRI data file (nii.gz) on which MELODIC should be run
  22. outDir: Full path of the output directory
  23. melDirIn: Full path of the MELODIC directory in case it has been run before, otherwise define empty string
  24. mask: Full path of the mask to be applied during MELODIC
  25. dim: Dimensionality of ICA
  26. TR: TR (in seconds) of the fMRI data
  28. Output (within the requested output directory)
  29. ---------------------------------------------------------------------------------
  30. melodic.ica MELODIC directory
  31. melodic_IC_thr.nii.gz merged file containing the mixture modeling thresholded Z-statistical maps located in melodic.ica/stats/ """
  33. # Import needed modules
  34. import os
  35. import subprocess
  37. # Define the 'new' MELODIC directory and predefine some associated files
  38. melDir = os.path.join(outDir, 'melodic.ica')
  39. melIC = os.path.join(melDir, 'melodic_IC.nii.gz')
  40. melICmix = os.path.join(melDir, 'melodic_mix')
  41. melICthr = os.path.join(outDir, 'melodic_IC_thr.nii.gz')
  42. melICthrIn = os.path.join(melDirIn, 'melodic_IC_thr.nii.gz')
  44. # When a MELODIC directory is specified,
  45. # check whether all needed files are present.
  46. # Otherwise... run MELODIC again
  47. if len(melDir) != 0 and os.path.isfile(os.path.join(melDirIn, 'melodic_IC.nii.gz')) and os.path.isfile(os.path.join(melDirIn, 'melodic_FTmix')) and os.path.isfile(os.path.join(melDirIn, 'melodic_mix')):
  49. print(' - The existing/specified MELODIC directory will be used.')
  51. # If a 'stats' directory is present (contains thresholded spatial maps)
  52. # create a symbolic link to the MELODIC directory.
  53. # Otherwise create specific links and
  54. # run mixture modeling to obtain thresholded maps.
  55. if os.path.isdir(os.path.join(melDirIn, 'stats')):
  56. os.symlink(melDirIn, melDir)
  57. else:
  58. print(' - The MELODIC directory does not contain the required \'stats\' folder. Mixture modeling on the Z-statistical maps will be run.')
  60. # Create symbolic links to the items in the specified melodic directory
  61. os.makedirs(melDir)
  62. for item in os.listdir(melDirIn):
  63. os.symlink(os.path.join(melDirIn, item),
  64. os.path.join(melDir, item))
  66. # Run mixture modeling
  67. os.system(' '.join([os.path.join(fslDir, 'melodic'),
  68. '--in=' + melIC,
  69. '--ICs=' + melIC,
  70. '--mix=' + melICmix,
  71. '--outdir=' + melDir,
  72. '--Ostats --mmthresh=0.5']))
  74. else:
  75. # If a melodic directory was specified, display that it did not contain all files needed for ICA-AROMA (or that the directory does not exist at all)
  76. if len(melDirIn) != 0:
  77. if not os.path.isdir(melDirIn):
  78. print(' - The specified MELODIC directory does not exist. MELODIC will be run seperately.')
  79. else:
  80. print(' - The specified MELODIC directory does not contain the required files to run ICA-AROMA. MELODIC will be run seperately.')
  82. # Run MELODIC
  83. os.system(' '.join([os.path.join(fslDir, 'melodic'),
  84. '--in=' + inFile,
  85. '--outdir=' + melDir,
  86. '--mask=' + mask,
  87. '--dim=' + str(dim),
  88. '--Ostats --nobet --mmthresh=0.5 --report',
  89. '--tr=' + str(TR)]))
  91. # Get number of components
  92. cmd = ' '.join([os.path.join(fslDir, 'fslinfo'),
  93. melIC,
  94. '| grep dim4 | head -n1 | awk \'{print $2}\''])
  95. nrICs = int(float(subprocess.getoutput(cmd)))
  97. if preICA_scme:
  98. shutil.copyfile(melICthrIn, melICthr)
  99. else:
  100. # Merge mixture modeled thresholded spatial maps. Note! In case that mixture modeling did not converge, the file will contain two spatial maps. The latter being the results from a simple null hypothesis test. In that case, this map will have to be used (first one will be empty).
  101. for i in range(1, nrICs + 1):
  102. # Define thresholded zstat-map file
  103. zTemp = os.path.join(melDir, 'stats', 'thresh_zstat' + str(i) + '.nii.gz')
  104. cmd = ' '.join([os.path.join(fslDir, 'fslinfo'),
  105. zTemp,
  106. '| grep dim4 | head -n1 | awk \'{print $2}\''])
  107. lenIC = int(float(subprocess.getoutput(cmd)))
  109. # Define zeropad for this IC-number and new zstat file
  110. cmd = ' '.join([os.path.join(fslDir, 'zeropad'),
  111. str(i),
  112. '4'])
  113. ICnum = subprocess.getoutput(cmd)
  114. zstat = os.path.join(outDir, 'thr_zstat' + ICnum)
  116. # Extract last spatial map within the thresh_zstat file
  117. os.system(' '.join([os.path.join(fslDir, 'fslroi'),
  118. zTemp, # input
  119. zstat, # output
  120. str(lenIC - 1), # first frame
  121. '1'])) # number of frames
  123. # Merge and subsequently remove all mixture modeled Z-maps within the output directory
  124. os.system(' '.join([os.path.join(fslDir, 'fslmerge'),
  125. '-t', # concatenate in time
  126. melICthr, # output
  127. os.path.join(outDir, 'thr_zstat????.nii.gz')])) # inputs
  129. os.system('rm ' + os.path.join(outDir, 'thr_zstat????.nii.gz'))
  131. # Apply the mask to the merged file (in case a melodic-directory was predefined and run with a different mask)
  132. os.system(' '.join([os.path.join(fslDir, 'fslmaths'),
  133. melICthr,
  134. '-mas ' + mask,
  135. melICthr]))
  138. def register2MNI(fslDir, inFile, outFile, affmat, warp):
  139. """ This function registers 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 registerd 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).
  141. Parameters
  142. ---------------------------------------------------------------------------------
  143. fslDir: Full path of the bin-directory of FSL
  144. inFile: Full path to the data file (nii.gz) which has to be registerd to MNI152 T1 2mm
  145. outFile: Full path of the output file
  146. affmat: Full path of the mat file describing the linear registration (if data is still in native space)
  147. warp: Full path of the warp file describing the non-linear registration (if data has not been registered to MNI152 space yet)
  149. Output (within the requested output directory)
  150. ---------------------------------------------------------------------------------
  151. melodic_IC_mm_MNI2mm.nii.gz merged file containing the mixture modeling thresholded Z-statistical maps registered to MNI152 2mm """
  154. # Import needed modules
  155. import os
  156. import subprocess
  158. # Define the MNI152 T1 2mm template
  159. fslnobin = fslDir.rsplit('/', 2)[0]
  160. ref = os.path.join(fslnobin, 'data', 'standard', 'MNI152_T1_2mm_brain.nii.gz')
  162. # 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
  163. if (len(affmat) == 0) and (len(warp) == 0):
  164. # Get 3D voxel size
  165. pixdim1 = float(subprocess.getoutput('%sfslinfo %s | grep pixdim1 | awk \'{print $2}\'' % (fslDir, inFile)))
  166. pixdim2 = float(subprocess.getoutput('%sfslinfo %s | grep pixdim2 | awk \'{print $2}\'' % (fslDir, inFile)))
  167. pixdim3 = float(subprocess.getoutput('%sfslinfo %s | grep pixdim3 | awk \'{print $2}\'' % (fslDir, inFile)))
  169. # If voxel size is not 2mm isotropic, resample the data, otherwise copy the file
  170. if (pixdim1 != 2) or (pixdim2 != 2) or (pixdim3 != 2):
  171. os.system(' '.join([os.path.join(fslDir, 'flirt'),
  172. ' -ref ' + ref,
  173. ' -in ' + inFile,
  174. ' -out ' + outFile,
  175. ' -applyisoxfm 2 -interp trilinear']))
  176. else:
  177. os.system('cp ' + inFile + ' ' + outFile)
  179. # 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
  180. elif (len(affmat) == 0) and (len(warp) != 0):
  181. # Apply warp
  182. os.system(' '.join([os.path.join(fslDir, 'applywarp'),
  183. '--ref=' + ref,
  184. '--in=' + inFile,
  185. '--out=' + outFile,
  186. '--warp=' + warp,
  187. '--interp=trilinear']))
  189. # If only a affmat-file has been specified perform affine registration to MNI
  190. elif (len(affmat) != 0) and (len(warp) == 0):
  191. os.system(' '.join([os.path.join(fslDir, 'flirt'),
  192. '-ref ' + ref,
  193. '-in ' + inFile,
  194. '-out ' + outFile,
  195. '-applyxfm -init ' + affmat,
  196. '-interp trilinear']))
  198. # If both a affmat- and warp-file have been defined, apply the warping accordingly
  199. else:
  200. os.system(' '.join([os.path.join(fslDir, 'applywarp'),
  201. '--ref=' + ref,
  202. '--in=' + inFile,
  203. '--out=' + outFile,
  204. '--warp=' + warp,
  205. '--premat=' + affmat,
  206. '--interp=trilinear']))
  208. def cross_correlation(a, b):
  209. """Cross Correlations between columns of two matrices"""
  210. assert a.ndim == b.ndim == 2
  211. _, ncols_a = a.shape
  212. # nb variables in columns rather than rows hence transpose
  213. # extract just the cross terms between cols in a and cols in b
  214. return np.corrcoef(a.T, b.T)[:ncols_a, ncols_a:]
  217. def feature_time_series(melmix, mc):
  218. """ This function extracts the maximum RP correlation feature scores.
  219. It determines the maximum robust correlation of each component time-series
  220. with a model of 72 realignment parameters.
  222. Parameters
  223. ---------------------------------------------------------------------------------
  224. melmix: Full path of the melodic_mix text file
  225. mc: Full path of the text file containing the realignment parameters
  227. Returns
  228. ---------------------------------------------------------------------------------
  229. maxRPcorr: Array of the maximum RP correlation feature scores for the components
  230. of the melodic_mix file"""
  232. # Import required modules
  233. import numpy as np
  234. import random
  236. # Read melodic mix file (IC time-series), subsequently define a set of squared time-series
  237. mix = np.loadtxt(melmix)
  239. # Read motion parameter file
  240. rp6 = np.loadtxt(mc)
  241. _, nparams = rp6.shape
  243. # Determine the derivatives of the RPs (add zeros at time-point zero)
  244. rp6_der = np.vstack((np.zeros(nparams),
  245. np.diff(rp6, axis=0)
  246. ))
  248. # Create an RP-model including the RPs and its derivatives
  249. rp12 = np.hstack((rp6, rp6_der))
  251. # Add the squared RP-terms to the model
  252. # add the fw and bw shifted versions
  253. rp12_1fw = np.vstack((
  254. np.zeros(2 * nparams),
  255. rp12[:-1]
  256. ))
  257. rp12_1bw = np.vstack((
  258. rp12[1:],
  259. np.zeros(2 * nparams)
  260. ))
  261. rp_model = np.hstack((rp12, rp12_1fw, rp12_1bw))
  263. # Determine the maximum correlation between RPs and IC time-series
  264. nsplits = 1000
  265. nmixrows, nmixcols = mix.shape
  266. nrows_to_choose = int(round(0.9 * nmixrows))
  268. # Max correlations for multiple splits of the dataset (for a robust estimate)
  269. max_correls = np.empty((nsplits, nmixcols))
  270. for i in range(nsplits):
  271. # Select a random subset of 90% of the dataset rows (*without* replacement)
  272. chosen_rows = random.sample(population=range(nmixrows),
  273. k=nrows_to_choose)
  275. # Combined correlations between RP and IC time-series, squared and non squared
  276. correl_nonsquared = cross_correlation(mix[chosen_rows],
  277. rp_model[chosen_rows])
  278. correl_squared = cross_correlation(mix[chosen_rows]**2,
  279. rp_model[chosen_rows]**2)
  280. correl_both = np.hstack((correl_squared, correl_nonsquared))
  282. # Maximum absolute temporal correlation for every IC
  283. max_correls[i] = np.abs(correl_both).max(axis=1)
  285. # Feature score is the mean of the maximum correlation over all the random splits
  286. # Avoid propagating occasional nans that arise in artificial test cases
  287. return np.nanmean(max_correls, axis=0)
  290. def feature_frequency(melFTmix, TR):
  291. """ This function extracts the high-frequency content feature scores.
  292. It determines the frequency, as fraction of the Nyquist frequency,
  293. at which the higher and lower frequencies explain half
  294. of the total power between 0.01Hz and Nyquist.
  296. Parameters
  297. ---------------------------------------------------------------------------------
  298. melFTmix: Full path of the melodic_FTmix text file
  299. TR: TR (in seconds) of the fMRI data (float)
  301. Returns
  302. ---------------------------------------------------------------------------------
  303. HFC: Array of the HFC ('High-frequency content') feature scores
  304. for the components of the melodic_FTmix file"""
  306. # Import required modules
  307. import numpy as np
  309. # Determine sample frequency
  310. Fs = old_div(1, TR)
  312. # Determine Nyquist-frequency
  313. Ny = old_div(Fs, 2)
  315. # Load melodic_FTmix file
  316. FT = np.loadtxt(melFTmix)
  318. # Determine which frequencies are associated with every row in the melodic_FTmix file (assuming the rows range from 0Hz to Nyquist)
  319. f = Ny * (np.array(list(range(1, FT.shape[0] + 1)))) / (FT.shape[0])
  321. # Only include frequencies higher than 0.01Hz
  322. fincl = np.squeeze(np.array(np.where(f > 0.01)))
  323. FT = FT[fincl, :]
  324. f = f[fincl]
  326. # Set frequency range to [0-1]
  327. f_norm = old_div((f - 0.01), (Ny - 0.01))
  329. # For every IC; get the cumulative sum as a fraction of the total sum
  330. fcumsum_fract = old_div(np.cumsum(FT, axis=0), np.sum(FT, axis=0))
  332. # Determine the index of the frequency with the fractional cumulative sum closest to 0.5
  333. idx_cutoff = np.argmin(np.abs(fcumsum_fract - 0.5), axis=0)
  335. # Now get the fractions associated with those indices index, these are the final feature scores
  336. HFC = f_norm[idx_cutoff]
  338. # Return feature score
  339. return HFC
  342. def feature_spatial(fslDir, tempDir, aromaDir, melIC):
  343. """ This function extracts the spatial feature scores. For each IC it determines the fraction of the mixture modeled thresholded Z-maps respecitvely located within the CSF or at the brain edges, using predefined standardized masks.
  345. Parameters
  346. ---------------------------------------------------------------------------------
  347. fslDir: Full path of the bin-directory of FSL
  348. tempDir: Full path of a directory where temporary files can be stored (called 'temp_IC.nii.gz')
  349. aromaDir: Full path of the ICA-AROMA directory, containing the mask-files (mask_edge.nii.gz, mask_csf.nii.gz & mask_out.nii.gz)
  350. melIC: Full path of the nii.gz file containing mixture-modeled threholded (p>0.5) Z-maps, registered to the MNI152 2mm template
  352. Returns
  353. ---------------------------------------------------------------------------------
  354. edgeFract: Array of the edge fraction feature scores for the components of the melIC file
  355. csfFract: Array of the CSF fraction feature scores for the components of the melIC file"""
  357. # Import required modules
  358. import numpy as np
  359. import os
  360. import subprocess
  362. # Get the number of ICs
  363. numICs = int(subprocess.getoutput('%sfslinfo %s | grep dim4 | head -n1 | awk \'{print $2}\'' % (fslDir, melIC) ))
  365. # Loop over ICs
  366. edgeFract = np.zeros(numICs)
  367. csfFract = np.zeros(numICs)
  368. for i in range(0, numICs):
  369. # Define temporary IC-file
  370. tempIC = os.path.join(tempDir, 'temp_IC.nii.gz')
  372. # Extract IC from the merged melodic_IC_thr2MNI2mm file
  373. os.system(' '.join([os.path.join(fslDir, 'fslroi'),
  374. melIC,
  375. tempIC,
  376. str(i),
  377. '1']))
  379. # Change to absolute Z-values
  380. os.system(' '.join([os.path.join(fslDir, 'fslmaths'),
  381. tempIC,
  382. '-abs',
  383. tempIC]))
  385. # Get sum of Z-values within the total Z-map (calculate via the mean and number of non-zero voxels)
  386. totVox = int(subprocess.getoutput(' '.join([os.path.join(fslDir, 'fslstats'),
  387. tempIC,
  388. '-V | awk \'{print $1}\''])))
  390. if not (totVox == 0):
  391. totMean = float(subprocess.getoutput(' '.join([os.path.join(fslDir, 'fslstats'),
  392. tempIC,
  393. '-M'])))
  394. else:
  395. print(' - The spatial map of component ' + str(i + 1) + ' is empty. Please check!')
  396. totMean = 0
  398. totSum = totMean * totVox
  400. # Get sum of Z-values of the voxels located within the CSF (calculate via the mean and number of non-zero voxels)
  401. csfVox = int(subprocess.getoutput(' '.join([os.path.join(fslDir, 'fslstats'),
  402. tempIC,
  403. '-k mask_csf.nii.gz',
  404. '-V | awk \'{print $1}\''])))
  406. if not (csfVox == 0):
  407. csfMean = float(subprocess.getoutput(' '.join([os.path.join(fslDir, 'fslstats'),
  408. tempIC,
  409. '-k mask_csf.nii.gz',
  410. '-M'])))
  411. else:
  412. csfMean = 0
  414. csfSum = csfMean * csfVox
  416. # Get sum of Z-values of the voxels located within the Edge (calculate via the mean and number of non-zero voxels)
  417. edgeVox = int(subprocess.getoutput(' '.join([os.path.join(fslDir, 'fslstats'),
  418. tempIC,
  419. '-k mask_edge.nii.gz',
  420. '-V | awk \'{print $1}\''])))
  421. if not (edgeVox == 0):
  422. edgeMean = float(subprocess.getoutput(' '.join([os.path.join(fslDir, 'fslstats'),
  423. tempIC,
  424. '-k mask_edge.nii.gz',
  425. '-M'])))
  426. else:
  427. edgeMean = 0
  429. edgeSum = edgeMean * edgeVox
  431. # Get sum of Z-values of the voxels located outside the brain (calculate via the mean and number of non-zero voxels)
  432. outVox = int(subprocess.getoutput(' '.join([os.path.join(fslDir, 'fslstats'),
  433. tempIC,
  434. '-k mask_out.nii.gz',
  435. '-V | awk \'{print $1}\''])))
  436. if not (outVox == 0):
  437. outMean = float(subprocess.getoutput(' '.join([os.path.join(fslDir, 'fslstats'),
  438. tempIC,
  439. '-k mask_out.nii.gz',
  440. '-M'])))
  441. else:
  442. outMean = 0
  444. outSum = outMean * outVox
  446. # Determine edge and CSF fraction
  447. if not (totSum == 0):
  448. edgeFract[i] = old_div((outSum + edgeSum), (totSum - csfSum))
  449. csfFract[i] = old_div(csfSum, totSum)
  450. else:
  451. edgeFract[i] = 0
  452. csfFract[i] = 0
  454. # Remove the temporary IC-file
  455. os.remove(tempIC)
  457. # Return feature scores
  458. return edgeFract, csfFract
  461. def classification(outDir, maxRPcorr, edgeFract, HFC, csfFract):
  462. """ This function classifies a set of components into motion and
  463. non-motion components based on four features;
  464. maximum RP correlation, high-frequency content, edge-fraction and CSF-fraction
  466. Parameters
  467. ---------------------------------------------------------------------------------
  468. outDir: Full path of the output directory
  469. maxRPcorr: Array of the 'maximum RP correlation' feature scores of the components
  470. edgeFract: Array of the 'edge fraction' feature scores of the components
  471. HFC: Array of the 'high-frequency content' feature scores of the components
  472. csfFract: Array of the 'CSF fraction' feature scores of the components
  474. Return
  475. ---------------------------------------------------------------------------------
  476. motionICs Array containing the indices of the components identified as motion components
  478. Output (within the requested output directory)
  479. ---------------------------------------------------------------------------------
  480. classified_motion_ICs.txt A text file containing the indices of the components identified as motion components """
  482. # Import required modules
  483. import numpy as np
  484. import os
  486. # Classify the ICs as motion or non-motion
  488. # Define criteria needed for classification (thresholds and hyperplane-parameters)
  489. thr_csf = 0.10
  490. thr_HFC = 0.35
  491. hyp = np.array([-19.9751070082159, 9.95127547670627, 24.8333160239175])
  493. # Project edge & maxRPcorr feature scores to new 1D space
  494. x = np.array([maxRPcorr, edgeFract])
  495. proj = hyp[0] + np.dot(x.T, hyp[1:])
  497. # Classify the ICs
  498. motionICs = np.squeeze(np.array(np.where((proj > 0) + (csfFract > thr_csf) + (HFC > thr_HFC))))
  500. # Put the feature scores in a text file
  501. np.savetxt(os.path.join(outDir, 'feature_scores.txt'),
  502. np.vstack((maxRPcorr, edgeFract, HFC, csfFract)).T)
  504. # Put the indices of motion-classified ICs in a text file
  505. txt = open(os.path.join(outDir, 'classified_motion_ICs.txt'), 'w')
  506. if motionICs.size > 1: # and len(motionICs) != 0: if motionICs is not None and
  507. txt.write(','.join(['{:.0f}'.format(num) for num in (motionICs + 1)]))
  508. elif motionICs.size == 1:
  509. txt.write('{:.0f}'.format(motionICs + 1))
  510. txt.close()
  512. # Create a summary overview of the classification
  513. txt = open(os.path.join(outDir, 'classification_overview.txt'), 'w')
  514. txt.write('\t'.join(['IC',
  515. 'Motion/noise',
  516. 'maximum RP correlation',
  517. 'Edge-fraction',
  518. 'High-frequency content',
  519. 'CSF-fraction']))
  520. txt.write('\n')
  521. for i in range(0, len(csfFract)):
  522. if (proj[i] > 0) or (csfFract[i] > thr_csf) or (HFC[i] > thr_HFC):
  523. classif = "True"
  524. else:
  525. classif = "False"
  526. txt.write('\t'.join(['{:d}'.format(i + 1),
  527. classif,
  528. '{:.2f}'.format(maxRPcorr[i]),
  529. '{:.2f}'.format(edgeFract[i]),
  530. '{:.2f}'.format(HFC[i]),
  531. '{:.2f}'.format(csfFract[i])]))
  532. txt.write('\n')
  533. txt.close()
  535. return motionICs
  538. def denoising(fslDir, inFile, outDir, melmix, denType, denIdx):
  539. """ This function classifies the ICs based on the four features;
  540. maximum RP correlation, high-frequency content, edge-fraction and CSF-fraction
  542. Parameters
  543. ---------------------------------------------------------------------------------
  544. fslDir: Full path of the bin-directory of FSL
  545. inFile: Full path to the data file (nii.gz) which has to be denoised
  546. outDir: Full path of the output directory
  547. melmix: Full path of the melodic_mix text file
  548. denType: Type of requested denoising ('aggr': aggressive, 'nonaggr': non-aggressive, 'both': both aggressive and non-aggressive
  549. denIdx: Indices of the components that should be regressed out
  551. Output (within the requested output directory)
  552. ---------------------------------------------------------------------------------
  553. denoised_func_data_<denType>.nii.gz: A nii.gz file of the denoised fMRI data"""
  555. # Import required modules
  556. import os
  557. import numpy as np
  559. # Check if denoising is needed (i.e. are there components classified as motion)
  560. check = denIdx.size > 0
  562. if check == 1:
  563. # Put IC indices into a char array
  564. if denIdx.size == 1:
  565. denIdxStrJoin = "%d"%(denIdx + 1)
  566. else:
  567. denIdxStr = np.char.mod('%i', (denIdx + 1))
  568. denIdxStrJoin = ','.join(denIdxStr)
  570. # Non-aggressive denoising of the data using fsl_regfilt (partial regression), if requested
  571. if (denType == 'nonaggr') or (denType == 'both'):
  572. os.system(' '.join([os.path.join(fslDir, 'fsl_regfilt'),
  573. '--in=' + inFile,
  574. '--design=' + melmix,
  575. '--filter="' + denIdxStrJoin + '"',
  576. '--out=' + os.path.join(outDir, 'denoised_func_data_nonaggr.nii.gz')]))
  578. # Aggressive denoising of the data using fsl_regfilt (full regression)
  579. if (denType == 'aggr') or (denType == 'both'):
  580. os.system(' '.join([os.path.join(fslDir, 'fsl_regfilt'),
  581. '--in=' + inFile,
  582. '--design=' + melmix,
  583. '--filter="' + denIdxStrJoin + '"',
  584. '--out=' + os.path.join(outDir, 'denoised_func_data_aggr.nii.gz'),
  585. '-a']))
  586. else:
  587. print(" - None of the components were classified as motion, so no denoising is applied (a symbolic link to the input file will be created).")
  588. if (denType == 'nonaggr') or (denType == 'both'):
  589. os.symlink(inFile, os.path.join(outDir, 'denoised_func_data_nonaggr.nii.gz'))
  590. if (denType == 'aggr') or (denType == 'both'):
  591. os.symlink(inFile, os.path.join(outDir, 'denoised_func_data_aggr.nii.gz'))

ICA_AROMA_functions_preICA.py at commit c16de31, no license · at the source

Overview

Authors: Zhishun Wang1,2, Feng Liu1,2, Rachel Marsh1,2, Gaurav H. Patel1,2, Jack Grinband1,2
ORCID iDs: Zhishun Wang, Jack Grinband
  1. The Department of Psychiatry, Vagelos College of Physicians and Surgeons, Columbia University, New York, NY, United States
  2. New York State Psychiatric Institute, New York, NY, United States
Institutions: New York State Psychiatric Institute; Columbia University (United States)
Journal: Imaging neuroscience (Cambridge, Mass.), volume 4, article IMAG.a.1198
Dates: received 7 October 2025; accepted 8 March 2026; published online 6 April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1162/imag.a.1198 · PMID 41952841 · PMCID PMC13055012 · OpenAlex W7136943186
Open access: diamond, a free copy (OpenAlex)
Status: code verified
Categories: fMRI (modality), methods / tools (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Connectivity, Graphs, fMRI & imaging
Keywords: multi-echo functional MRI (ME-fMRI), functional connectome, fMRIPrep, tedana, independent component analysis (ICA), Nipype
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Citations: not cited yet (Europe PMC); 73 references in the paper

Abstract

Multi-echo fMRI has emerged as a powerful strategy to mitigate head motion-related noise and minimize susceptibility-related signal loss in BOLD data. Multi-echo independent component analysis (ME-ICA) effectively distinguishes between BOLD-related (TE-dependent) signals and non-BOLD (TE-independent) noise, yielding substantial enhancements in performance compared to traditional echo-combination methods. We introduce a novel ICA-based denoising step, preICA, applied to raw multi-echo data before optimal T2*-weighted echo combination. This approach, combined with ME-ICA, yields substantial gains in data denoising. Our results show that preICA significantly enhances the efficacy of optimal echo combination and ME-ICA to reduce noise. To facilitate the reliable processing of multi-echo fMRI data, we integrated preICA and ME-ICA into fMRIPrep, resulting in the creation of a robust multi-echo processing pipeline, called MEPrep, offering flexibility in preprocessing options (with or without preICA and/or ME-ICA) beyond the echo combination approach offered by fMRIPrep. We validated MEPrep on an open resting-state multi-echo fMRI dataset, demonstrating that incorporating the preICA step leads to statistically significant improvements in denoising efficacy, as evidenced by (1) enhanced T2* exponential model fitting accuracy; (2) reduced motion-related BOLD fluctuations; (3) increased temporal signal-to-noise ratio; (4) improved spatial and temporal reliability of functional connectivity; and (5) increased Shannon entropy. MEPrep outperforms existing pipelines by synergistically integrating preICA and ME-ICA, achieving superior noise suppression while preserving the neurobiological complexity of denoised BOLD signals. By automating multi-echo preprocessing within a robust pipeline, MEPrep provides a scalable solution for high-quality multi-echo fMRI data preprocessing. The pipeline is openly available, ensuring reproducibility and accessibility for the neuroimaging community.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 6 matches between paragraphs and lines of code.

zswang-brainimaging/MEPrep

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: c16de3143eff7ed5e475a810adc5ca7a3a541ff5, 6 April 2026
Languages: Python (97), Shell (1)
Size: 651 files, 98 scripts
Software Heritage: not archived
Found in: “Data and Code Availability”
Holds: README, environment (MEPrep_Source/Dockerfile, MEPrep_Source/env.yml, MEPrep_Source/pyproject.toml, MEPrep_Source/requirements.txt, MEPrep_Source/MEPrep_ICA-AROMA/Dockerfile, MEPrep_Source/MEPrep_ICA-AROMA/requirements.txt), tests
Not found: license file, CITATION.cff, continuous integration, documentation
Tools: Nipype (50 files), NumPy (29 files), NiBabel (23 files), fMRIPrep (21 files), PyBIDS (19 files), FSL (14 files), TemplateFlow (9 files), FreeSurfer (8 files), pandas (4 files), SciPy (3 files), Matplotlib (2 files), seaborn (2 files), tedana (2 files), AFNI (1 file), BIDS Validator (1 file), h5py (1 file), Nilearn (1 file), scikit-image (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
99 files
At the source: github.com/zswang-brainimaging/MEPrep

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 98 scripts, each with its path and the digest of its content;
  • 6 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data and Code Availability

Our MEPrep pipeline is publicly available on Docker Hub. Users can download and install MEPrep on a Docker-enabled system by running a command under Docker: “docker pull zswang2020/meprep_final:latest”. In addition to the Docker image, the full MEPrep source code and documentation have been made available in the GitHub repository at https://github.com/zswang-brainimaging/MEPrep to facilitate extension, community contributions, and integration into other neuroimaging workflows. A quick user guide on how to install and use MEPrep is available in the Supplementary Methods (https://doi.org/10.1162/IMAG.a.1198#supplementary-data). The readers can test this pipeline using the same dataset as we used, the Siemens dataset (Multi-echo Cambridge), which can be downloaded from https://openfmri.org/dataset/ds000258/.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 5 authors, 6 keywords, 72 references.

Cite

This paper

Wang, Z., Liu, F., Marsh, R., Patel, G. H., & Grinband, J. (2026). MEPrep: A robust pipeline for multi-echo fMRI denoising and preprocessing. Imaging neuroscience (Cambridge, Mass.), 4, IMAG.a.1198. https://doi.org/10.1162/imag.a.1198

BibTeX

@article{wang2026meprep,
author = {Wang, Zhishun and Liu, Feng and Marsh, Rachel and Patel, Gaurav H. and Grinband, Jack},
title = {{MEPrep: A robust pipeline for multi-echo fMRI denoising and preprocessing}},
journal = {Imaging neuroscience (Cambridge, Mass.)},
year = {2026},
month = apr,
volume = {4},
pages = {IMAG.a.1198},
publisher = {MIT Press},
issn = {2837-6056},
doi = {10.1162/imag.a.1198},
url = {https://doi.org/10.1162/imag.a.1198},
pmid = {41952841},
pmcid = {PMC13055012}
}

RIS

TY - JOUR
AU - Wang, Zhishun
AU - Liu, Feng
AU - Marsh, Rachel
AU - Patel, Gaurav H.
AU - Grinband, Jack
TI - MEPrep: A robust pipeline for multi-echo fMRI denoising and preprocessing
T2 - Imaging neuroscience (Cambridge, Mass.)
J2 - Imaging Neurosci (Camb)
PY - 2026
DA - 2026/04/06
VL - 4
SP - IMAG.a.1198
SN - 2837-6056
PB - MIT Press
DO - 10.1162/imag.a.1198
UR - https://doi.org/10.1162/imag.a.1198
LA - en
ER -

CSL-JSON

{
"id": "10.1162/imag.a.1198",
"type": "article-journal",
"title": "MEPrep: A robust pipeline for multi-echo fMRI denoising and preprocessing",
"container-title": "Imaging neuroscience (Cambridge, Mass.)",
"author": [
{
"family": "Wang",
"given": "Zhishun"
},
{
"family": "Liu",
"given": "Feng"
},
{
"family": "Marsh",
"given": "Rachel"
},
{
"family": "Patel",
"given": "Gaurav H."
},
{
"family": "Grinband",
"given": "Jack"
}
],
"container-title-short": "Imaging Neurosci (Camb)",
"volume": "4",
"page": "IMAG.a.1198",
"DOI": "10.1162/imag.a.1198",
"PMID": "41952841",
"PMCID": "PMC13055012",
"ISSN": "2837-6056",
"publisher": "MIT Press",
"URL": "https://doi.org/10.1162/imag.a.1198",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
6
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1162/netn.a.547 [code]
An evaluation of the efficacy of single-echo and multi-echo fMRI denoising strategies.
Journal: Network neuroscience (Cambridge, Mass.)
In common: tedana, fMRIPrep, Nipype, 9 other tools, fMRI, 21 references
[2] doi:10.1002/hbm.70621 [code]
Optimising 7T-fMRI for Imaging Regions of Magnetic Susceptibility.
Journal: Human brain mapping
In common: tedana, fMRIPrep, AFNI, 1 other tool, fMRI, 17 references
[3] doi:10.1162/imag.a.1262 [code]
Frame-wise multi-echo distortion correction for superior functional MRI.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: PyBIDS, FreeSurfer, FSL, 8 other tools, fMRI, 8 references
[4] doi:10.1038/s41398-026-04157-5 [code]
Association of glymphatic function with 40-Hz neural oscillations, systemic metabolic markers, and cognitive performance in healthy aging adults: An EEG and MRI study.
Journal: Translational psychiatry
In common: TemplateFlow, fMRIPrep, PyBIDS, 11 other tools, fMRI, 1 reference
[5] doi:10.1162/imag.a.1362 [code]
Human fMRI at 11.7T: Assessing feasibility, stability, and reliability on the Iseult scanner.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: fMRIPrep, Nipype, FreeSurfer, 8 other tools, fMRI, methods / tools, 6 references
[6] doi:10.1162/imag.a.1347 [code]
Neural and behavioural correlates of theory of mind reasoning in five-year-old children born preterm.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: tedana, fMRIPrep, PyBIDS, 9 other tools, fMRI, 3 references
[7] doi:10.1162/imag.a.1252 [code]
Does the brain's E:I balance really shape long-range temporal correlations? Lessons learned from 3T MRI.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: tedana, fMRIPrep, AFNI, 8 other tools, fMRI, 6 references
[8] doi:10.1016/j.neuron.2026.04.011 [code]
Precision fMRI reveals densely interdigitated network patches with conserved motifs in the lateral prefrontal cortex.
Journal: Neuron
In common: tedana, AFNI, FreeSurfer, 6 other tools, fMRI, 8 references
[9] doi:10.1038/s41467-026-71151-2 [code]
Common and distinct neural correlates of social interaction processing and theory of mind in narratives.
Journal: Nature communications
In common: fMRIPrep, Nipype, AFNI, 9 other tools, 5 references
[10] doi:10.1016/j.dcn.2026.101765 [code]
Fusiform face area development correlates with development in higher-order social brain regions.
Journal: Developmental cognitive neuroscience
In common: tedana, fMRIPrep, PyBIDS, 9 other tools, fMRI, 2 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.