From 391c6e1117981dfb921b823f5bb77bf49e5ae9ed Mon Sep 17 00:00:00 2001 From: "Manjari Narayan (SU)" Date: Thu, 29 Aug 2019 17:32:46 -0700 Subject: [PATCH] Saving modified denoiser module --- .../02-denoiser/denoiser/.circleci/config.yml | 94 +++ .../.circleci/test_data_downloader.sh | 8 + src/data/02-denoiser/denoiser/.gitignore | 102 +++ src/data/02-denoiser/denoiser/LICENSE | 201 ++++++ src/data/02-denoiser/denoiser/README.rst | 17 + .../02-denoiser/denoiser/report_template.html | 85 +++ .../02-denoiser/denoiser/requirements.txt | 13 + src/data/02-denoiser/denoiser/run.py | 1 + src/data/02-denoiser/denoiser/run_denoise.py | 678 ++++++++++++++++++ 9 files changed, 1199 insertions(+) create mode 100644 src/data/02-denoiser/denoiser/.circleci/config.yml create mode 100644 src/data/02-denoiser/denoiser/.circleci/test_data_downloader.sh create mode 100644 src/data/02-denoiser/denoiser/.gitignore create mode 100644 src/data/02-denoiser/denoiser/LICENSE create mode 100644 src/data/02-denoiser/denoiser/README.rst create mode 100644 src/data/02-denoiser/denoiser/report_template.html create mode 100644 src/data/02-denoiser/denoiser/requirements.txt create mode 100644 src/data/02-denoiser/denoiser/run.py create mode 100755 src/data/02-denoiser/denoiser/run_denoise.py diff --git a/src/data/02-denoiser/denoiser/.circleci/config.yml b/src/data/02-denoiser/denoiser/.circleci/config.yml new file mode 100644 index 0000000..ae3e925 --- /dev/null +++ b/src/data/02-denoiser/denoiser/.circleci/config.yml @@ -0,0 +1,94 @@ +# Python CircleCI 2.0 configuration file +# +# Check https://circleci.com/docs/2.0/language-python/ for more details +# +version: 2 +jobs: + build: + docker: + # specify the version you desire here + # use `-browsers` prefix for selenium tests, e.g. `3.6.1-browsers` + - image: circleci/python:3.6.1 + + # Specify service dependencies here if necessary + # CircleCI maintains a library of pre-built images + # documented at https://circleci.com/docs/2.0/circleci-images/ + # - image: circleci/postgres:9.4 + + working_directory: ~/repo + + steps: + - checkout + + # Download and cache dependencies + - restore_cache: + keys: + - v1-dependencies-{{ checksum "requirements.txt" }} + # fallback to using the latest cache if no exact match is found + - v1-dependencies- + + + - run: + name: Install dependencies + command: | + python3 -m venv venv + . venv/bin/activate + pip install -r requirements.txt + + - save_cache: + paths: + - ./venv + key: v1-dependencies-{{ checksum "requirements.txt" }} + + - restore_cache: + keys: + - test-data + + - run: + name: Download test data + command: | + bash .circleci/test_data_downloader.sh + + - save_cache: + paths: + - ./test_data + key: test-data + + # run tests! + - run: + name: Smoke test - print help + command: | + . venv/bin/activate + ./run_denoise.py --help + + - run: + name: Smoke test - basic denoising + command: | + . venv/bin/activate + mkdir -p test_outputs/basic + ./run_denoise.py test_data/sub-01_task-rhymejudgment_bold_space-MNI152NLin2009cAsym_preproc.nii.gz test_data/sub-01_task-rhymejudgment_bold_confounds.tsv test_outputs/basic + + - run: + name: Smoke test - select subsest of columns + command: | + . venv/bin/activate + mkdir -p test_outputs/subset + ./run_denoise.py test_data/sub-01_task-rhymejudgment_bold_space-MNI152NLin2009cAsym_preproc.nii.gz test_data/sub-01_task-rhymejudgment_bold_confounds.tsv test_outputs/subset --col_names FramewiseDisplacement GlobalSignal + + - run: + name: Smoke test - high pass + command: | + . venv/bin/activate + mkdir -p test_outputs/hp + ./run_denoise.py test_data/sub-01_task-rhymejudgment_bold_space-MNI152NLin2009cAsym_preproc.nii.gz test_data/sub-01_task-rhymejudgment_bold_confounds.tsv test_outputs/hp --hp_filter 0.009 + + - run: + name: Smoke test - low pass + command: | + . venv/bin/activate + mkdir -p test_outputs/lp + ./run_denoise.py test_data/sub-01_task-rhymejudgment_bold_space-MNI152NLin2009cAsym_preproc.nii.gz test_data/sub-01_task-rhymejudgment_bold_confounds.tsv test_outputs/lp --lp_filter 0.1 + + - store_artifacts: + path: test_outputs + destination: test_outputs diff --git a/src/data/02-denoiser/denoiser/.circleci/test_data_downloader.sh b/src/data/02-denoiser/denoiser/.circleci/test_data_downloader.sh new file mode 100644 index 0000000..047b461 --- /dev/null +++ b/src/data/02-denoiser/denoiser/.circleci/test_data_downloader.sh @@ -0,0 +1,8 @@ +#!/usr/bin/env bash +set -x +set -e +if [ ! -e test_data ]; then + mkdir test_data + wget http://openneuro.outputs.s3.amazonaws.com/091b91a257f6b517790f2fb82a784c8e/3e22d7c1-c7e8-4cc7-bfdd-c89dd6681982/fmriprep/sub-01/func/sub-01_task-rhymejudgment_bold_space-MNI152NLin2009cAsym_preproc.nii.gz -O test_data/sub-01_task-rhymejudgment_bold_space-MNI152NLin2009cAsym_preproc.nii.gz + wget http://openneuro.outputs.s3.amazonaws.com/091b91a257f6b517790f2fb82a784c8e/3e22d7c1-c7e8-4cc7-bfdd-c89dd6681982/fmriprep/sub-01/func/sub-01_task-rhymejudgment_bold_confounds.tsv -O test_data/sub-01_task-rhymejudgment_bold_confounds.tsv +fi diff --git a/src/data/02-denoiser/denoiser/.gitignore b/src/data/02-denoiser/denoiser/.gitignore new file mode 100644 index 0000000..c49da9d --- /dev/null +++ b/src/data/02-denoiser/denoiser/.gitignore @@ -0,0 +1,102 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +env/ +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +*.egg-info/ +.installed.cfg +*.egg + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +.hypothesis/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# pyenv +.python-version + +# celery beat schedule file +celerybeat-schedule + +# SageMath parsed files +*.sage.py + +# dotenv +.env + +# virtualenv +.venv +venv/ +ENV/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.idea/ \ No newline at end of file diff --git a/src/data/02-denoiser/denoiser/LICENSE b/src/data/02-denoiser/denoiser/LICENSE new file mode 100644 index 0000000..261eeb9 --- /dev/null +++ b/src/data/02-denoiser/denoiser/LICENSE @@ -0,0 +1,201 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/src/data/02-denoiser/denoiser/README.rst b/src/data/02-denoiser/denoiser/README.rst new file mode 100644 index 0000000..c0695b5 --- /dev/null +++ b/src/data/02-denoiser/denoiser/README.rst @@ -0,0 +1,17 @@ +Denoiser: A nuisance regression tool for fMRI BOLD data +======================================================= + +``Denoiser`` is a tool for removing sources of noise from and performing temporal filtering +of functional MRI data. It also provides visualization of the content of nuisance signals +(including *motion* information, if provided), allowing the user to get a quick sense of +data quality before and after noise removal. + +Denoiser acts on 4D fMRI data (takes either a nifti file path or an already loaded nibabel +object as input). Nuisance signal removal and temporal filtering are performed on a +voxel-wise level, and a 'cleaned' 4D Nifti file (``_NR`` file)/ nibabel object are created +as outputs. + +The specific noise signals to be removed are specified by the user (contained in a tsv file). +This tool should be used only **after** BOLD data are *minimally preprocessed* (for example, +after preprocessing the data using ``fmriprep``, which creates a .tsv file containing nuisance signals). + diff --git a/src/data/02-denoiser/denoiser/report_template.html b/src/data/02-denoiser/denoiser/report_template.html new file mode 100644 index 0000000..679eec9 --- /dev/null +++ b/src/data/02-denoiser/denoiser/report_template.html @@ -0,0 +1,85 @@ + + + + + + + +

Denoising Report for:

+

{{ img_file }}

+ +
+ +

Output 'Cleaned' file:

+

{{ save_img_file }}

+ +
+

{{ Ntrs }} TRs detected

+ +

Nuisance signal file:

+

{{ tsv_file }}

+ +

Signals removed:

+ + +{% if hp_filter %} + +

High-pass filter applied [frequencies < {{ hp_filter }} removed]

+{% else %} +

***NO additional high-pass filter applied***

+

Note that if Nuisance regressors contain low-frequency DRIFT and if low-frequency cosines were not included as Nuisance regressors - drift-like content will be present in 'cleaned' output data!!

+{% endif %} + +{% if lp_filter %} +

Low-pass filter applied [frequencies > {{ lp_filter }} removed]

+{% else %} +

NO Low-pass filter applied

+{% endif %} + +
+ +

Nuisance Regressor Information

+Nuisance Design Matrix + +Nuisance Corr Matrix + +
+ +

Motion Information

+

FD timeseries plot

+ +
+ +

Data (voxel time series) before and after Nuisance Regression

+

carpet plots

+ +
+ +

Nuisance Regressor Stat Maps

+ +

Total Nuisance R2 Map

+

R squared map

+
+ +{% for item in file_tstat %} + +

{{ item.name }} Tstat Map

+

{{ item.name }} Map

+
+{% endfor %} + + + + \ No newline at end of file diff --git a/src/data/02-denoiser/denoiser/requirements.txt b/src/data/02-denoiser/denoiser/requirements.txt new file mode 100644 index 0000000..4253029 --- /dev/null +++ b/src/data/02-denoiser/denoiser/requirements.txt @@ -0,0 +1,13 @@ +git+https://github.com/chrisfilo/nilearn.git@enh/carpet_plot +nistats +nibabel +numpy +pandas +matplotlib +scipy +sklearn +patsy +seaborn +jinja2 +git+https://github.com/poldracklab/niworkflows +nitime \ No newline at end of file diff --git a/src/data/02-denoiser/denoiser/run.py b/src/data/02-denoiser/denoiser/run.py new file mode 100644 index 0000000..89cceda --- /dev/null +++ b/src/data/02-denoiser/denoiser/run.py @@ -0,0 +1 @@ +print("I'm running!") \ No newline at end of file diff --git a/src/data/02-denoiser/denoiser/run_denoise.py b/src/data/02-denoiser/denoiser/run_denoise.py new file mode 100755 index 0000000..0e6950c --- /dev/null +++ b/src/data/02-denoiser/denoiser/run_denoise.py @@ -0,0 +1,678 @@ +#!/usr/bin/env python3 +from nistats import regression +from nistats import reporting +from nistats.design_matrix import make_first_level_design_matrix +import nibabel as nb +import numpy as np +import os, pandas, sys, pdb, argparse, copy, scipy, jinja2 +from os.path import join as pjoin +from nilearn import plotting +from nilearn import _utils +from nilearn.signal import butterworth, clean +from nilearn.input_data import NiftiMasker +from functools import reduce +from collections import deque, OrderedDict + +from matplotlib import pylab as plt +from matplotlib import cm +from matplotlib import axes +import seaborn as sns +from nilearn._utils.niimg import load_niimg, _safe_get_data +from nipype.algorithms import confounds as nac + +parser = argparse.ArgumentParser(description='Function for performing nuisance regression. Saves resulting output ' + 'nifti file, information about nuisance regressors and motion (html ' + 'report), and outputs nibabel object containing clean data') +parser.add_argument('img_file', help='4d nifti img: file path or nibabel object loaded into memory') +parser.add_argument('tsv_file', help='tsv file containing nuisance regressors to be removed') +parser.add_argument('out_path', help='output directory for saving new data file') +parser.add_argument('--col_names', + help='which columns of TSV file to include as nuisance regressors. defaults to ALL columns.', + nargs="+") +parser.add_argument('--hp_filter', help='frequency cut-off for high pass filter (removing low frequencies). Recommend ' + '.009 Hz') +parser.add_argument('--lp_filter', help='frequency cut-off for low pass filter (removing high frequencies). Recommend ' + '.1 Hz for non-task data') +parser.add_argument('--out_figure_path', + help='output directory for saving figures. Defaults to location of out_path + _figures') + +args = parser.parse_args() + +img_file = args.img_file +tsv_file = args.tsv_file +out_path = args.out_path +col_names = args.col_names +hp_filter = args.hp_filter +lp_filter = args.lp_filter +out_figure_path = args.out_figure_path + +def plot_carpet(img, mask_img=None, detrend=True, output_file=None, + figure=None, axes=None, title=None): + """Plot an image representation of voxel intensities across time also know + as the "carpet plot" or "Power plot". See Jonathan Power Neuroimage + 2017 Jul 1; 154:150-158. + Parameters + ---------- + img : Niimg-like object + See http://nilearn.github.io/manipulating_images/input_output.html + 4D input image + mask_img : Niimg-like object, optional + See http://nilearn.github.io/manipulating_images/input_output.html + Limit plotted voxels to those inside the provided mask. If not + specified a new mask will be derived from data. + detrend : boolean, optional + Detrend and standardize the data prior to plotting. + output_file : string, or None, optional + The name of an image file to export the plot to. Valid extensions + are .png, .pdf, .svg. If output_file is not None, the plot + is saved to a file, and the display is closed. + figure : matplotlib figure, optional + Matplotlib figure used. If None is given, a + new figure is created. + axes : matplotlib axes, optional + The axes used to display the plot. If None, the complete + figure is used. + title : string, optional + The title displayed on the figure. + """ + img_nii = _utils.check_niimg_4d(img, dtype='auto') + img_data = _safe_get_data(img_nii, ensure_finite=True) + + # Define TR and number of frames + tr = img_nii.header.get_zooms()[-1] + ntsteps = img_nii.shape[-1] + + if not mask_img: + nifti_masker = NiftiMasker(mask_strategy='epi', standardize=False) + nifti_masker.fit(img) + mask_data = nifti_masker.mask_img_.get_data().astype(bool) + else: + mask_nii = _utils.check_niimg_3d(mask_img, dtype='auto') + mask_data = _safe_get_data(mask_nii, ensure_finite=True) + + data = img_data[mask_data > 0].reshape(-1, ntsteps) + # Detrend data + if detrend: + data = clean(data.T, t_r=tr).T + + if not figure: + if not axes: + figure = plt.figure() + else: + figure = axes.figure + + if not axes: + axes = figure.add_subplot(1, 1, 1) + else: + assert axes.figure is figure, ("The axes passed are not " + "in the figure") + + # Avoid segmentation faults for long acquisitions by decimating the input + # data + long_cutoff = 800 + if data.shape[1] > long_cutoff: + data = data[:, ::2] + else: + data = data[:, :] + + axes.imshow(data, interpolation='nearest', + aspect='auto', cmap='gray', vmin=-2, vmax=2) + + axes.grid(False) + axes.set_yticks([]) + axes.set_yticklabels([]) + + # Set 10 frame markers in X axis + interval = max( + (int(data.shape[-1] + 1) // 10, int(data.shape[-1] + 1) // 5, 1)) + xticks = list(range(0, data.shape[-1])[::interval]) + axes.set_xticks(xticks) + + axes.set_xlabel('time (s)') + axes.set_ylabel('voxels') + if title: + axes.set_title(title) + labels = tr * (np.array(xticks)) + if data.shape[1] > long_cutoff: + labels *= 2 + axes.set_xticklabels(['%.02f' % t for t in labels.tolist()]) + + # Remove and redefine spines + for side in ["top", "right"]: + # Toggle the spine objects + axes.spines[side].set_color('none') + axes.spines[side].set_visible(False) + + axes.yaxis.set_ticks_position('left') + axes.xaxis.set_ticks_position('bottom') + axes.spines["bottom"].set_position(('outward', 20)) + axes.spines["left"].set_position(('outward', 20)) + + if output_file is not None: + figure.savefig(output_file) + figure.close() + figure = None + + return figure + +def temporal_derivatives(order, variables, data): + """ + Compute temporal derivative terms by the method of backwards differences. + Parameters + ---------- + order: range or list(int) + A list of temporal derivative terms to include. For instance, [1, 2] + indicates that the first and second derivative terms should be added. + To retain the original terms, 0 *must* be included in the list. + variables: list(str) + List of variables for which temporal derivative terms should be + computed. + data: pandas DataFrame object + Table of values of all observations of all variables. + Returns + ------- + variables_deriv: list + A list of variables to include in the final data frame after adding + the specified derivative terms. + data_deriv: pandas DataFrame object + Table of values of all observations of all variables, including any + specified derivative terms. + Credit: niworkflow. + """ + variables_deriv = OrderedDict() + data_deriv = OrderedDict() + if 0 in order: + data_deriv[0] = data[variables] + variables_deriv[0] = variables + order = set(order) - set([0]) + for o in order: + variables_deriv[o] = ['{}_derivative{}'.format(v, o) + for v in variables] + data_deriv[o] = np.tile(np.nan, data[variables].shape) + data_deriv[o][o:, :] = np.diff(data[variables], n=o, axis=0) + variables_deriv = reduce((lambda x, y: x + y), variables_deriv.values()) + data_deriv = pandas.DataFrame(columns=variables_deriv, + data=np.concatenate([*data_deriv.values()], + axis=1)) + + return (variables_deriv, data_deriv) + + +def exponential_terms(order, variables, data): + """ + Compute exponential expansions. + Parameters + ---------- + order: range or list(int) + A list of exponential terms to include. For instance, [1, 2] + indicates that the first and second exponential terms should be added. + To retain the original terms, 1 *must* be included in the list. + variables: list(str) + List of variables for which exponential terms should be computed. + data: pandas DataFrame object + Table of values of all observations of all variables. + Returns + ------- + variables_exp: list + A list of variables to include in the final data frame after adding + the specified exponential terms. + data_exp: pandas DataFrame object + Table of values of all observations of all variables, including any + specified exponential terms. + Credit: niworkflow + """ + variables_exp = OrderedDict() + data_exp = OrderedDict() + if 1 in order: + data_exp[1] = data[variables] + variables_exp[1] = variables + order = set(order) - set([1]) + for o in order: + variables_exp[o] = ['{}_power{}'.format(v, o) for v in variables] + data_exp[o] = data[variables]**o + variables_exp = reduce((lambda x, y: x + y), variables_exp.values()) + data_exp = pandas.DataFrame(columns=variables_exp, + data=np.concatenate([*data_exp.values()], axis=1)) + return (variables_exp, data_exp) + + +def denoise(img_file, tsv_file, out_path, col_names=False, hp_filter=False, lp_filter=False, out_figure_path=False): + nii_ext = '.nii.gz' + FD_thr = [.5] + sc_range = np.arange(-1, 3) + constant = 'constant' + + # read in files + img = load_niimg(img_file) + # get file info + img_name = os.path.basename(img.get_filename()) + file_base = img_name[0:img_name.find('.')] + save_img_file = pjoin(out_path, file_base + \ + '_NR' + nii_ext) + data = img.get_data() + df_orig = pandas.read_csv(tsv_file, '\t', na_values='n/a') + df = copy.deepcopy(df_orig) + Ntrs = df.as_matrix().shape[0] + print('# of TRs: ' + str(Ntrs)) + assert (Ntrs == data.shape[len(data.shape) - 1]) + + # select columns to use as nuisance regressors + if col_names: + df = df[col_names] + str_append = ' [SELECTED regressors in CSV]' + else: + col_names = df.columns.tolist() + str_append = ' [ALL regressors in CSV]' + + # trans_x trans_y trans_z rot_x rot_y rot_z + motion_names = ['trans_x','trans_y','trans_z','rot_x','rot_y','rot_z'] + motion = df[motion_names] + ## square movement variables and add (Friston 12P) + ## Edit: Removing below and using niworkflow functions instead. + # motion_squared=motion**2 + # motion_squared.columns=['%s**2'%i for i in motion_squared.columns] + # df = pandas.concat([df, motion_squared], axis=1) + ## Create Friston 24P + motion_dd_names,motion_dd = temporal_derivatives( + [0,1], + motion_names, + motion) + motion_ddsq_names,motion_ddsq = exponential_terms( + [2], + motion_dd_names, + motion_dd) + motion_dd = motion_dd[motion_dd_names[6:12]] + motion_ddsq = pandas.concat([motion_dd, motion_ddsq], axis=1) + df = pandas.concat([df, motion_ddsq], axis=1) + + # fill in missing nuisance values with mean for that variable + for col in df.columns: + if sum(df[col].isnull()) > 0: + print('Filling in ' + str(sum(df[col].isnull())) + ' NaN value for ' + col) + df[col] = df[col].fillna(np.mean(df[col])) + print('# of Confound Regressors: ' + str(len(df.columns)) + str_append) + + # implement HP filter in regression + TR = img.header.get_zooms()[-1] + frame_times = np.arange(Ntrs) * TR + if hp_filter: + hp_filter = float(hp_filter) + assert (hp_filter > 0) + period_cutoff = 1. / hp_filter + df = make_first_level_design_matrix(frame_times, period_cut=period_cutoff, add_regs=df.as_matrix(), + add_reg_names=df.columns.tolist()) + # fn adds intercept into dm + + hp_cols = [col for col in df.columns if 'drift' in col] + print('# of High-pass Filter Regressors: ' + str(len(hp_cols))) + else: + # add in intercept column into data frame + df[constant] = 1 + print('No High-pass Filter Applied') + + dm = df.as_matrix() + + # prep data + data = np.reshape(data, (-1, Ntrs)) + data_mean = np.mean(data, axis=1) + Nvox = len(data_mean) + + # setup and run regression + model = regression.OLSModel(dm) + results = model.fit(data.T) + if not hp_filter: + results_orig_resid = copy.deepcopy(results.resid) # save for rsquared computation + + # apply low-pass filter + if lp_filter: + # input to butterworth fn is time x voxels + low_pass = float(lp_filter) + Fs = 1. / TR + if low_pass >= Fs / 2: + raise ValueError('Low pass filter cutoff if too close to the Nyquist frequency (%s)' % (Fs / 2)) + + temp_img_file = pjoin(out_path, file_base + \ + '_temp' + nii_ext) + temp_img = nb.Nifti1Image(np.reshape(results.resid.T + np.reshape(data_mean, (Nvox, 1)), img.shape).astype('float32'), + img.affine, header=img.header) + temp_img.to_filename(temp_img_file) + results.resid = butterworth(results.resid, sampling_rate=Fs, low_pass=low_pass, high_pass=None) + print('Low-pass Filter Applied: < ' + str(low_pass) + ' Hz') + + # add mean back into data + clean_data = results.resid.T + np.reshape(data_mean, (Nvox, 1)) # add mean back into residuals + + # save out new data file + print('Saving output file...') + clean_data = np.reshape(clean_data, img.shape).astype('float32') + new_img = nb.Nifti1Image(clean_data, img.affine, header=img.header) + new_img.to_filename(save_img_file) + + ######### generate Rsquared map for confounds only + if hp_filter: + # first remove low-frequency information from data + hp_cols.append(constant) + model_first = regression.OLSModel(df[hp_cols].as_matrix()) + results_first = model_first.fit(data.T) + results_first_resid = copy.deepcopy(results_first.resid) + del results_first, model_first + + # compute sst - borrowed from matlab + sst = np.square(np.linalg.norm(results_first_resid - + np.mean(results_first_resid, axis=0), axis=0)) + + # now regress out 'true' confounds to estimate their Rsquared + nr_cols = [col for col in df.columns if 'drift' not in col] + model_second = regression.OLSModel(df[nr_cols].as_matrix()) + results_second = model_second.fit(results_first_resid) + + # compute sse - borrowed from matlab + sse = np.square(np.linalg.norm(results_second.resid, axis=0)) + + del results_second, model_second, results_first_resid + + elif not hp_filter: + # compute sst - borrowed from matlab + sst = np.square(np.linalg.norm(data.T - + np.mean(data.T, axis=0), axis=0)) + + # compute sse - borrowed from matlab + sse = np.square(np.linalg.norm(results_orig_resid, axis=0)) + + del results_orig_resid + + # compute rsquared of nuisance regressors + zero_idx = scipy.logical_and(sst == 0, sse == 0) + sse[zero_idx] = 1 + sst[zero_idx] = 1 # would be NaNs - become rsquared = 0 + rsquare = 1 - np.true_divide(sse, sst) + rsquare[np.isnan(rsquare)] = 0 + + ######### Visualizing DM & outputs + fontsize = 12 + fontsize_title = 14 + def_img_size = 8 + + if not out_figure_path: + out_figure_path = save_img_file[0:save_img_file.find('.')] + '_figures' + + if not os.path.isdir(out_figure_path): + os.mkdir(out_figure_path) + png_append = '_' + img_name[0:img_name.find('.')] + '.png' + print('Output directory: ' + out_figure_path) + + # DM corr matrix + cm = df[df.columns[0:-1]].corr() + curr_sz = copy.deepcopy(def_img_size) + if cm.shape[0] > def_img_size: + curr_sz = curr_sz + ((cm.shape[0] - curr_sz) * .3) + mtx_scale = curr_sz * 100 + + mask = np.zeros_like(cm, dtype=np.bool) + mask[np.triu_indices_from(mask)] = True + + fig, ax = plt.subplots(figsize=(curr_sz, curr_sz)) + cmap = sns.diverging_palette(220, 10, as_cmap=True) + sns.heatmap(cm, mask=mask, cmap=cmap, center=0, vmax=cm[cm < 1].max().max(), vmin=cm[cm < 1].min().min(), + square=True, linewidths=.5, cbar_kws={"shrink": .6}) + ax.set_xticklabels(ax.get_xticklabels(), rotation=60, ha='right', fontsize=fontsize) + ax.set_yticklabels(cm.columns.tolist(), rotation=-30, va='bottom', fontsize=fontsize) + ax.set_title('Nuisance Corr. Matrix', fontsize=fontsize_title) + plt.tight_layout() + file_corr_matrix = 'Corr_matrix_regressors' + png_append + fig.savefig(pjoin(out_figure_path, file_corr_matrix)) + plt.close(fig) + del fig, ax + + # DM of Nuisance Regressors (all) + tr_label = 'TR (Volume #)' + fig, ax = plt.subplots(figsize=(curr_sz - 4.1, def_img_size)) + x_scale_html = ((curr_sz - 4.1) / def_img_size) * 890 + reporting.plot_design_matrix(df, ax=ax) + ax.set_title('Nuisance Design Matrix', fontsize=fontsize_title) + ax.set_xticklabels(ax.get_xticklabels(), rotation=60, ha='right', fontsize=fontsize) + ax.set_yticklabels(ax.get_yticklabels(), fontsize=fontsize) + ax.set_ylabel(tr_label, fontsize=fontsize) + plt.tight_layout() + file_design_matrix = 'Design_matrix' + png_append + fig.savefig(pjoin(out_figure_path, file_design_matrix)) + plt.close(fig) + del fig, ax + + # FD timeseries plot + FD = 'FD' + poss_names = ['FramewiseDisplacement', FD, 'framewisedisplacement', 'fd', 'framewise_displacement'] + fd_idx = [df_orig.columns.__contains__(i) for i in poss_names] + if np.sum(fd_idx) > 0: + FD_name = np.ravel(poss_names)[np.where(fd_idx)[0]][0] + if sum(df_orig[FD_name].isnull()) > 0: + df_orig[FD_name] = df_orig[FD_name].fillna(np.mean(df_orig[FD_name])) + y = df_orig[FD_name].as_matrix() + Nremove = [] + sc_idx = [] + for thr_idx, thr in enumerate(FD_thr): + idx = y >= thr + sc_idx.append(copy.deepcopy(idx)) + for iidx in np.where(idx)[0]: + for buffer in sc_range: + curr_idx = iidx + buffer + if curr_idx >= 0 and curr_idx <= len(idx): + sc_idx[thr_idx][curr_idx] = True + Nremove.append(np.sum(sc_idx[thr_idx])) + + Nplots = len(FD_thr) + sns.set(font_scale=1.5) + sns.set_style('ticks') + fig, axes = plt.subplots(Nplots, 1, figsize=(def_img_size * 1.5, def_img_size / 2), squeeze=False) + sns.despine() + bound = .4 + fd_mean = np.mean(y) + for curr in np.arange(0, Nplots): + axes[curr, 0].plot(y) + axes[curr, 0].plot((-bound, Ntrs + bound), FD_thr[curr] * np.ones((1, 2))[0], '--', color='black') + axes[curr, 0].scatter(np.arange(0, Ntrs), y, s=20) + + if Nremove[curr] > 0: + info = scipy.ndimage.measurements.label(sc_idx[curr]) + for cluster in np.arange(1, info[1] + 1): + temp = np.where(info[0] == cluster)[0] + axes[curr, 0].axvspan(temp.min() - bound, temp.max() + bound, alpha=.5, color='red') + + axes[curr, 0].set_ylabel('Framewise Disp. (' + FD + ')') + axes[curr, 0].set_title(FD + ': ' + str(100 * Nremove[curr] / Ntrs)[0:4] + + '% of scan (' + str(Nremove[curr]) + ' volumes) would be scrubbed (FD thr.= ' + + str(FD_thr[curr]) + ')') + plt.text(Ntrs + 1, FD_thr[curr] - .01, FD + ' = ' + str(FD_thr[curr]), fontsize=fontsize) + plt.text(Ntrs, fd_mean - .01, 'avg = ' + str(fd_mean), fontsize=fontsize) + axes[curr, 0].set_xlim((-bound, Ntrs + 8)) + + plt.tight_layout() + axes[curr, 0].set_xlabel(tr_label) + file_fd_plot = FD + '_timeseries' + png_append + fig.savefig(pjoin(out_figure_path, file_fd_plot)) + plt.close(fig) + del fig, axes + print(FD + ' timeseries plot saved') + + else: + print(FD + ' not found: ' + FD + ' timeseries not plotted') + file_fd_plot = None + + # Carpet and DVARS plots - before & after nuisance regression + + # need to create mask file to input to DVARS function + mask_file = pjoin(out_figure_path, 'mask_temp.nii.gz') + nifti_masker = NiftiMasker(mask_strategy='epi', standardize=False) + nifti_masker.fit(img) + nifti_masker.mask_img_.to_filename(mask_file) + + # create 2 or 3 carpet plots, depending on if LP filter is also applied + Ncarpet = 2 + total_sz = int(16) + carpet_scale = 840 + y_labels = ['Input (voxels)', 'Output \'cleaned\''] + imgs = [img, new_img] + img_files = [img_file, save_img_file] + color = ['red', 'salmon'] + labels = ['input', 'cleaned'] + if lp_filter: + Ncarpet = 3 + total_sz = int(20) + carpet_scale = carpet_scale * (9/8) + y_labels = ['Input', 'Clean Pre-LP', 'Clean LP'] + imgs.insert(1, temp_img) + img_files.insert(1, temp_img_file) + color.insert(1, 'firebrick') + labels.insert(1, 'clean pre-LP') + labels[-1] = 'clean LP' + + dvars = [] + print('Computing dvars...') + for in_file in img_files: + temp = nac.compute_dvars(in_file=in_file, in_mask=mask_file)[1] + dvars.append(np.hstack((temp.mean(), temp))) + del temp + + small_sz = 2 + fig = plt.figure(figsize=(def_img_size * 1.5, def_img_size + ((Ncarpet - 2) * 1))) + row_used = 0 + if np.sum(fd_idx) > 0: # if FD data is available + row_used = row_used + small_sz + ax0 = plt.subplot2grid((total_sz, 1), (0, 0), rowspan=small_sz) + ax0.plot(y) + ax0.scatter(np.arange(0, Ntrs), y, s=10) + curr = 0 + if Nremove[curr] > 0: + info = scipy.ndimage.measurements.label(sc_idx[curr]) + for cluster in np.arange(1, info[1] + 1): + temp = np.where(info[0] == cluster)[0] + ax0.axvspan(temp.min() - bound, temp.max() + bound, alpha=.5, color='red') + ax0.set_ylabel(FD) + + for side in ["top", "right", "bottom"]: + ax0.spines[side].set_color('none') + ax0.spines[side].set_visible(False) + + ax0.set_xticks([]) + ax0.set_xlim((-.5, Ntrs - .5)) + ax0.spines["left"].set_position(('outward', 10)) + + ax_d = plt.subplot2grid((total_sz, 1), (row_used, 0), rowspan=small_sz) + for iplot in np.arange(len(dvars)): + ax_d.plot(dvars[iplot], color=color[iplot], label=labels[iplot]) + ax_d.set_ylabel('DVARS') + for side in ["top", "right", "bottom"]: + ax_d.spines[side].set_color('none') + ax_d.spines[side].set_visible(False) + ax_d.set_xticks([]) + ax_d.set_xlim((-.5, Ntrs - .5)) + ax_d.spines["left"].set_position(('outward', 10)) + ax_d.legend(fontsize=fontsize - 2) + row_used = row_used + small_sz + + st = 0 + carpet_each = int((total_sz - row_used) / Ncarpet) + for idx, img_curr in enumerate(imgs): + ax_curr = plt.subplot2grid((total_sz, 1), (row_used + st, 0), rowspan=carpet_each) + fig = plot_carpet(img_curr, mask_img=nifti_masker.mask_img_, figure=fig, axes=ax_curr) + ax_curr.set_ylabel(y_labels[idx]) + for side in ["bottom", "left"]: + ax_curr.spines[side].set_position(('outward', 10)) + + if idx < len(imgs)-1: + ax_curr.spines["bottom"].set_visible(False) + ax_curr.set_xticklabels('') + ax_curr.set_xlabel('') + st = st + carpet_each + + file_carpet_plot = 'Carpet_plots' + png_append + fig.savefig(pjoin(out_figure_path, file_carpet_plot)) + plt.close() + del fig, ax0, ax_curr, ax_d, dvars + os.remove(mask_file) + print('Carpet/DVARS plots saved') + if lp_filter: + os.remove(temp_img_file) + del temp_img + + # Display T-stat maps for nuisance regressors + # create mean img + img_size = (img.shape[0], img.shape[1], img.shape[2]) + mean_img = nb.Nifti1Image(np.reshape(data_mean, img_size), img.affine) + mx = [] + for idx, col in enumerate(df.columns): + if not 'drift' in col and not constant in col: + con_vector = np.zeros((1, df.shape[1])) + con_vector[0, idx] = 1 + con = results.Tcontrast(con_vector) + mx.append(np.max(np.absolute([con.t.min(), con.t.max()]))) + mx = .8 * np.max(mx) + t_png = 'Tstat_' + file_tstat = [] + for idx, col in enumerate(df.columns): + if not 'drift' in col and not constant in col: + con_vector = np.zeros((1, df.shape[1])) + con_vector[0, idx] = 1 + con = results.Tcontrast(con_vector) + m_img = nb.Nifti1Image(np.reshape(con, img_size), img.affine) + + title_str = col + ' Tstat' + fig = plotting.plot_stat_map(m_img, mean_img, threshold=3, colorbar=True, display_mode='z', vmax=mx, + title=title_str, + cut_coords=7) + file_temp = t_png + col + png_append + fig.savefig(pjoin(out_figure_path, file_temp)) + file_tstat.append({'name': col, 'file': file_temp}) + plt.close() + del fig, file_temp + print(title_str + ' map saved') + + # Display R-sq map for nuisance regressors + m_img = nb.Nifti1Image(np.reshape(rsquare, img_size), img.affine) + title_str = 'Nuisance Rsq' + mx = .95 * rsquare.max() + fig = plotting.plot_stat_map(m_img, mean_img, threshold=.2, colorbar=True, display_mode='z', vmax=mx, + title=title_str, + cut_coords=7) + file_rsq_map = 'Rsquared' + png_append + fig.savefig(pjoin(out_figure_path, file_rsq_map)) + plt.close() + del fig + print(title_str + ' map saved') + + ######### html report + templateLoader = jinja2.FileSystemLoader(searchpath="/") + templateEnv = jinja2.Environment(loader=templateLoader) + + templateVars = {"img_file": img_file, + "save_img_file": save_img_file, + "Ntrs": Ntrs, + "tsv_file": tsv_file, + "col_names": col_names, + "hp_filter": hp_filter, + "lp_filter": lp_filter, + "file_design_matrix": file_design_matrix, + "file_corr_matrix": file_corr_matrix, + "file_fd_plot": file_fd_plot, + "file_rsq_map": file_rsq_map, + "file_tstat": file_tstat, + "x_scale": x_scale_html, + "mtx_scale": mtx_scale, + "file_carpet_plot": file_carpet_plot, + "carpet_scale": carpet_scale + } + + TEMPLATE_FILE = pjoin(os.getcwd(), "denoiser", "report_template.html") + template = templateEnv.get_template(TEMPLATE_FILE) + + outputText = template.render(templateVars) + + html_file = pjoin(out_figure_path, img_name[0:img_name.find('.')] + '.html') + with open(html_file, "w") as f: + f.write(outputText) + + print('') + print('HTML report: ' + html_file) + return new_img + +denoise(img_file, tsv_file, out_path, col_names, hp_filter, lp_filter, out_figure_path)